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.