alg082.txt

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

> restart;
> # CHEBYSHEV RATIONAL APPROXIMATION ALGORITHM 8.2
> #
> # To obtain the rational approximation
> #
> # rT(x) = (p0*T0 + p1*T1 +...+ pn*Tn) / (q0*T0 + q1*T1 +...+ qm*Tm)
> #
> # for a given function f(x):
> #
> # INPUT  nonnegative integers m and n.
> #
> # OUTPUT  coefficients q0, q1, ... , qm, p0, p1, ... , pn.
> #
> # The coefficients of the Chebyshev expansion a0, a1, ..., aN could
> # be calculated instead of input as is assumed in this program.
> alg082 := proc() local OK, LM, LN, BN, FLAG, I, AA, AAA, NAME, INP, N, M, NROW, NN, Q, J, A, PP, IMAX, AMAX, JJ, IP, JP, NCOPY, I1, J1, XM, K, N1, N2, SUM, KK, LL, P, OUP;
> printf(`This is Chebyshev Approximation.\n\n`);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input m and n separated by a space\n`);
> LM := scanf(`%d`)[1];
> LN := scanf(`%d`)[1];
> BN := LM+LN;
> if LM >= 0 and LN >= 0 then
> OK := TRUE;
> else printf(`m and n must both be nonnegative.\n`);
> fi;
> if LM = 0 and LN = 0 then
> OK := FALSE;
> printf(`Not both m and n can be zero\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`The Chebyshev coefficients a[0], a[1], ... , a[N+m]\n`);
> printf(`are to be input.\n`);
> 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(`Choose 1 or 2 please\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 1 or FLAG = 2 then
> OK := TRUE;
> fi;
> od;
> if FLAG = 1 then
> printf(`Input in order a[0] to a[N+m]\n`);
> for I from 0 to BN+LM do
> printf(`Input A[%d]\n`, I);
> AA[I] := scanf(`%f`)[1];
> od;
> fi;
> if FLAG = 2 then
> printf(`The text file may contain as many entries\n`);
> printf(`per line as desired each separated by blank.\n`);
> printf(`Has such a text file been created?\n`);
> printf(`Enter Y or N\n`);
> AAA := scanf(`\n%c`)[1];
> if AAA = "Y" or AAA = "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);
> for I from 0 to BN+LM do
> AA[I] := fscanf(INP, `%f`)[1];
> od;fclose(INP);
> else
> printf(`Please create the input file.\n`);
> printf(`The program will end so the input file can `);
> printf(`be created.\n`);
> OK := FALSE;
> fi;
> fi;
> if OK = TRUE then
> # Step 1
> N := BN;
> M := N+1;
> # Step 2 - performed during input
> for I from 0 to N do
> NROW[I] := I;
> od;
> # Initialize row pointer
> NN := N-1;
> # Step 3
> Q[0] := 1.0;
> # Step 4
> # Set up a linear system with A used in place of B
> for I from 0 to N do
> # Step 5
> for J from 0 to I do
> if J <= LN then
> A[I,J] := 0;
> fi;
> od;
> # Step 6
> if I <= LN then
> A[I,I] := 1.0;
> fi;
> # Step 7
> for J from I+1 to LN do
> A[I,J] := 0;
> od;
> # Step 8
> for J from LN+1 to N do
> if I <> 0 then
> PP := I-J+LN;
> if PP < 0 then
> PP := -PP;
> fi;
> A[I,J] := -(AA[I+J-LN]+AA[PP])/2.0;
> else
> A[I,J] := -AA[J-LN]/2.0;
> fi;
> od;
> A[I,N+1] := AA[I];
> od;
> # Step 9
> A[0,N+1] := A[0,N+1]/2.0;
> # Steps 10 - 21 solve the linear system using partial pivoting
> I := LN+1;
> # Step 10
> while OK = TRUE and I <= NN do
> # Step 11
> IMAX := NROW[I];
> AMAX := abs(A[IMAX,I]);
> IMAX := I;
> JJ := I+1;
> for IP from JJ to N do
> JP := NROW[IP];
> if abs(A[JP,I]) > AMAX then
> AMAX := abs(A[JP,I]);
> IMAX := IP;
> fi;
> od;
> # Step 12
> if AMAX <= 1.0e-20 then
> OK := false;
> else
> # Step 13
> # Simulate row interchange
> if NROW[I] <> NROW[IMAX] then
> NCOPY := NROW[I];
> NROW[I] := NROW[IMAX];
> NROW[IMAX] := NCOPY;
> fi;
> I1 := NROW[I];
> # Step 14
> # Perform elimination
> for J from JJ to M do
> J1 := NROW[J];
> # Step 15
> XM := A[J1,I]/A[I1,I];
> # Step 16
> for K from JJ to M do
> A[J1,K] := A[J1,K]-XM*A[I1,K];
> od;
> # Step 17
> A[J1,I] := 0;
> od;
> fi;
> I := I+1;
> od;
> if OK = TRUE then
> # Step 18
> N1 := NROW[N];
> if abs(A[N1,N]) <= 1.0e-20 then
> OK := false;
> # System has no unique solution
> else
> # Step 19
> # Start backward substitution
> if LM > 0 then
> Q[LM] := A[N1,M]/A[N1,N];
> A[N1,M] := Q[LM];
> fi;
> PP := 1;
> # Step 20
> for K from LN+1 to NN do
> I := NN-K+LN+1;
> JJ := I+1;
> N2 := NROW[I];
> SUM := A[N2,N+1];
> for KK from JJ to N do
> LL := NROW[KK];
> SUM := SUM - A[N2,KK] * A[LL,M];
> od;
> A[N2,M] := SUM / A[N2,I];
> Q[LM-PP] := A[N2,M];
> PP := PP+1;
> od;
> # Step 21
> for K from 0 to LN do
> I := LN-K;
> N2 := NROW[I];
> SUM := A[N2,N+1];
> for KK from LN+1 to N do
> LL := NROW[KK];
> SUM := SUM-A[N2,KK]*A[LL,M];
> od;
> A[N2,M] := SUM;
> P[LN-K] := A[N2,M];
> od;
> # Step 22
> # Procedure completed
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to 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, `CHEBYSHEV RATIONAL APPROXIMATION\n\n`);
> fprintf(OUP, `Denominator Coefficients Q[0], ..., Q[M] \n`);
> for I from 0 to LM do
> fprintf(OUP, ` %11.8f`, Q[I]);
> od;
> fprintf(OUP, `\n`);
> fprintf(OUP, `Numerator Coefficients P[0], ..., P[N]\n`);
> for I from 0 to LN do
> fprintf(OUP, ` %11.8f`, P[I]);
> od;
> fprintf(OUP, `\n`);
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> fi;
> if OK = FALSE then
> printf(`System has no unique solution\n`);
> fi;
> fi;
> RETURN(0);
> end;
> alg082();