program alg116; { CUBIC SPLINE RAYLEIGH-RITZ METHOD 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 To change problems do the following: 1. Change functions P, Q and F in the subprograms named P,Q,F 2. Change dimensions in main program to values given by dimension A(N+2,N+3),X(N+2),C(N+2),CO(N+2,4,4), DCO(N+2,4,3),AF(N+1),BF(N+1),CF(N+1),DF(N+1), AP(N+1),BP(N+1),CP(N+1),DP(N+1), AQ(N+1),BQ(N+1),CQ(N+1),DQ(N+1) 3. Change dimensions in subroutine COEF to dimension AA(N+2),BB(N+2),CC(N+2),DD(N+2), H(N+2),XA(N+2),XL(N+2),XU(N+2),XZ(N+2) 4. Change call statements for subroutine COEF to reflect the derivatives of P,Q,F at X(1)=0 and X(N+2)=1 GENERAL OUTLINE 1. Nodes labelled X(I)=(I-1)*H, 1 <= I <= N+2, where H=1/(N+1) so that zero subscripts are avoided 2. The functions PHI(I) and PHI'(I) are shifted so that PHI(1) and PHI'(1) are centered at X(1), PHI(2) and PHI'(2) are centered at X(2), . . . , PHI(N+2) and PHI'(N+2) are centered at (X(N+2)---for example, PHI(3) = S((X-X(3))/H) = S(X/H + 2) 3. The functions PHI(I) are represented in terms of their coefficients in the following way: (PHI(I))(X) = CO(I,K,1) + CO(I,K,2)*X + CO(I,K,3)*X**2 CO(I,K,4) *X**3 for X(J) <= X <= X(J+1) where K=1 IF J=I-2, K=2 IF J=I-1, K=3 IF J=I, K=4 IF J=I+1 since PHI(I) is nonzero only between X(I-2) and X(I+2) unless I = 1, 2, N+1 or N+2 (see subroutine PHICO) 4. The derivative of PHI(I) denoted PHI'(I) is represented as in 3. By its coefficients DCO(I,K,L), L = 1, 2, 3 (See subroutine DPHICO). 5. The functions P,Q and F are represented by their cubic spline interpolants using clamped boundary conditions (see Algorithm 3.5). Thus, for X(I) <= X <= X(I+1) we use AF(I) + BF(I)*X + CF(I)*X**2 + DF(I)*X**3 to represent F(X). Similarly, AP,BP,CP,DP are used for P and AQ,BQ,CQ,DQ are used for Q. (see subroutine COEF). 6. The integrands in STEPS 6 and 9 are replaced by products of cubic polynomial approximations on each subinterval of length H and the integrals of the resulting polynomials are computed exactly. (see subroutine XINT). } type VEC = array [1..20] of real; var A : array [1..21,1..22] of real; CO : array [1..21,1..4,1..4] of real; DCO : array [1..21,1..4,1..3] of real; X, C : array [1..21] of real; AF, BF, CF, DF, AP, BP, CP, DP, AQ, BQ, CQ, DQ : VEC; T,H, XU, XL, A1, B1, C1, D1, A2, B2, C2, D2 : real; A3, B3, C3, D3, A4, B4, C4, D4, CC, S, SS : real; FPL,FPR,PPL,PPR,QPL,QPR : real; N, N1, N2, N3, I, J, K, II, JJ, J0, J1, J2 : integer; FLAG,KK,K2, K3, JJ1, JJ2 : integer; OK : boolean; NAME : string[14]; AA : char; OUP : text; function F (X : real) : real; begin F := 2*pi*pi*sin(pi*X) end; function P (x : real) : real; begin P := 1 end; function Q (x : real) : real; begin Q := pi*pi end; procedure INPUT; begin writeln('This is the Cubic Spline Rayleigh-Ritz Method.'); OK := false; writeln ('Have functions P, Q and F been created '); writeln ('immediately preceding the INPUT procedure? '); writeln ('Answer Y or N. '); readln ( AA ); if ( AA = 'Y' ) or ( AA = 'y' ) then begin while ( not OK ) do begin write ('Input a positive integer n, where x(0) = 0, '); writeln ('..., x(n+1) = 1. '); readln ( N ); if ( N <= 0 ) then writeln ('Number must be a positive integer. ') else OK := true end; writeln('Input the derivative of f at 0 and at 1'); readln(FPL,FPR); writeln('Input the derivative of p at 0 and at 1'); readln(PPL,PPR); writeln('Input the derivative of q at 0 and at 1'); readln(QPL,QPR); end else writeln ('The program will end so that P, Q, F can be created. ') end; procedure OUTPUT; begin writeln ('Choice of output method: '); writeln ('1. Output to screen '); writeln ('2. Output to text file '); writeln ('Please enter 1 or 2. '); readln ( FLAG ); if ( FLAG = 2 ) then begin writeln ('Input the file name in the form - drive:name.ext, '); writeln('for example: A:OUTPUT.DTA'); readln ( NAME ); assign ( OUP, NAME ) end else assign ( OUP, 'CON' ); rewrite ( OUP ); writeln(OUP,'CUBIC SPLINE RAYLEIGH-RITZ METHOD'); writeln(OUP) end; function MINO ( I,J : integer) : integer; begin MINO := I; if (J < I) then MINO := J end; function MAXO (I,J : integer) : integer; begin MAXO := I; if (J > I) then MAXO := J end; function INTE (J,JJ : integer) : integer; begin INTE := JJ - J + 3 end; procedure DPHICO ( I,J : integer; var A,B,C : real); var EE,E, AA, BB, CC : real; OK : boolean; begin A := 0.0; B := 0.0; C := 0.0; E := I - 1.0; OK := true; if ( (J < I-2) or (J >= I+2) ) then OK := false; if ( (I = 1) and (J < I) ) then OK := false; if ( (I = 2) and (J < I-1) ) then OK := false; if ( (I = N+1) and (J > N+1) ) then OK := false; if ( (I = N+2) and (J >= N+2) ) then OK := false; if ( OK ) then begin if (J <= I-2) then begin A := ((E-4.0)*E+4.0)/(8.0*H); B := (-E+2.0)/(4.0*H*H); C := 1.0/(8.0*H*H*H); OK := false end else begin if ( J > I) then begin A := ((-E-4.0)*E-4.0)/(8.0*H); B := (E+2.0)/(4.0*H*H); C := -1.0/(8.0*H*H*H); OK := false end else begin if (J > I-1) then begin A := (3.0*E+4.0)*E/(8.0*H); B := (-3.0*E-2.0)/(4.0*H*H); C := 3.0/(8.0*H*H*H); if ( (I <> 1) and (I <> N+1) ) then OK := false end else begin A := (-3.0*E+4.0)*E/(8.0*H); B := (3.0*E-2.0)/(4.0*H*H); C := -3.0/(8.0*H*H*H); if ( (I <> 2) and (I <> N+2) ) then OK := false end end end end; if ( OK ) then begin if ( I <= 2 ) then begin AA := -1.0/(8.0*H); BB := 1.0/(4.0*H*H); CC := -1.0/(8.0*H*H*H); if ( I = 2 ) then begin A := A - AA; B := B - BB; C := C - CC end else begin A := A - 4.0*AA; B := B - 4.0*BB; C := C - 4.0*CC end end else begin EE := N+2.0; AA := ((EE-4.0)*EE+4.0)/(8.0*H); BB := (-EE+2.0)/(4.0*H*H); CC := 1.0/(8.0*H*H*H); if ( I = N+1 ) then begin A := A - AA; B := B - BB; C := C - CC end else begin A := A - 4.0*AA; B := B - 4.0*BB; C := C - 4.0*CC end end end end; procedure PHICO ( I,J : integer; var A,B,C,D : real); var EE,E, AA, BB, CC,DD : real; OK : boolean; begin A := 0.0; B := 0.0; C := 0.0; D := 0.0; E := I - 1.0; OK := true; if ( (J < I-2) or (J >= I+2) ) then OK := false; if ( (I = 1) and (J < I) ) then OK := false; if ( (I = 2) and (J < I-1) ) then OK := false; if ( (I = N+1) and (J > N+1) ) then OK := false; if ( (I = N+2) and (J >= N+2) ) then OK := false; if ( OK ) then begin if (J <= I-2) then begin A := (((-E+6.0)*E-12.0)*E+8.0)/24.0; B := ((E-4.0)*E+4.0)/(8.0*H); C := (-E+2.0)/(8.0*H*H); D := 1.0/(24.0*H*H*H); OK := false end else begin if (J > I) then begin A := (((E+6.0)*E+12.0)*E+8.0)/24.0; B := ((-E-4.0)*E-4.0)/(8.0*H); C := (E+2.0)/(8.0*H*H); D := -1.0/(24.0*H*H*H); OK := false end else begin if ( J > I-1) then begin A := ((-3.0*E-6.0)*E*E+4.0)/24.0; B := (3.0*E+4.0)*E/(8.0*H); C := (-3.0*E-2.0)/(8.0*H*H); D := 1.0/(8.0*H*H*H); if ( (I <> 1) and (I <> N+1) ) then OK := false end else begin A := ((3.0*E-6.0)*E*E+4.0)/24.0; B := (-3.0*E+4.0)*E/(8.0*H); C := (3.0*E-2.0)/(8.0*H*H); D := -1.0/(8.0*H*H*H); if (( I <> 2) and (I <> N+2)) then OK := false end end end end; if ( OK ) then begin if ( I <= 2 ) then begin AA := 1.0/24.0; BB := -1.0/(8.0*H); CC := 1.0/(8.0*H*H); DD := -1.0/(24.0*H*H*H); if ( I = 2 ) then begin A := A - AA; B := B - BB; C := C - CC; D := D - DD end else begin A := A - 4.0*AA; B := B - 4.0*BB; C := C - 4.0*CC; D := D - 4.0*DD end end else begin EE := N+2.0; AA := (((-EE+6.0)*EE-12.0)*EE+8.0)/24.0; BB := ((EE-4.0)*EE+4.0)/(8.0*H); CC := (-EE+2.0)/(8.0*H*H); DD := 1.0/(24.0*H*H*H); if ( I = N+1 ) then begin A := A - AA; B := B - BB; C := C - CC; D := D - DD end else begin A := A - 4.0*AA; B := B - 4.0*BB; C := C - 4.0*CC; D := D - 4.0*DD end end end end; function XINT(XU,XL,A1,B1,C1,D1,A2,B2,C2,D2,A3,B3,C3,D3 : real) : real; var C : array [1..20] of real; AA,BB,CC,DD,EE,FF,GG,XHIGH,XLOW : real; I : integer; begin 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[10] := AA*A3; C[9] := (AA*B3+BB*A3)/2.0; C[8] := (AA*C3+BB*B3+CC*A3)/3.0; C[7] := (AA*D3+BB*C3+CC*B3+DD*A3)/4.0; C[6] := (BB*D3+CC*C3+DD*B3+EE*A3)/5.0; C[5] := (CC*D3+DD*C3+EE*B3+FF*A3)/6.0; C[4] := (DD*D3+EE*C3+FF*B3+GG*A3)/7.0; C[3] := (EE*D3+FF*C3+GG*B3)/8.0; C[2] := (FF*D3+GG*C3)/9.0; C[1] := (GG*D3)/10.0; XHIGH := 0.0; XLOW := 0.0; for I := 1 to 10 do begin XHIGH := (XHIGH + C[I])*XU; XLOW := (XLOW + C[I])*XL end; XINT := XHIGH - XLOW end; procedure COEF(L,N,M : integer; FPO, FPN : real; var A,B,C,D : VEC); var AA,BB,CC,DD,XA,XL,XU,XZ : array[1..25] of real; I, J : integer; begin for I := 1 to N do case L of 1 : AA[I] := F(X[I]); 2 : AA[I] := P(X[I]); 3 : AA[I] := Q(X[I]) end; XA[1] := 3.0*(AA[2]-AA[1])/H - 3.0*FPO; XA[N] := 3.0*FPN - 3.0*(AA[N]-AA[N-1])/H; XL[1] := 2.0*H; XU[1] := 0.5; XZ[1] := XA[1]/XL[1]; for I := 2 to M do begin XA[I] := 3.0*(AA[I+1]-2.0*AA[I]+AA[I-1])/H; XL[I] := H*(4.0-XU[I-1]); XU[I] := H/XL[I]; XZ[I] := (XA[I]-H*XZ[I-1])/XL[I] end; XL[N] := H*(2.0-XU[N-1]); XZ[N] := (XA[N]-H*XZ[N-1])/XL[N]; CC[N] := XZ[N]; for I := 1 to M do begin J := N-I; CC[J] := XZ[J]-XU[J]*CC[J+1]; BB[J] := (AA[J+1]-AA[J])/H -H*(CC[J+1]+2.0*CC[J])/3.0; DD[J] := (CC[J+1]-CC[J])/(3.0*H) end; for I := 1 to M do begin A[I] := ((-DD[I]*X[I]+CC[I])*X[I]-BB[I])*X[I]+AA[I]; B[I] := (3.0*DD[I]*X[I]-2.0*CC[I])*X[I]+BB[I]; C[I] := CC[I]-3.0*DD[I]*X[I]; D[I] := DD[I] end end; begin INPUT; if (OK) then begin OUTPUT; { STEP 1 } H := 1.0/(N+1.0); N1 := N+1; N2 := N+2; N3 := N+3; { Initialize matrix A at zero, not that A[I, N+3] = B[I] } for I := 1 to N2 do for J := 1 to N3 do A[I,J] := 0.0; { STEP 2 } { X[1]=0,...,X[I] = (I-1)*H,...,X[N+1] = 1 - H, X[N+2] = 1 } for I := 1 to N2 do X[I] := (I-1.0)*H; { STEPS 3 and 4 are implemented in what follows. Initialize coefficients CO[I,J,K], DCO[I,J,K] } for I := 1 to N2 do for J := 1 to 4 do begin for K := 1 to 4 do begin CO[I,J,K] := 0.0; if (K <> 4) then DCO[I,J,K] := 0.0 end; { JJ corresponds the coefficients of phi and phi' to the proper interval involving J } JJ := I+J-3; PHICO(I,JJ,CO[I,J,1],CO[I,J,2],CO[I,J,3],CO[I,J,4]); DPHICO(I,JJ,DCO[I,J,1],DCO[I,J,2],DCO[I,J,3]) end; { Output the basis functions. } writeln(OUP,'Basis Function: A + B*X + C*X**2 + D*X**3'); writeln(OUP); writeln(OUP,' A',' B' ,' C',' D'); writeln(OUP); for I := 1 to N2 do begin writeln(OUP,'phi( ',I,' )'); writeln(OUP); for J := 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 begin JJ1 := I+J-3; JJ2 := I+J-2; write(OUP,'On (X( ',JJ1,' ), X( ',JJ2,' ) '); for K := 1 to 4 do write(OUP,' ',CO[I,J,K]:12:8,' '); writeln(OUP); writeln(OUP) end end; { Obtain coefficients for F, P, Q - NOTE _ the fourth and fifth arguements in COEF are the derivatives of F, P, or Q evaluated at 0 and 1 respectively. } COEF(1,N2,N1,FPL,FPR,AF,BF,CF,DF); COEF(2,N2,N1,PPL,PPR,AP,BP,CP,DP); COEF(3,N2,N1,QPL,QPR,AQ,BQ,CQ,DQ); { STEPS 5 - 9 are implemented in what follows } for I := 1 to N2 do { indices for limits of integration for A[I,I] and B[I] } begin J1 := MINO(I+2,N+2); J0 := MAXO(I-2,1); J2 := J1 - 1; { integrate over each subinterval where phi(I) nonzero } for JJ := J0 to J2 do begin { limits of integration for each call } XU := X[JJ+1]; XL := X[JJ]; { coefficients of bases } K := INTE(I,JJ); A1 := DCO[I,K,1]; B1 := DCO[I,K,2]; C1 := DCO[I,K,3]; D1 := 0.0; A2 := CO[I,K,1]; B2 := CO[I,K,2]; C2 := CO[I,K,3]; D2 := CO[I,K,4]; { call subprogram for integrations } A[I,I] := A[I,I]+ XINT(XU,XL,AP[JJ],BP[JJ],CP[JJ],DP[JJ], A1,B1,C1,D1,A1,B1,C1,D1)+ XINT(XU,XL,AQ[JJ],BQ[JJ],CQ[JJ],DQ[JJ], A2,B2,C2,D2,A2,B2,C2,D2); A[I,N+3] := A[I,N+3] + XINT(XU,XL,AF[JJ],BF[JJ],CF[JJ],DF[JJ], A2,B2,C2,D2,1.0,0.0,0.0,0.0) end; { compute A[I,J] for J = I+1, ..., Min(I+3,N+2) } K3 := I+1; if (K3 <= N2) then begin K2 := MINO(I+3,N+2); for J := K3 to K2 do begin J0 := MAXO(J-2,1); for JJ := J0 to J2 do begin XU := X[JJ+1]; XL := X[JJ]; K := INTE(I,JJ); A1 := DCO[I,K,1]; B1 := DCO[I,K,2]; C1 := DCO[I,K,3]; D1 := 0.0; A2 := CO[I,K,1]; B2 := CO[I,K,2]; C2 := CO[I,K,3]; D2 := CO[I,K,4]; K := INTE(J,JJ); A3 := DCO[J,K,1]; B3 := DCO[J,K,2]; C3 := DCO[J,K,3]; D3 := 0.0; A4 := CO[J,K,1]; B4 := CO[J,K,2]; C4 := CO[J,K,3]; D4 := CO[J,K,4]; A[I,J] := A[I,J] + XINT(XU,XL,AP[JJ],BP[JJ], CP[JJ],DP[JJ],A1,B1,C1,D1, A3,B3,C3,D3)+ XINT(XU,XL,AQ[JJ],BQ[JJ], CQ[JJ],DQ[JJ],A2,B2,C2,D2, A4,B4,C4,D4) end; A[J,I] := A[I,J] end end end; { STEP 10 } for I := 1 to N1 do begin for J := I+1 to N2 do begin CC := A[J,I]/A[I,I]; for K := I+1 to N3 do A[J,K] := A[J,K]-CC*A[I,K]; A[J,I] := 0.0 end end; C[N2] := A[N2,N3]/A[N2,N2]; for I := 1 to N1 do begin J := N1-I+1; C[J] := A[J,N3]; for KK := J+1 to N2 do C[J] := C[J]-A[J,KK]*C[KK]; C[J] := C[J]/A[J,J] end; { STEP 11 output coefficients } writeln(OUP); writeln(OUP,'Coefficients: c(1), c(2), ... , c(n)'); writeln(OUP); for I := 1 to N1 do writeln(OUP,' ',C[I]:12,' '); writeln(OUP); { compute and output value of the approximation at the nodes } writeln(OUP,'The approximation evaluated at the nodes: '); writeln(OUP); writeln(OUP,' Node Value'); writeln(OUP); for I := 1 to N2 do begin S := 0.0; for J := 1 to N2 do begin J0 := MAXO(J-2,1); J1 := MINO(J+2,N+2); SS := 0.0; if ((I < J0) or (I >= J1)) then S := S + C[J]*SS else begin K := INTE(J,I); SS := ((CO[J,K,4]*X[I]+CO[J,K,3])*X[I]+ CO[J,K,2])*X[I]+CO[J,K,1]; S := S + C[J]*SS end end; writeln(OUP,X[I]:12:8,' ',S:12:8) end; { STEP 12 } close(OUP) end end.