> restart;
> Digits := 20;
> # CUBIC SPLINE RAYLEIGH-RITZ ALGORITHM 11.6
> #
> # To approximate the solution to the boundary-value problem
> #
> #    -D(P(X)Y')/DX + Q(X)Y = F(X), 0 <= X <= 1, Y(0)=Y(1)=0
> #
> # with a sum of cubic splines:
> #
> # INPUT:  Integer n
> #
> # OUTPUT: Coefficients C(0),...,C(n+1) of the basis functions
> alg116 := proc() global F,Q,P;
> INTE := proc(J,JJ) local inte;
> inte := JJ-J+3;
> RETURN(inte);
> end;
> XINT := proc(XU,XL,A1,B1,C1,D1,A2,B2,C2,D2,A3,B3,C3,D3)
> AA := A1*A2;
> BB := A1*B2+A2*B1;
> CC := A1*C2+B1*B2+C1*A2;
> DD := A1*D2+B1*C2+C1*B2+D1*A2;
> EE := B1*D2+C1*C2+D1*B2;
> FF := C1*D2+D1*C2;
> GG := D1*D2;
> C[9] := AA*A3;
> C[8] := (AA*B3+BB*A3)/2;
> C[7] := (AA*C3+BB*B3+CC*A3)/3;
> C[6] := (AA*D3+BB*C3+CC*B3+DD*A3)/4;
> C[5] := (BB*D3+CC*C3+DD*B3+EE*A3)/5;
> C[4] := (CC*D3+DD*C3+EE*B3+FF*A3)/6;
> C[3] := (DD*D3+EE*C3+FF*B3+GG*A3)/7;
> C[2] := (EE*D3+FF*C3+GG*B3)/8;
> C[1] := (FF*D3+GG*C3)/9;
> C[0] := (GG*D3)/10;
> XHIGH := 0;
> XLOW := 0;
> for I1 from 1 to 10 do
> XHIGH := (XHIGH+C[I1-1])*XU;
> XLOW := (XLOW+C[I1-1])*XL;
> od;
> xint := XHIGH-XLOW;
> RETURN(xint);
> end;
> printf(`This is the Cubic Spline Rayleigh-Ritz Method.\n`);
> OK := FALSE;
> printf(`Input functions P(X), Q(X), and F(X) in terms of x 
> separated by a space.\n`);
> printf(`For example: \n`);
> printf(`1 3.141592654^2 2*3.141592654^2*sin(3.1415926548*x)\n`);
> printf(`separated by at least one space.\n`);
> Pf := scanf(`%a`)[1];
> Qf := scanf(`%a`)[1];
> Ff := scanf(`%a`)[1];
> FPL := evalf(subs(x=0,diff(Ff,x)));
> FPR := evalf(subs(x=1,diff(Ff,x)));
> QPL := evalf(subs(x=0,diff(Qf,x)));
> QPR := evalf(subs(x=1,diff(Qf,x)));
> PPL := evalf(subs(x=0,diff(Pf,x)));
> PPR := evalf(subs(x=1,diff(Pf,x)));
> F := unapply(Ff,x,y,z);
> Q := unapply(Qf,x,y,z);
> P := unapply(Pf,x,y,z);
> while OK = FALSE do
> printf(`Input a positive integer n, where x(0) = 0, `);
> printf(`..., x(n+1) = 1.\n`);
> N := scanf(`%d`)[1];
> if N <= 0 then
> printf(`Number must be a positive integer.\n`);
> else
> OK := TRUE;
> fi;
> od;
> if OK = TRUE then
> 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, `CUBIC SPLINE RAYLEIGH-RITZ METHOD\n\n`);
> # Step 1
> H := 1/(N+1);
> N1 := N+1;
> N2 := N+2;
> N3 := N+3;
> # Initialize matrix A at 0, note that A[I,N+3] = B[I]
> for I from 1 to N2 do
> for J from 1 to N3 do
> A[I-1,J-1] := 0;
> od;
> od;
> # Step 2
> # X[1] = 0, ... , X[I] = (I-1)*H, ... , X[N+1] = 1-H, X[N+2] = 1
> for I from 1 to N2 do
> X[I-1] := (I-1)*H;
> od;
> # STEPS 3 and 4 are implemented in what follows. Initialize coefficients
> # CO[I,J,K], DCO[I,J,K] 
> for I from 1 to N2 do
> for J from 1 to 4 do
> # JJ corresponds the coefficients of phi and phi' to the proper interval 
> # involving J 
> JJ := I+J-3;
> CO[I-1,J-1,0] := 0;
> CO[I-1,J-1,1] := 0;
> CO[I-1,J-1,2] := 0;
> CO[I-1,J-1,3] := 0;
> E := I-1;
> OK := TRUE;
> if JJ < I-2 or JJ >= I+2 then
> OK := FALSE;
> fi;
> if I = 1 and JJ < I then
> OK := FALSE;
> fi;
> if I = 2 and JJ < I-1 then
> OK := FALSE;
> fi;
> if I = N+1 and JJ > N+1 then
> OK := FALSE;
> fi;
> if I = N+2 and JJ >= N+2 then
> OK := FALSE;
> fi;
> if OK = TRUE then
> if JJ <= I-2 then
> CO[I-1,J-1,0] := (((-E+6)*E-12)*E+8)/24;
> CO[I-1,J-1,1] := ((E-4)*E+4)/(8*H);
> CO[I-1,J-1,2] := (-E+2)/(8*H^2);
> CO[I-1,J-1,3] := 1/(24*H^3);
> OK := FALSE;
> else
> if JJ > I then
> CO[I-1,J-1,0] := (((E+6)*E+12)*E+8)/24;
> CO[I-1,J-1,1] := ((-E-4)*E-4)/(8*H);
> CO[I-1,J-1,2] := (E+2)/(8*H^2);
> CO[I-1,J-1,3] := -1/(24*H^3);
> OK := FALSE;
> else
> if JJ > I-1 then 
> CO[I-1,J-1,0] := ((-3*E-6)*E*E+4)/24;
> CO[I-1,J-1,1] := (3*E+4)*E/(8*H);
> CO[I-1,J-1,2] := (-3*E-2)/(8*H^2);
> CO[I-1,J-1,3] := 1/(8*H^3);
> if I <> 1 and I <> N+1 then
> OK := FALSE;
> fi;
> else
> CO[I-1,J-1,0] := ((3*E-6)*E*E+4)/24;
> CO[I-1,J-1,1] := (-3*E+4)*E/(8*H);
> CO[I-1,J-1,2] := (3*E-2)/(8*H^2);
> CO[I-1,J-1,3] := -1/(8*H^3);
> if I <> 2 and I <> N+2 then
> OK := FALSE;
> fi;
> fi;
> fi;
> fi;
> fi;
> if OK = TRUE then
> if I <= 2 then
> AA := 1/24;
> BB := -1/(8*H);
> CC := 1/(8*H^2);
> DD := -1/(24*H^3);
> if I = 2 then
> CO[I-1,J-1,0] := CO[I-1,J-1,0]-AA;
> CO[I-1,J-1,1] := CO[I-1,J-1,1]-BB;
> CO[I-1,J-1,2] := CO[I-1,J-1,2]-CC;
> CO[I-1,J-1,3] := CO[I-1,J-1,3]-DD;
> else
> CO[I-1,J-1,0] := CO[I-1,J-1,0]-4*AA;
> CO[I-1,J-1,1] := CO[I-1,J-1,1]-4*BB;
> CO[I-1,J-1,2] := CO[I-1,J-1,2]-4*CC;
> CO[I-1,J-1,3] := CO[I-1,J-1,3]-4*DD;
> fi;
> else
> EE := N+2;
> AA := (((-EE+6)*EE-12)*EE+8)/24;
> BB := ((EE-4)*EE+4)/(8*H);
> CC := (-EE+2)/(8*H^2);
> DD := 1/(24*H^3);
> if I = N+1 then
> CO[I-1,J-1,0] := CO[I-1,J-1,0]-AA;
> CO[I-1,J-1,1] := CO[I-1,J-1,1]-BB;
> CO[I-1,J-1,2] := CO[I-1,J-1,2]-CC;
> CO[I-1,J-1,3] := CO[I-1,J-1,3]-DD;
> else
> CO[I-1,J-1,0] := CO[I-1,J-1,0]-4*AA;
> CO[I-1,J-1,1] := CO[I-1,J-1,1]-4*BB;
> CO[I-1,J-1,2] := CO[I-1,J-1,2]-4*CC;
> CO[I-1,J-1,3] := CO[I-1,J-1,3]-4*DD;
> fi;
> fi;
> fi;
> DCO[I-1,J-1,0] := 0;
> DCO[I-1,J-1,1] := 0;
> DCO[I-1,J-1,2] := 0;
> E := I-1;
> OK := TRUE;
> if JJ < I-2 or JJ >= I+2 then
> OK := FALSE;
> fi;
> if I = 1 and JJ < I then
> OK := FALSE;
> fi;
> if I = 2 and JJ < I-1 then
> OK := FALSE;
> fi;
> if I = N+1 and JJ > N+1 then
> OK := FALSE;
> fi;
> if I = N+2 and JJ >= N+2 then
> OK := FALSE;
> fi;
> if OK = TRUE then
> if JJ <= I-2 then
> DCO[I-1,J-1,0] := ((E-4)*E+4)/(8*H);
> DCO[I-1,J-1,1] := (-E+2)/(4*H^2);
> DCO[I-1,J-1,2] := 1/(8*H^3);
> OK := FALSE;
> else
> if JJ > I then
> DCO[I-1,J-1,0] := ((-E-4)*E-4)/(8*H);
> DCO[I-1,J-1,1] := (E+2)/(4*H^2);
> DCO[I-1,J-1,2] := -1/(8*H^3);
> OK := FALSE;
> else
> if JJ > I-1 then
> DCO[I-1,J-1,0] := (3*E+4)*E/(8*H);
> DCO[I-1,J-1,1] := (-3.0*E-2.0)/(4.0*H^2);
> DCO[I-1,J-1,2] := 3/(8*H^3);
> if I <> 1 and I <> N+1 then
> OK := FALSE;
> fi;
> else
> DCO[I-1,J-1,0] := (-3*E+4)*E/(8*H);
> DCO[I-1,J-1,1] := (3*E-2)/(4*H^2);
> DCO[I-1,J-1,2] := -3/(8*H^3);
> if I <> 2 and I <> N+2 then
> OK := FALSE;
> fi;
> fi;
> fi;
> fi;
> fi;
> if OK = TRUE then
> if I <= 2 then
> AA := -1/(8*H);
> BB := 1/(4*H^2);
> CC := -1/(8*H^3);
> if I = 2 then
> DCO[I-1,J-1,0] := DCO[I-1,J-1,0]-AA;
> DCO[I-1,J-1,1] := DCO[I-1,J-1,1]-BB;
> DCO[I-1,J-1,2] := DCO[I-1,J-1,2]-CC;
> else
> DCO[I-1,J-1,0] := DCO[I-1,J-1,0]-4*AA;
> DCO[I-1,J-1,1] := DCO[I-1,J-1,1]-4*BB;
> DCO[I-1,J-1,2] := DCO[I-1,J-1,2]-4*CC;
> fi;
> else
> EE := N+2;
> AA := ((EE-4)*EE+4)/(8*H);
> BB := (-EE+2)/(4*H^2);
> CC := 1/(8*H^3);
> if I = N+1 then
> DCO[I-1,J-1,0] := DCO[I-1,J-1,0]-AA;
> DCO[I-1,J-1,1] := DCO[I-1,J-1,1]-BB;
> DCO[I-1,J-1,2] := DCO[I-1,J-1,2]-CC;
> else
> DCO[I-1,J-1,0] := DCO[I-1,J-1,0]-4*AA;
> DCO[I-1,J-1,1] := DCO[I-1,J-1,1]-4*BB;
> DCO[I-1,J-1,2] := DCO[I-1,J-1,2]-4*CC;
> fi;
> fi;
> fi;
> od;
> od;
> # Output the basis functions. 
> fprintf(OUP, `Basis Function: A + B*X + C*X**2 + D*X**3\n\n`);
> fprintf(OUP, `  A  B  C  D\n\n`);
> for I from 1 to N2 do
> fprintf(OUP, `phi( %d )\n\n`, I);
> for J from 1 to 4 do
> if I <> 1 or (J <> 1 and J <> 2) then
> if I <> 2 or J <> 2 then
> if I <> N1 or J <> 4 then
> if I <> N2 or (J <> 3 and J <> 4) then
> JJ1 := I+J-3;
> JJ2 := I+J-2;
> fprintf(OUP, `On (X( %d ), X( %d )) `, JJ1, JJ2);
> for K from 1 to 4 do
> fprintf(OUP, ` %12.8f `, CO[I-1,J-1,K-1]);
> od;
> fprintf(OUP, `\n`);
> fi;
> fi;
> fi;
> fi;
> od;
> od;
> # Obtain coefficients for F, P, Q  
> for I from 1 to N2 do
> AA[I-1] := evalf(F(X[I-1]));
> od;
> XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*FPL;
> XA[N2-1] := 3.0*FPR-3.0*(AA[N2-1]-AA[N2-2])/H;
> XL[0] := 2.0*H;
> XU[0] := 0.5;
> XZ[0] := XA[0]/XL[0];
> for I from 2 to N1 do
> XA[I-1] := 3.0*(AA[I]-2.0*AA[I-1]+AA[I-2])/H;
> XL[I-1] := H*(4.0-XU[I-2]);
> XU[I-1] := H/XL[I-1];
> XZ[I-1] := (XA[I-1]-H*XZ[I-2])/XL[I-1];
> od;
> XL[N2-1] := H*(2.0-XU[N2-2]);
> XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1];
> CC[N2-1] := XZ[N2-1];
> for I from 1 to N1 do
> J := N2-I;
> CC[J-1] := XZ[J-1]-XU[J-1]*CC[J];
> BB[J-1] := (AA[J]-AA[J-1])/H-H*(CC[J]+2.0*CC[J-1])/3.0;
> DD[J-1] := (CC[J]-CC[J-1])/(3.0*H);
> od;
> for I from 1 to N1 do
> AF[I-1] := ((-DD[I-1]*X[I-1]+CC[I-1])*X[I-1]-BB[I-1])*X[I-1]+AA[I-1];
> BF[I-1] := (3.0*DD[I-1]*X[I-1]-2.0*CC[I-1])*X[I-1]+BB[I-1];
> CF[I-1] := CC[I-1]-3.0*DD[I-1]*X[I-1];
> DF[I-1] := DD[I-1];
> od;
> for I from 1 to N2 do
> AA[I-1] := evalf(P(X[I-1]));
> od;
> XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*PPL;
> XA[N2-1] := 3.0*PPR-3.0*(AA[N2-1]-AA[N2-2])/H;
> XL[0] := 2.0*H;
> XU[0] := 0.5;
> XZ[0] := XA[0]/XL[0];
> for I from 2 to N1 do
> XA[I-1] := 3.0*(AA[I]-2.0*AA[I-1]+AA[I-2])/H;
> XL[I-1] := H*(4.0-XU[I-2]);
> XU[I-1] := H/XL[I-1];
> XZ[I-1] := (XA[I-1]-H*XZ[I-2])/XL[I-1];
> od;
> XL[N2-1] := H*(2.0-XU[N2-2]);
> XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1];
> CC[N2-1] := XZ[N2-1];
> for I from 1 to N1 do
> J := N2-I;
> CC[J-1] := XZ[J-1]-XU[J-1]*CC[J];
> BB[J-1] := (AA[J]-AA[J-1])/H -H*(CC[J]+2.0*CC[J-1])/3.0;
> DD[J-1] := (CC[J]-CC[J-1])/(3.0*H);
> od;
> for I from 1 to N1 do
> AP[I-1] := ((-DD[I-1]*X[I-1]+CC[I-1])*X[I-1]-BB[I-1])*X[I-1]+AA[I-1];
> BP[I-1] := (3.0*DD[I-1]*X[I-1]-2.0*CC[I-1])*X[I-1]+BB[I-1];
> CP[I-1] := CC[I-1]-3.0*DD[I-1]*X[I-1];
> DP[I-1] := DD[I-1];
> od;
> for I from 1 to N2 do
> AA[I-1] := evalf(Q(X[I-1]));
> od;
> XA[0] := 3.0*(AA[1]-AA[0])/H-3.0*QPL;
> XA[N2-1] := 3.0*QPR-3.0*(AA[N2-1]-AA[N2-2])/H;
> XL[0] := 2.0*H;
> XU[0] := 0.5;
> XZ[0] := XA[0]/XL[0];
> for I from 2 to N1 do
> XA[I-1] := 3.0*(AA[I]-2.0*AA[I-1]+AA[I-2])/H;
> XL[I-1] := H*(4.0-XU[I-2]);
> XU[I-1] := H/XL[I-1];
> XZ[I-1] := (XA[I-1]-H*XZ[I-2])/XL[I-1];
> od;
> XL[N2-1] := H*(2.0-XU[N2-2]);
> XZ[N2-1] := (XA[N2-1]-H*XZ[N2-2])/XL[N2-1];
> CC[N2-1] := XZ[N2-1];
> for I from 1 to N1 do
> J := N2-I;
> CC[J-1] := XZ[J-1]-XU[J-1]*CC[J];
> BB[J-1] := (AA[J]-AA[J-1])/H -H*(CC[J]+2.0*CC[J-1])/3.0;
> DD[J-1] := (CC[J]-CC[J-1])/(3.0*H);
> od;
> for I from 1 to N1 do
> AQ[I-1] := ((-DD[I-1]*X[I-1]+CC[I-1])*X[I-1]-BB[I-1])*X[I-1]+AA[I-1];
> BQ[I-1] := (3.0*DD[I-1]*X[I-1]-2.0*CC[I-1])*X[I-1]+BB[I-1];
> CQ[I-1] := CC[I-1]-3.0*DD[I-1]*X[I-1];
> DQ[I-1] := DD[I-1];
> od;
> # Steps 5-9 are implemented in what follows
> for I from 1 to N2 do
> # indices for limits of integration for A[I,J] and B[I] 
> J1 := min(I+2,N+2);
> J0 := max(I-2,1);
> J2 := J1-1;
> # Integrate over each subinterval where phi(I) is nonzero
> for JJ from J0 to J2 do
> # limits of integration for each call
> XU := X[JJ];
> XL := X[JJ-1];
> # coefficients of bases
> K := INTE(I,JJ);
> A1 := DCO[I-1,K-1,0];
> B1 := DCO[I-1,K-1,1];
> C1 := DCO[I-1,K-1,2];
> D1 := 0;
> A2 := CO[I-1,K-1,0];
> B2 := CO[I-1,K-1,1];
> C2 := CO[I-1,K-1,2];
> D2 := CO[I-1,K-1,3];
> # call subprogram for integrations
> A[I-1,I-1] := A[I-1,I-1]+XINT(XU,XL,AP[JJ-1],BP[JJ-1],CP[JJ-1],DP[JJ-1],A1,B1,C1,D1,A1,B1,C1,D1)+XINT(XU,XL,AQ[JJ-1],BQ[JJ-1],CQ[JJ-1],DQ[JJ-1],A2,B2,C2,D2,A2,B2,C2,D2);
> A[I-1,N+2]:=A[I-1,N+2]+XINT(XU,XL,AF[JJ-1],BF[JJ-1],CF[JJ-1],DF[JJ-1],A2,B2,C2,D2,1,0,0,0);
> od;
> # compute A[I,J] for J = I+1, ..., Min(I+3,N+2)
> K3 := I+1;
> if K3 <= N2 then
> K2 := min(I+3,N+2);
> for J from K3 to K2 do
> J0 := max(J-2,1);
> for JJ from J0 to J2 do
> XU := X[JJ];
> XL := X[JJ-1];
> K := INTE(I,JJ);
> A1 := DCO[I-1,K-1,0];
> B1 := DCO[I-1,K-1,1];
> C1 := DCO[I-1,K-1,2];
> D1 := 0;
> A2 := CO[I-1,K-1,0];
> B2 := CO[I-1,K-1,1];
> C2 := CO[I-1,K-1,2];
> D2 := CO[I-1,K-1,3];
> K := INTE(J,JJ);
> A3 := DCO[J-1,K-1,0];
> B3 := DCO[J-1,K-1,1];
> C3 := DCO[J-1,K-1,2];
> D3 := 0;
> A4 := CO[J-1,K-1,0];
> B4 := CO[J-1,K-1,1];
> C4 := CO[J-1,K-1,2];
> D4 := CO[J-1,K-1,3];
> A[I-1,J-1] := A[I-1,J-1]+XINT(XU,XL,AP[JJ-1],BP[JJ-1],CP[JJ-1],DP[JJ-
> 1],A1,B1,C1,D1,A3,B3,C3,D3)+XINT(XU,XL,AQ[JJ-1],BQ[JJ-1],CQ[JJ-1],DQ[JJ-
> 1],A2,B2,C2,D2,A4,B4,C4,D4);
> od;
> A[J-1,I-1] := A[I-1,J-1];
> od;
> fi;
> od;
> # Step 10
> for I from 1 to N1 do
> for J from I+1 to N2 do
> CC := A[J-1,I-1]/A[I-1,I-1];
> for K from I+1 to N3 do
> A[J-1,K-1] := A[J-1,K-1]-CC*A[I-1,K-1];
> od;
> A[J-1,I-1] := 0;
> od;
> od;
> C[N2-1] := A[N2-1,N3-1]/A[N2-1,N2-1];
> for I from 1 to N1 do
> J := N1-I+1;
> C[J-1] := A[J-1,N3-1];
> for KK from J+1 to N2 do
> C[J-1] := C[J-1]-A[J-1,KK-1]*C[KK-1];
> od;
> C[J-1] := C[J-1]/A[J-1,J-1];
> od;
> # Step 14
> # Output coefficients
> fprintf(OUP, `\nCoefficients:  c(1), c(2), ... , c(n+1)\n\n`);
> for I from 1 to N1 do
> fprintf(OUP, `  %12.6e \n`, C[I-1]);
> od;
> fprintf(OUP, `\n`);
> # compute and output value of the approximation at the nodes
> fprintf(OUP, `The approximation evaluated at the nodes:\n\n`);
> fprintf(OUP, `  Node            Value\n\n`);
> for I from 1 to N2 do
> S := 0;
> for J from 1 to N2 do
> J0 := max(J-2,1);
> J1 := min(J+2,N+2);
> SS := 0;
> if I < J0 or I >= J1 then
> S := S + C[J-1]*SS;
> else
> K := INTE(J,I);
> SS := 
> ((CO[J-1,K-1,3]*X[I-1]+CO[J-1,K-1,2])*X[I-1]+CO[J-1,K-1,1])*X[I-1]+CO[J-1,
> K-1,0];
> S := S + C[J-1]*SS;
> fi;
> od;
> fprintf(OUP, `%12.8f %12.8f\n`, X[I-1], S);
> od;
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> RETURN(0);
> end;
> alg116();