alg035.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> restart;
> # CLAMPED CUBIC SPLINE ALGORITHM 3.5
> #
> # 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)) = f'(x(0)) and S'(x(n)) = f'(x(n)):
> #
> # 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;
> #          FP0 = f'(x(0)); FPN = f'(x(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)
> alg035 := proc() local OK, FLAG, N, I, X, A, AA, NAME, INP, F, FP0, FPN, M, H, XA, XL, XU, XZ, C, J, B, D, OUP;
> printf(`This is Clamped 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 a 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 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
> printf(`Enter F'(X(0)) and F'(X(N)) separated by a space\n`);
> FP0 := scanf(`%f`)[1];
> FPN := scanf(`%f`)[1];
> 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 instead of ALPHA
> XA[0] := 3.0 * (A[1] - A[0]) / H[0] - 3.0 * FP0;
> XA[N] := 3.0 * FPN - 3.0 * (A[N] - A[N-1]) / H[N-1];
> # Step 3
> 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 4
> # STEPS 4, 5, 6 and part of 7 solve the tridiagonal system using 
> # Algorithm 6.7.
> # use XL, XU, XZ in place of L, MU, Z resp.
> XL[0] := 2.0 * H[0];
> XU[0] := 0.5;
> XZ[0] := XA[0] / XL[0];
> # Step 5
> for I from 1 to M do
> 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];
> od;
> # Step 6
> XL[N] := H[N-1] * (2.0 - XU[N-1]);
> XZ[N] := (XA[N] - H[N-1] * XZ[N-1]) / XL[N];
> C[N] := XZ[N];
> # Step 7
> for I from 1 to N do
> J := N - 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 8
> 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, `CLAMPED 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;
> alg035();