> restart;
> # HERMITE INTERPOLATION ALGORITHM 3.3
> #
> # TO OBTAIN THE COEFFICIENTS OF THE HERMITE INTERPOLATING
> # POLYNOMIAL H ON THE (N+1) DISTINCT NUMBERS X(0), ..., X(N)
> # FOR THE FUNCTION F:
> #
> # INPUT:   NUMBERS X(0), X(1), ..., X(N); VALUES F(X(0)), 
> #          F(X(1)),..., F(X(N)) AND F'(X(0)), F'(X(1)), ..., 
> #          F'(X(N)).
> #
> # OUTPUT:  NUMBERS Q(0,0), Q(1,1), ..., Q(2N + 1,2N + 1) WHERE
> #
> #          H(X) = Q(0,0) + Q(1,1) * ( X - X(0) ) + Q(2,2) *
> #                 ( X - X(0) )**2 + Q(3,3) * ( X - X(0) )**2 *
> #                 ( X - X(1) ) + Q(4,4) * ( X - X(0) )**2 *
> #                 ( X - X(1) )**2 + ... + Q(2N + 1,2N + 1) *
> #                 ( X - X(0) )**2 * ( X - X(1) )**2 * ... *
> #                 ( X - X(N - 1) )**2 * (X - X(N) ).
> alg033 := proc() local OK, FLAG, N, I, X, Q, A, NAME, INP, F, FP, Z, K, J, OUP, XX, S;
> printf(`This is Hermite 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\n`);
> printf(`Choose 1, 2, or 3 please\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 1 or FLAG = 2 or FLAG = 3 then
> OK := TRUE;
> fi;
> od;
> if FLAG = 1 then
> OK := FALSE;         
> while OK = FALSE do
> printf(`Input the number of data points minus 1\n`);            
> N := scanf(`%d`)[1];
> if N > 0 then
> OK := TRUE;               
> for I from 0 to N do                  
> printf(`Input X(%d), F(X(%d)), and `, I, I);
> printf(`F'(X(%d)) separated by space\n `, I);                  
> X[I] := scanf(`%f`)[1];
> Q[2*I,0] := scanf(`%f`)[1];
> Q[2*I+1,1] := 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 three columns?\n`);
> printf(`Enter Y or N\n`);
> A := scanf(`\n%c`)[1];
> if A = "Y" or A = "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 the number of data points minus 1\n`);
> N := scanf(`%d`)[1];
> if N > 0 then
> for I from 0 to N do
> X[I] := fscanf(INP, `%f`)[1];
> Q[2*I,0] := fscanf(INP, `%f`)[1];
> Q[2*I+1,1] := 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 three column `);
> printf(`form with the X values, F(X), and\n`);
> printf(`F'(X) values in the corresponding columns.\n`);
> printf(`The program will end so the input file can `);
> printf(`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: sin(x)\n`);
> F := scanf(`%a`)[1];
> FP := diff(F,x);
> F := unapply(F,x);
> FP := unapply(FP,x);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the number of data points minus 1\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];
> Q[2*I,0] := F(X[I]);
> Q[2*I+1,1] := FP(X[I]);
> od;
> OK := TRUE;
> else
> printf(`Number must be a positive integer\n`);
> fi;
> od;
> fi;
> if OK = TRUE then
> # Step 1
> for I from 0 to N do
> # Step 2
> Z[2*I] := X[I];
> Z[2*I+1] := X[I];
> Q[2*I+1,0] := Q[2*I,0];
> # Step 3
> if I <> 0 then
> Q[2*I,1] := (Q[2*I,0]-Q[2*I-1,0])/(Z[2*I]-Z[2*I-1]);
> fi;
> od;
> # Step 4
> K := 2*N+1;
> for I from 2 to K do
> for J from 2 to I do
> Q[I,J] := (Q[I,J-1]-Q[I-1,J-1])/(Z[I]-Z[I-J]);
> od;
> od;
> # Step 5
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please 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, `HERMITE INTERPOLATING POLYNOMIAL\n\n`);
> fprintf(OUP, `The input data follows:\n`);
> fprintf(OUP, `  X, F(X), F'(X)\n`);
> for I from 0 to N do
> fprintf(OUP, `  %12.10e %12.10e %12.10e\n`, X[I], Q[2*I,0], Q[2*I+1,1]);
> od;
> fprintf(OUP, `\nThe Coefficients of the Hermite Interpolation `);
> fprintf(OUP, `Polynomial\n`);
> fprintf(OUP, `in order of increasing exponent follow:\n\n`);
> for I from 0 to K do
> fprintf(OUP, `  %12.10e\n`, Q[I,I]);
> od;
> printf(`Do you wish to evaluate this polynomial?\n`);
> printf(`Enter Y or N\n`);
> A := scanf(`\n%c`)[1];
> if A = "Y" or A = "y" then
> printf(`Enter a point at which to evaluate\n`);
> XX := scanf(`%f`)[1];
> S := Q[K,K]*(XX-Z[K-1]);
> for I from 2 to K do
> J := K-I+1;
> S := (S+Q[J,J])*(XX-Z[J-1]);
> od;
> S := S + Q[0,0];
> fprintf(OUP, `x-value and interpolated-value\n`);
> fprintf(OUP, `  %12.10e %12.10e\n`, XX, S);
> fi;
> if OUP <> default then
> fclose(OUP);
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg033();