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