Demonstration program for using CAS with PPL

[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 >=