> restart;
> # PIECEWISE LINEAR RAYLEIGH-RITZ ALGORITHM 11.5
> #
> # To approximate the solution of the boundary-value problem
> #
> #      -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1,
> #             Y(0) = Y(1) = 0,
> #
> # with a piecewise linear function:
> #
> # INPUT:   integer N; mesh points X(0) = 0 < X(1) < ...
> #          < X(N) < X(N+1) = 1
> #
> # OUTPUT:  coefficients C(1),...,C(N) of the basis functions
> alg115 := proc() local SIMPSON, AA, OK, N, X, FLAG, HC, J, H, NAME, INP, N1, Q, ALPHA, BETA, B, A, ZETA, Z, C, J1, OUP; global  F, QQ, P;
> SIMPSON := proc(FN,A,B) local H, I, Y, Z, simpson;
> H := (B-A)/4;
> for I from 0 to 4 do
> Y := A+I*H;
> if FN = 1 then
> Z[I] := (4-I)*I*H*H*QQ(Y);
> fi;
> if FN = 2 then
> Z[I] := (I*H)*(I*H)*QQ(Y);
> fi;
> if FN = 3 then
> Z[I] := (H*(4-I))*(H*(4-I))*QQ(Y);
> fi;
> if FN = 4 then
> Z[I] := P(Y);
> fi;
> if FN = 5 then
> Z[I] := I*H*F(Y);
> fi;
> if FN = 6 then
> Z[I] := (4-I)*H*F(Y);
> fi;
> od;
> simpson := (Z[0]+Z[4]+2*Z[2]+4*(Z[1]+Z[3]))*H/3;
> RETURN(simpson);
> end; 
> printf(`This is the Piecewise Linear Rayleigh-Ritz Method.\n`);
> printf(`Input F(X), Q(X), and P(X) in terms of x separated by a 
> space.\n`);
> printf(`For example:\n`);   
> printf(`2*3.141592654^2*sin(3.141592654*x) 3.141592654^2 1\n`);
> F := scanf(`%a`)[1];
> QQ := scanf(`%a`)[1];
> P := scanf(`%a`)[1];
> F := unapply(F,x);
> QQ := unapply(QQ,y);
> P := unapply(P,z);
> printf(`X(0), ..., X(N+1) are to be supplied.\n`);
> printf(`Are the preparations complete? Answer Y or N.\n`);
> AA := scanf(`\n%c`)[1];
> OK := FALSE;
> if AA = "Y" or AA = "y" then
> OK := FALSE;
> while OK = FALSE do
> printf(`Input integer N where X(0) = 0, X(N+1) = 1.\n`);
> N := scanf(`%d`)[1];
> if N <= 1 then
> printf(`N must be greater than one.\n`);
> else
> OK := TRUE;
> fi;
> od;
> X[0] := 0;
> X[N+1] := 1;
> printf(`Choice of method to input X(1), ..., X(N):\n`);
> printf(`1.  Input from keyboard at the prompt\n`);
> printf(`2.  Equally spaced nodes to be calculated\n`);
> printf(`3.  Input from text file\n`);
> printf(`Please enter 1, 2, or 3.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> HC := 1/(N+1);
> for J from 1 to N do
> X[J] := evalf(J*HC);
> H[J-1] := HC;
> od;
> H[N] := HC;
> else
> if FLAG = 3 then
> printf(`Has the input file been created?  `);
> printf(`Enter Y or N.\n`);
> AA := scanf(`\n%c`)[1];
> if AA = "Y" or AA = "y" then
> printf(`Enter the input file name using the format\n`);
> printf(` - drive:\\name.ext,\n`);
> printf(`for example:  A:\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> INP := fopen(NAME,READ,TEXT);
> for J from 1 to N do
> X[J] := evalf(fscanf(INP, `%f`)[1]);
> od;
> for J from 0 to N do
> H[J] := X[J+1]-X[J];
> od;
> fclose(INP);
> else
> printf(`The program will end so that the input `);
> printf(`file can be created.\n`);
> OK := FALSE;
> fi;
> else
> for J from 1 to N do
> printf(`Input X(%d).\n`, J);
> X[J] := evalf(scanf(`%f`)[1]);
> H[J-1] := X[J]-X[J-1];
> od;
> H[N] := X[N+1]-X[N];
> fi;
> fi;
> else
> printf(`The program will end so that the preparations\n`);
> printf(`can be completed.\n`);
> OK := FALSE;
> fi;
> # Step 1 is done in the input procedure.
> if OK = TRUE then
> N1 := N-1;
> # Step 3
> for J from 1 to N1 do
> Q[0,J-1] := SIMPSON(1,X[J],X[J+1])/((H[J])*(H[J]));
> Q[1,J-1] := SIMPSON(2,X[J-1],X[J])/((H[J-1])*(H[J-1]));
> Q[2,J-1] := SIMPSON(3,X[J],X[J+1])/((H[J])*(H[J]));
> Q[3,J-1] := SIMPSON(4,X[J-1],X[J])/((H[J-1])*(H[J-1]));
> Q[4,J-1] := SIMPSON(5,X[J-1],X[J])/H[J-1] ;
> Q[5,J-1] := SIMPSON(6,X[J],X[J+1])/H[J];
> od;
> Q[1,N-1] := SIMPSON(2,X[N-1],X[N])/((H[N-1])*(H[N-1]));
> Q[2,N-1] := SIMPSON(3,X[N],X[N+1])/((H[N])*(H[N]));
> Q[3,N-1] := SIMPSON(4,X[N-1],X[N])/((H[N-1])*(H[N-1]));
> Q[3,N] := SIMPSON(4,X[N],X[N+1])/((H[N])*(H[N]));
> Q[4,N-1] := SIMPSON(5,X[N-1],X[N])/H[N-1];
> Q[5,N-1] := SIMPSON(6,X[N],X[N+1])/H[N];
> # Step 4
> for J from 1 to N1 do
> ALPHA[J-1] := Q[3,J-1]+Q[3,J]+Q[1,J-1]+Q[2,J-1];
> BETA[J-1] := Q[0,J-1]-Q[3,J];
> B[J-1] := Q[4,J-1]+Q[5,J-1];
> od;
> # Step 5
> ALPHA[N-1] := Q[3,N-1]+Q[3,N]+Q[1,N-1]+Q[2,N-1];
> B[N-1] := Q[4,N-1]+Q[5,N-1];
> # Steps 6 - 10 solve a tridiagonal linear system using Algorithm 6.7
> # Step 6
> A[0] := ALPHA[0];
> ZETA[0] := BETA[0]/ALPHA[0];
> Z[0] := B[0]/A[0];
> # Step 7
> for J from 2 to N1 do
> A[J-1] := ALPHA[J-1]-BETA[J-2]*ZETA[J-2];
> ZETA[J-1] := BETA[J-1]/A[J-1];
> Z[J-1] := (B[J-1]-BETA[J-2]*Z[J-2])/A[J-1];
> od;
> # Step 8
> A[N-1] := ALPHA[N-1] - BETA[N-2] * ZETA[N-2];
> Z[N-1] := (B[N-1]-BETA[N-2]*Z[N-2])/A[N-1];
> # Step 9
> C[N-1] := evalf(Z[N-1]);
> for J from 1 to N1 do
> J1 := N - J;
> C[J1-1] := evalf(Z[J1-1]-ZETA[J1-1]*C[J1]);
> od;
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example  A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `PIECEWISE LINEAR RAYLEIGH-RITZ METHOD\n\n`);
> fprintf(OUP, ` I      X(I-1)       X(I)       X(I+1)        C(I)\n\n`);
> for J from 1 to N do
> fprintf(OUP, `%3d %11.8f %11.8f %11.8f %13.8f\n`, J, X[J-1],X[J], X[J+1], C[J-1]);
> od;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg115();