alg081.txt

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

> restart;
> # PADE' RATIONAL APPROXIMATION ALGORITHM 8.1
> #
> # To obtain the rational approximation
> #
> #     r(x) = p(x) / q(x)
> #          = (p0 + p1*x + ... + pn*x^n) / (q0 + q1*x + ... + qm*x^m)
> #
> # for a given function f(x):
> #
> # INPUT  nonnegative integers m and n.
> #
> # OUTPUT  coefficients q0, q1, ... , qm, p0, p1, ... , pn.
> #
> # The coefficients of the Maclaurin polynomial a0, a1, ... , aN could
> # be calculated instead of input as is assumed in this program.
> alg081 := proc() local OK, LM, LN, BN, FLAG, I, AA, AAA, NAME, INP, N, M, NROW, NN, Q, P, J, A, IMAX, AMAX, JJ, IP, JP, NCOPY, I1, J1, XM, K, N1, PP, N2, SUM, KK, LL, OUP;
> printf(`This is Pade' Approximation.\n\n`);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input m and n\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 MacLaurin coefficients a[0], a[1], ... , a[N]\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]\n`);
> for I from 0 to BN do
> printf(`Input A[ %d ] \n`,I);
> AA[I] := scanf(`%f`)[1];
> od;
> fi;
> if FLAG = 2 then
> printf(`As many entries as desired can be placed\n`);
> printf(`on each line of the file 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 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 1 to N do
> NROW[I-1] := I;
> od;
> # Initialize row pointer for linear system
> NN := N-1;
> # Step 3
> Q[0] := 1;
> P[0] := AA[0];
> # Step 4
> # Set up a linear system but use A[i,j] in place of B[i,j]
> for I from 1 to N do
> # Step 5
> for J from 1 to I-1 do
> if J <= LN then
> A[I-1,J-1] := 0;
> fi;
> od;
> # Step 6
> if I <= LN then
> A[I-1,I-1] := 1;
> fi;
> # Step 7
> for J from I+1 to N do
> A[I-1,J-1] := 0;
> od;
> # Step 8
> for J from 1 to I do
> if  J <= LM then
> A[I-1,LN+J-1] := -AA[I-J];
> fi;
> od;
> # Step 9
> for J from LN+I+1 to N do
> A[I-1,J-1] := 0;
> od;
> # Step 10
> A[I-1,N] := AA[I];
> od;
> # Solve the linear system using partial pivoting.
> I := LN+1;
> # Step 11
> while OK = TRUE and I <= NN do
> # Step 12
> IMAX := NROW[I-1];
> AMAX := abs(A[IMAX-1,I-1]);
> IMAX := I;
> JJ := I+1;
> for IP from JJ to N do
> JP := NROW[IP-1];
> if abs(A[JP-1,I-1]) > AMAX then
> AMAX := abs(A[JP-1,I-1]);
> IMAX := IP;
> fi;
> od;
> # Step 13
> if AMAX <= 1.0e-20 then
> OK := false;
> else
> # Step 14
> # Simulate row interchange
> if NROW[I-1] <> NROW[IMAX-1] then
> NCOPY := NROW[I-1];
> NROW[I-1] := NROW[IMAX-1];
> NROW[IMAX-1] := NCOPY;
> fi;
> I1 := NROW[I-1];
> # Step 15
> # Perform elimination
> for J from JJ to M do
> J1 := NROW[J-1];
> # Step 16
> XM := A[J1-1,I-1]/A[I1-1,I-1];
> # Step 17
> for K from JJ to M do
> A[J1-1,K-1] := A[J1-1,K-1]-XM * A[I1-1,K-1];
> od;
> # Step 18
> A[J1-1,I-1] := 0;
> od;
> fi;
> I := I+1;
> od;
> if OK = TRUE then
> # Step 19
> N1 := NROW[N-1];
> if abs(A[N1-1,N-1]) <= 1.0e-20 then
> OK := FALSE;
> # System has no unique solution
> else
> # Step 20
> # Start backward substitution
> if LM > 0 then
> Q[LM] := A[N1-1,M-1]/A[N1-1,N-1];
> A[N1-1,M-1] := Q[LM];
> fi;
> PP := 1;
> # Step 21
> for K from LN+1 to NN do
> I := NN-K+LN+1;
> JJ := I+1;
> N2 := NROW[I-1];
> SUM := A[N2-1,N];
> for KK from JJ to N do
> LL := NROW[KK-1];
> SUM := SUM-A[N2-1,KK-1]*A[LL-1,M-1];
> od;
> # Step 22
> A[N2-1,M-1] := SUM/A[N2-1,I-1];
> Q[LM-PP] := A[N2-1,M-1];
> PP := PP+1;
> od;
> for K from 1 to LN do
> I := LN-K+1;
> N2 := NROW[I-1];
> SUM := A[N2-1,N];
> for KK from LN+1 to N do
> LL := NROW[KK-1];
> SUM := SUM-A[N2-1,KK-1]*A[LL-1,M-1];
> od;
> # Step 23
> # Procedure is completed.
> A[N2-1,M-1] := SUM;
> P[LN-K+1] := A[N2-1,M-1];
> od;
> 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, `PADE' 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;
> alg081();