program ALG034; { NATURAL CUBIC SPLINE To construct the cubic spline interpolant S for the function f, defined at the numbers x(0) < x(1) < ... < x(n), satisfying S''(x(0)) = S''(x(n)) = 0: INPUT: n; x(0), x(1), ..., x(N); either generate A(I) = f(x(I)) for I = 0, 1, ..., n or input A(I) for I = 0, 1, ..., n. OUTPUT: A(J), B(J), C(J), D(J) for J = 0, 1, ..., n - 1. NOTE: S(x) = A(J) + B(J) * ( x - x(J) ) + C(J) * ( x - x(J) )**2 + D(J) * ( x - x(J) )**3 for x(J) <= x < x(J + 1) } var X,A,B,C,D,H,XA,XL,XU,XZ : array [ 0..25 ] of real; XX, S, Y : real; FLAG,N,I,J,M : integer; OK : boolean; AA : char; INP,OUP : text; NAME : string [ 14 ]; { Change function F to solve a new problem } function F ( X : real ) : real; begin F := exp(2.0*X) end; procedure INPUT; begin writeln('This is the natural cubic spline interpolation.'); OK := false; while ( not OK ) do begin writeln ('Choice of input method: '); writeln ('1. Input entry by entry from keyboard '); writeln ('2. Input data from a text file '); write ('3. Generate data using a function F with nodes entered '); writeln ('from keyboard '); write ('4. Generate data using a function F with nodes from '); writeln ('a text file '); writeln ('Choose 1, 2, 3, or 4 please '); readln ( FLAG ); if ( FLAG >= 1 ) and ( FLAG <= 4 ) then OK := true end; case FLAG of 1 : begin OK := false; while ( not OK ) do begin writeln ('Input n '); readln ( N ); if ( N > 0 ) then begin OK := true; for I := 0 to N do begin write ('Input X(',I,') and F(X(',I,')) '); writeln ('separated by space '); readln ( X[I] , A[I] ) end end else writeln ('Number must be a positive integer ') end end; 2 : begin write ('Has a text file been created with data in two '); writeln ('columns? '); writeln ('Enter Y or N '); readln ( AA ); if ( AA = 'Y' ) or ( AA = 'y' ) then begin write ('Input the file name in the form - '); writeln ('drive:name.ext '); writeln ('for example: A:DATA.DTA '); readln ( NAME ); assign ( INP , NAME ); reset ( INP ); OK := false; while ( not OK ) do begin writeln ('Input n '); readln ( N ); if ( N > 0 ) then begin OK := true; for I := 0 to N do readln ( INP , X[I] , A[I] ); close ( INP ) end else writeln ('Number must be a positive integer ') end end else begin write ('Please create the input file in two column '); writeln ('form with the '); write ('X values and F(X) values in the '); writeln ('corresponding columns '); write ('The program will end so the input file can '); writeln ('be created. '); OK := false end end; 3 : begin write ('Has the function F been created in the program '); writeln ('immediately proceding the INPUT procedure? '); writeln ('Enter Y or N '); readln ( AA ); if ( AA = 'Y' ) or ( AA = 'y' ) then begin OK := false; while ( not OK ) do begin writeln ('Input n '); readln ( N ); if ( N > 0 ) then begin OK := true; for I := 0 to N do begin writeln ('Input X(',I,') '); readln ( X[I] ); A[I] := F( X[I] ) end end else writeln ('Number must be a positive integer ') end end else begin write ('The program will end so that the function F '); writeln ('can be created '); OK := false end end; 4 : begin write ('Has the text file with X-values been created and '); writeln ('has the function F been created in the program '); writeln ('immediately preceding the INPUT procedure? '); writeln ('Enter Y or N '); readln ( AA ); if ( AA = 'Y' ) or ( AA = 'y') then begin write ('Input the file name in the form - '); writeln ('drive:name.ext '); writeln ('for example: A:DATA.DTA '); readln ( NAME ); assign ( INP , NAME ); reset ( INP ); OK := false; while ( not OK ) do begin writeln ('Input n '); readln ( N ); if ( N > 0 ) then begin OK := true; for I := 0 to N do begin readln ( INP , X[I] ); A[I] := F( X[I] ) end; close ( INP ) end else writeln ('Number must be a positive integer ') end end else begin write ('Please create the input file with one entry '); writeln ('per row and create the function F. '); write ('The program will end so the input file and '); writeln ('F can be created. '); OK := false end end end end; procedure OUTPUT; begin writeln ('Select output destination: '); writeln ('1. Screen '); writeln ('2. Text file '); writeln ('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,'NATURAL CUBIC SPLINE INTERPOLATION'); writeln(OUP); writeln ( OUP, 'The numbers X(0), ..., X(N) are: '); for I := 0 to N do write ( X[I]:12:8 ); writeln ( OUP ); writeln ( OUP ); writeln(OUP,'The coefficients of the spline on the subintervals are:'); writeln ( OUP, 'for I = 0, ..., N-1 '); write(OUP, ' A(I) B(I) C(I) '); writeln(OUP,' D(I) '); for I := 0 to M do writeln ( OUP, A[I]:14:8, B[I]:14:8, C[I]:14:8, D[I]:14:8 ); close ( OUP ) end; begin INPUT; if ( OK ) then begin M := N - 1; { STEP 1 } for I := 0 to M do H[I] := X[I+1] - X[I]; { STEP 2 } for I := 1 to M do XA[I] := 3.0 * ( A[I+1] * H[I-1] - A[I] * ( X[I+1] - X[I-1] ) + A[I-1] * H[I] ) / ( H[I] * H[I-1 ] ); { STEP 3 } { STEPS 3, 4, 5 and part of 6 solve the tridiagonal system using Algorithm 6.7 use XL, XU, XZ in place of L, MU, Z resp. } XL[0] := 1.0; XU[0] := 0.0; XZ[0] := 0.0; { STEP 4 } for I := 1 to M do begin XL[I] := 2.0 * ( X[I+1] - X[I-1] ) - H[I-1] * XU[I-1]; XU[I] := H[I] / XL[I]; XZ[I] := ( XA[I] - H[I-1] * XZ[I-1] ) / XL[I] end; { STEP 5 } XL[N] := 1.0; XZ[N] := 0.0; C[N] := XZ[N]; { STEP 6 } for I := 0 to M do begin J := M - I; C[J] := XZ[J] - XU[J] * C[J+1]; B[J] := ( A[J+1] - A[J] ) / H[J] - H[J] * ( C[J+1] + 2.0 * C[J] ) / 3.0; D[J] := ( C[J+1] - C[J] ) / (3.0 * H[J] ) end; { STEP 7 } OUTPUT end end.