Demonstration program for using CAS with PPL

[Michael de Estrada]

The program is probably of little interest to most people on this forum, and is meant to demonstrate usage of CAS functions and operations in PPL programs.

1. The best way I've found to incorporate CAS functions, commands and operations inside a PPL program is to first develop and debug them inside the CAS view and then enclose them in a CAS shell, i.e. <CAS Expression> becomes CAS(" <CAS Expression> ") inside a PPL program.

2. Make sure any variables that are used with a CAS function are global, either by using predefined variables such as A, B, C etc. or by defining them as either global User variables or CAS variables. Never use LOCAL variables or dummy parameters.

3. Avoid using CAS functions if a non-CAS alternate is available. For example, instead of using sqrt(), use ()^.5.

PHP Code:

//Eigenvalue solver for root based on guess "G". The guess should be a positive
//number slightly larger than the lowest non-trivial (zero) root. If it is too low, 
//the result will be a very small number representing zero. Once the lowest root is
//determined, the program can be run again w/o entering the beam data to determine
//additional higher roots if needed, which will be stored in global variable X.
//The natural frequency "F" is proportional to the square of the eigenvalue.
//
EXPORT Beam2_Solve(guess)
BEGIN
//Save angle setting of current active App.
LOCAL Save_AAngle:=AAngle;
//Angle mode must be set to radians for this program to run properly.
AAngle:=1;
G:=guess;
CAS("Z0:=nSolve(DET(m0)=0,x=G)");
//Restore angle setting of current active App. 
AAngle:=Save_AAngle;
X:=ABS(Z0);
//Print eigenvalue.
RETURN ROUND(X,6);
END;

//Program to solve for natural frequency of a beam with flexible end supports
//
// k1 = left support translational spring constant
// k2 = left support rotational spring constant
// k3 = right support translational spring constant
// k4 = right support rotational spring constant
// EI = product of beam material modulus of elasticity and bending moment of inertia
// m = beam mass per unit length
// L = beam length
// guess = guess for eigenvalue solver
//
EXPORT Beam2_Main(k1,k2,k3,k4,EI,m,L,guess)
BEGIN
//Dimensionless coefficients - a negative value designates a rigid support.
A:=k1*L^3/EI;
B:=k2*L/EI;
C:=k3*L^3/EI;
D:=k4*L/EI;
//Create CAS matrix elements.
IF A>=0 THEN
CAS("a11_:=x^3");
CAS("a12_:=-A");
CAS("a13_:=-x^3");
CAS("a14_:=-A");
ELSE
CAS("a11_:=0");
CAS("a12_:=-1");
CAS("a13_:=0");
CAS("a14_:=-1");
END;
IF B>=0 THEN
CAS("a21_:=-B");
CAS("a22_:=-x");
CAS("a23_:=-B");
CAS("a24_:=x");
ELSE
CAS("a21_:=-1");
CAS("a22_:=0");
CAS("a23_:=-1");
CAS("a24_:=0");
END;
IF C>=0 THEN
CAS("a31_:=-x^3*COS(x)-C*SIN(x)");
CAS("a32_:=x^3*SIN(x)-C*COS(x)");
CAS("a33_:=x^3*COSH(x)-C*SINH(x)");
CAS("a34_:=x^3*SINH(x)-C*COSH(x)");
ELSE
CAS("a31_:=-SIN(x)");
CAS("a32_:=-COS(x)");
CAS("a33_:=-SINH(x)");
CAS("a34_:=-COSH(x)");
END;
IF D>=0 THEN
CAS("a41_:=-x*SIN(x)+D*COS(x)");
CAS("a42_:=-x*COS(x)-D*SIN(x)");
CAS("a43_:=x*SINH(x)+D*COSH(x)");
CAS("a44_:=x*COSH(x)+D*SINH(x)");
ELSE
CAS("a41_:=COS(x)");
CAS("a42_:=-SIN(x)");
CAS("a43_:=COSH(x)");
CAS("a44_:=SINH(x)");
END;
//Create CAS eigenvalue matrix.
CAS("m0:=[[a11_,a12_,a13_,a14_],[a21_,a22_,a23_,a24_],
[a31_,a32_,a33_,a34_],[a41_,a42_,a43_,a44_]]");
//Solve for eigenvalue based on guess.
Beam2_Solve(guess);
//Calculate natural frequency in Hz.
F:=(X/L)^2*(EI/m)^.5/(2*pi);
//Print eigenvalue and natural frequency.
RETURN {ROUND(X,6),ROUND(F,3)};
END; 


Revised code to correct misassignment of complex output from nSolve to real global variable X. Replaced X with global complex variable Z0 and then assigned its absolute value (scalar magnitude) to X. Note that the imaginary part of Z0 is approximately zero.

Revised code to use AAngle instead of HAngle to set angle mode in program and then restore AAngle setting in current active App. AAngle has priority over HAngle.

Revised code to add limiting case solution for rigid supports (k-->infinity). A rigid support is flagged by entering the stiffness value k as a negative number. Previously, the stiffness value k had to be entered as a very large positive value.

[Han]

Here's a modification of your code to use only local variables. Moreover, after the program is finished, you won't be left with a lot of temporary CAS variables of the form \( a_{ij\underline{ }} \) active in the CAS view. Lastly, no global variables are altered in the process.

The idea is that since the command will ultimately be delimited by double quotes, then we can use local variables and string addition to do any intermediate calculations which do not actually involve CAS commands. The string manipulation can then be saved as a local variable. So using CAS(locavariable) is fine because the CAS() command itself is a non-CAS command, so it can handle local variables. If localvariable is a string, then all is well.

PHP Code:

EXPORT Beam2_Solve(matrix,guess)
BEGIN
  LOCAL Save_AAngle:=AAngle;
  LOCAL soln;

  // nSolve(DET(matrix)=0,x=guess)
  matrix:="nSolve(DET(" + matrix + ")=0,x=" + guess + ")";

  AAngle:=1; // set radian mode

  // ensure that user's angle mode is always restored, even if error occurs
  IFERR soln:=CAS(matrix); THEN  END;
  AAngle:=Save_AAngle;
  RETURN ROUND(ABS(soln),6);
END;

//Program to solve for natural frequency of a beam with flexible end supports
//
// k1 = left support translational spring constant
// k2 = left support rotational spring constant
// k3 = right support translational spring constant
// k4 = right support rotational spring constant
// EI = product of beam material modulus of elasticity and bending moment of inertia
// m = beam mass per unit length
// L = beam length
// guess = guess for eigenvalue solver
//
EXPORT Beam2_Main(k1,k2,k3,k4,EI,m,L,guess)
BEGIN

  LOCAL a,b,c,d,s1,s2,s;
  LOCAL matname="m0";

  // initialize matrix rows to form that does not involve pre-computed values
  LOCAL row1="[[0,-1,0,-1],";
  LOCAL row2="[-1,0,-1,0],";
  LOCAL row3="[-SIN(x),-COS(x),-SINH(x),-COSH(x)],";
  LOCAL row4="[COS(x),-SIN(x),COSH(x),SINH(x)]]";

  a:=k1*L^3/EI;  b:=k2*L/EI;
  c:=k3*L^3/EI;  d:=k4*L/EI;

  // rebuild rows only as necessary
  IF a>=0 THEN
    row1:="[[x^3,-" + a + ",-x^3,-" + a + "],";
  END;
  
  IF b>=0 THEN
    row2:="[-" + b + ",-x,-" + b + ",x],";
  END;

  IF c>=0 THEN
    row3:="[-x^3*COS(x)-" + c +  "*SIN(x),x^3*SIN(x)-" + c + "*COS(x),x^3*COSH(x)-" + c + "*SINH(x),x^3*SINH(x)-" + c + "*COSH(x)],";
  END;

  IF d>=0 THEN
    row4:="[-x*SIN(x)+" + d + "*COS(x),-x*COS(x)-" + d + "*SIN(x),x*SINH(x)+" + d + "*COSH(x),x*COSH(x)+" + d + "*SINH(x)]]";
  END;

  // ask user to save matrix
  s:=INPUT(matname,"Save Matrix?","Name=","Enter CAS name for matrix",matname);
  if s AND TYPE(matname)==2 then
    matname:=matname+":="+row1+row2+row3+row4;
    iferr CAS(matname); then msgbox("Unable to save matrix."); end;
  end;

  // send subroutine both the string to pass to CAS() and guess
  s1:=Beam2_Solve(row1+row2+row3+row4,guess);
  s2:=(s1/L)^2*(EI/m)^.5/(2*pi);
  RETURN {ROUND(s1,6),ROUND(s2,3)};
END; 


[Michael de Estrada]

Your modification has one serious drawback; namely you can't run Beam2_Solve directly, since the parameter matrix has to be evaluated each time the program is run. The reason I wrote the program that way and exported Beam2_Solve was to permit initial entry of the matrix parameters through Beam2_Main to obtain an initial root, and then run Beam2_Solve repeatedly to obtain additional roots w/o having to re-enter all the beam parameters, which don't change. This is why I chose to make the matrix a global CAS variable. I agree, however, that it wasn't necessary to make the individual matrix coefficients into permanent CAS variables.

[Han]

Oh, I misunderstood the use of m0, then. I thought it was just a temporary result, so I took it to the extreme and made everything local variables. :-)

I adjusted the code to enable users to save the matrix as any CAS variable in the main program, or they can choose to also not save the result (by pressing [Esc]). Beam2_Solve() will still work whether or not the parameter 'matrix' is a string. I suppose it still has a minor drawback in that now you have to also specify which matrix to use when doing repeated calls.

[Michael de Estrada]

(01-21-2014 05:05 AM)Han Wrote:  I adjusted the code to enable users to save the matrix as any CAS variable in the main program, or they can choose to also not save the result (by pressing [Esc]). Beam2_Solve() will still work whether or not the parameter 'matrix' is a string. I suppose it still has a minor drawback in that now you have to also specify which matrix to use when doing repeated calls.

When running Beam2_Main, I get the error "bad argument type" after entering the matrix name from the command line. The program halts, and the results are not output. The matrix does get created, so the problem occurs someplace after that point. The matrix name is entered with double quotes, and if I remove them the error becomes "invalid dimension." Also, you cannot enter the matrix name from the program console, because it won't accept a string and I get "invalid input." It does work, however, if I run Beam2_Solve from the command line and enter the existing matrix name with double quotes.

Edit: I found the problem. You added the string variable "m" for the matrix, which conflicted with the mass dummy parameter "m". When computing the variable s2, this caused the bad argument type error. Changing the dummy parameter "m" to "mass" resolves this problem.

[Han]

Quote:Edit: I found the problem. You added the string variable "m" for the matrix, which conflicted with the mass dummy parameter "m". When computing the variable s2, this caused the bad argument type error. Changing the dummy parameter "m" to "mass" resolves this problem.

My apologies -- I went ahead and changed the source to use 'matname'

[Michael de Estrada]

(01-21-2014 04:38 PM)Han Wrote:  

Quote:Edit: I found the problem. You added the string variable "m" for the matrix, which conflicted with the mass dummy parameter "m". When computing the variable s2, this caused the bad argument type error. Changing the dummy parameter "m" to "mass" resolves this problem.

My apologies -- I went ahead and changed the source to use 'matname'

No problem. I think you forgot to change m to matname in a few places.

[Michael de Estrada]

I also revised my original code to replace the 16 individual matrix coefficients aij_ with four row variables row1 through row4, which reduces the number of CAS variables from 16 to 4 and simplifies the code:

PHP Code:

//Eigenvalue solver for root based on guess "G". The guess should be a positive
//number slightly larger than the lowest non-trivial (zero) root. If it is too low, 
//the result will be a very small number representing zero. Once the lowest root is
//determined, the program can be run again w/o entering the beam data to determine
//additional higher roots if needed, which will be stored in global variable X.
//The natural frequency "F" is proportional to the square of the eigenvalue.

EXPORT Beam3_Solve(guess)
BEGIN

//Save angle setting of current active App.

LOCAL Save_AAngle:=AAngle;

//Angle mode must be set to radians for this program to run properly.

AAngle:=1;

G:=guess;

CAS("Z0:=nSolve(DET(matev)=0,x=G)");

//Restore angle setting of current active App.
 
AAngle:=Save_AAngle;

X:=ABS(Z0);

//Print eigenvalue.

RETURN ROUND(X,6);
END;

//Program to solve for natural frequency of a beam with flexible end supports

// k1 = left support translational spring constant
// k2 = left support rotational spring constant
// k3 = right support translational spring constant
// k4 = right support rotational spring constant
// EI = product of beam material modulus of elasticity and bending moment of inertia
// m = beam mass per unit length
// L = beam length
// guess = guess for eigenvalue solver

EXPORT Beam3_Main(k1,k2,k3,k4,EI,m,L,guess)
BEGIN

//Dimensionless coefficients - a negative value designates a rigid support.

A:=k1*L^3/EI;
B:=k2*L/EI;
C:=k3*L^3/EI;
D:=k4*L/EI;

//Create CAS matrix rows.

IF A>=0 THEN CAS("row1:=[x^3,-A,-x^3,-A]");
ELSE CAS("row1:=[0,-1,0,-1]");
END;

IF B>=0 THEN CAS("row2:=[-B,-x,-B,x]");
ELSE CAS("row2:=[-1,0,-1,0]");
END;

IF C>=0 THEN CAS("row3:=[-x^3*COS(x)-C*SIN(x),x^3*SIN(x)-C*COS(x),
  x^3*COSH(x)-C*SINH(x),x^3*SINH(x)-C*COSH(x)]");
ELSE CAS("row3:=[-SIN(x),-COS(x),-SINH(x),-COSH(x)]");
END;

IF D>=0 THEN CAS("row4:=[-x*SIN(x)+D*COS(x),-x*COS(x)-D*SIN(x),x*SINH(x)+D*COSH(x),
  x*COSH(x)+D*SINH(x)]");
ELSE CAS("row4:=[COS(x),-SIN(x),COSH(x),SINH(x)]");
END;

//Create CAS eigenvalue matrix.

CAS("matev:=[row1,row2,row3,row4]");

//Solve for eigenvalue based on guess.

Beam3_Solve(guess);

//Calculate natural frequency in Hz.

F:=(X/L)^2*(EI/m)^.5/(2*pi);

//Print eigenvalue and natural frequency.

RETURN {ROUND(X,6),ROUND(F,3)};
END; 


Ironically, this is the way I formulated my original RPL solution:

EQ <<COEF 'MAT' {1,1} ROW1 ROW2 ROW3 ROW4 MAT DET ROT ROT DROP2>>

where the coefficients A,B,C,D are created by:

COEF <<LAM SIN 'SSL' STO LAM COS 'CCL' STO LAM SINH 'SHL' STO
LAM COSH 'CHL' STO LAM 3 ^ 'L3' STO>>

and I used user flags 1 through 4 to identify the rigid support cases, for example:

ROW1 << IF 1 FC? THEN L3 NEG PUTI K1 PUTI L3 PUTI K1 PUTI ELSE 0 PUTI 1 PUTI 0 PUTI 1 PUTI END >>

I'm changing the RPL code to use negative values instead of flags to designate rigid supports e.g. 1 FC? becomes K1 0 >=

[Han]

(01-21-2014 04:42 PM)Michael de Estrada Wrote:  

(01-21-2014 04:38 PM)Han Wrote:  My apologies -- I went ahead and changed the source to use 'matname'

No problem. I think you forgot to change m to matname in a few places.

Indeed my sloppiness knows no bounds ;_;