restart:# NATURAL CUBIC SPLINE ALGORITHM 3.4 ## 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) print(`This is the natural cubic spline interpolation.\134n`):OK := FALSE:while OK = FALSE doprint(`Choice of input method: `):print(`1. Input entry by entry from keyboard `):print(`2. Input data from a text file `):print(`3. Generate data using a function F with nodes entered `):print(`from keyboard `):print(`4. Generate data using a function F with nodes from `):print(`a text file `):print(`Choose 1, 2, 3, or 4 please `):FLAG := scanf(`%d`)[1]:print(`Entry = `):print(FLAG):if FLAG >= 1 and FLAG <= 4 thenOK := TRUE:fi:od:if FLAG = 1 thenOK := FALSE:while OK = FALSE doprint(`Input n `):N := scanf(`%d`)[1]:print(`n = `):print(N):if N > 0 thenOK := TRUE:for I1 from 0 to N doprint(`Input X(i) and F(X(i)) for i = `):print(I1):print(`separated by space `):X[I1] := scanf(`%f`)[1]:A[I1] := scanf(`%f`)[1]: print(`Entries are `):print(X[I1],A[I1]):od:else print(`Number must be a positive integer `):fi:od:fi:if FLAG = 2 thenprint(`Has a text file been created with the data in two `):print(`columns ? `):print(`Enter Y or N `):AA := scanf(` %c`)[1]:print(`Entry = `):print(AA):if AA = "Y" or AA = "y" thenprint(`Input the file name in the form - `):print(`drive:\134\134name.ext `):print(`For example: A:\134\134DATA.DTA `):NAME := scanf(`%s`)[1]:print(`Entry = `):print(NAME):INP := fopen(NAME,READ,TEXT):OK := FALSE:while OK = FALSE doprint(`Input n `):N := scanf(`%d`)[1]:print(`n = `):print(N):if N > 0 thenfor I1 from 0 to N doX[I1] := fscanf(INP, `%f`)[1]:A[I1] := fscanf(INP, `%f`)[1]:od:fclose(INP):OK := TRUE:else print(`Number must be a positive integer `):fi:od:elseprint(`Please create the input file in two column `):print(`form with the `):print(`X values and F(X) values in the corresponding columns. `):print(`The program will end so the input file can be created. `):OK := FALSE:fi:fi:if FLAG = 3 thenprint(`Input the function F(x) in terms of x. `):print(`For example: cos(x) `):F := scanf(`%a`)[1]:print(`F(x) = `):print(F):F := unapply(F,x):OK := FALSE:while OK = FALSE doprint(`Input n `):N := scanf(`%d`)[1]:print(`n = `):print(N):if N > 0 thenfor I1 from 0 to N doprint(`Input X(i) for i =`):print(I1):X[I1] := scanf(`%f`)[1]:print(`Entry = `):print(X[I1]):A[I1] := F(X[I1]):od:OK := TRUE:elseprint(`Number must be a positive integer `):fi:od:fi:if FLAG = 4 thenprint(`Has the text file with X-values been created? `):print(`Enter Y or N `):AA := scanf(` %c`)[1]:print(`Entry = `):print(AA):if AA = "Y" or AA = "y" thenprint(`Input the file name in the form - `):print(`drive:\134\134name.ext `):print(`For example: A:\134\134DATA.DTA `):NAME := scanf(`%s`)[1]:print(`Entry = `):print(NAME):INP := fopen(NAME,READ,TEXT):print(`Input the function F(x) in terms of x. `):print(`For example: cos(x) `):F := scanf(`%a`)[1]:print(`F(x) = `):print(F):F := unapply(F,x):OK := FALSE:while OK = FALSE doprint(`Input n `):N := scanf(`%d`)[1]:print(`n = `):print(N):if N > 0 thenOK := TRUE:for I1 from 0 to N doX[I1] := fscanf(INP, `%f`)[1]:A[I1] := F(X[I1]):od:fclose(INP):else print(`Number must be a positive integer `):fi:od:else print(`The program will end so the input file can be created`):OK := FALSE:fi:fi:if OK = TRUE thenM := N - 1:# Step 1for I1 from 0 to M doH[I1] := X[I1+1] - X[I1]:od:# Step 2# Use XA in place of ALPHAfor I1 from 1 to M doXA[I1] := 3.0*(A[I1+1]*H[I1-1]-A[I1]*(X[I1+1]-X[I1-1])+A[I1-1]*H[I1])/(H[I1]*H[I1-1]):od:# 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:XU[0] := 0:XZ[0] := 0:# Step 4for I1 from 1 to M doXL[I1] := 2*(X[I1+1]-X[I1-1])-H[I1-1]*XU[I1-1]:XU[I1] := H[I1]/XL[I1]:XZ[I1] := (XA[I1]-H[I1-1]*XZ[I1-1])/XL[I1]:od:# Step 5XL[N] := 1:XZ[N] := 0:C[N] := XZ[N]:# Step 6for I1 from 0 to M doJ := M-I1: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:D1[J] := (C[J+1] - C[J]) / (3.0 * H[J]):od:# Step 7print(`Select output destination `):print(`1. Screen `):print(`2. Text file `):print(`Enter 1 or 2 `):FLAG := scanf(`%d`)[1]:print(`Entry = `):print(FLAG):if FLAG = 2 thenprint(`Input the file name in the form - drive:\134\134name.ext `):print(`For example: A:\134\134OUTPUT.DTA `):NAME := scanf(`%s`)[1]:print(`Entry = `):print(NAME):OUP := fopen(NAME,WRITE,TEXT):else OUP := default:fi:fprintf(OUP, `NATURAL CUBIC SPLINE INTERPOLATION\134n\134n`):fprintf(OUP, `The numbers X(0), ..., X(N) are:\134n`):for I1 from 0 to N dofprintf(OUP, ` %11.8f`, X[I1]):od:fprintf(OUP, `\134n\134nThe coefficients of the spline on the subintervals `): fprintf(OUP, `are:\134n`):fprintf(OUP, `for I = 0, ..., N-1\134n`):fprintf(OUP, ` A(I) B(I) C(I) D(I)\134n`):for I1 from 0 to M dofprintf(OUP, ` %13.8f %13.8f %13.8f %13.8f \134n`, A[I1], B[I1], C[I1], D1[I1]):od:if OUP <> default thenfclose(OUP):print(`Output file `,NAME,` created successfully`):fi:fi:JSFH