Download Pade Rational Approximation-Numerical Analysis-MATLAB Code and more Exercises Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
% PADE RATIONAL APPROXIMATION ALGORITHM 8.
% To obtain the rational approximation % % r(x) = p(x) / q(x) % = (p0 + p1x + ... + Pnx^n) / (q0 + q1x + ... + qmx^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, ... could % be calculated instead of input as is assumed in this program. syms('OK', 'LM', 'LN', 'BN', 'FLAG', 'I', 'AA', 'AAA'); syms('NAME', 'INP', 'N', 'M', 'NROW', 'NN', 'Q', 'P', 'J'); syms('A', 'IMAX', 'AMAX', 'JJ', 'IP', 'JP', 'NCOPY', 'I1'); syms('J1', 'XM', 'K', 'N1', 'PP', 'N2', 'SUM', 'KK', 'LL', 'OUP'); TRUE = 1; FALSE = 0; fprintf(1,'This is Pade Approximation.\n\n'); OK = FALSE; while OK == FALSE fprintf(1,'Input m and n on separate lines.\n'); LM = input(' '); LN = input(' '); BN = LM+LN; if LM >= 0 & LN >= 0 OK = TRUE; else fprintf(1,'m and n must both be nonnegative.\n'); end; if LM == 0 & LN == 0 OK = FALSE; fprintf(1,'Not both m and n can be zero\n'); end; end; OK = FALSE; while OK == FALSE fprintf(1,'The MacLaurin coefficients a(0), a(1), ... , a(N)\n'); fprintf(1,'are to be input.\n'); fprintf(1,'Choice of input method:\n'); fprintf(1,'1. Input entry by entry from keyboard\n'); fprintf(1,'2. Input data from a text file\n'); fprintf(1,'Choose 1 or 2 please\n'); FLAG = input(' '); if FLAG == 1 | FLAG == 2 OK = TRUE; end; end; AA = zeros(1,BN+1); if FLAG == 1 fprintf(1,'Input in order a(0) to a(N)\n');
for I = 0 : BN fprintf(1,'Input A( %d ) \n',I); AA(I+1) = input(' '); end; end; if FLAG == 2 fprintf(1,'As many entries as desired can be placed\n'); fprintf(1,'on each line of the file each separated by blank.\n'); fprintf(1,'Has such a text file been created?\n'); fprintf(1,'Enter Y or N\n'); AAA = input(' ','s'); if AAA == 'Y' | AAA == 'y' fprintf(1,'Input the file name in the form - '); fprintf(1,'drive:\name.ext\n'); fprintf(1,'for example: A:\DATA.DTA\n'); NAME = input(' ','s'); INP = fopen(NAME,'rt'); for I = 0 : BN AA(I+1) = fscanf(INP, '%f',1); end; fclose(INP); else fprintf(1,'Please create the input file.\n'); fprintf(1,'The program will end so the input file can '); fprintf(1,'be created.\n'); OK = FALSE; end; end; if OK == TRUE % STEP 1 N = BN; M = N+1; % STEP 2 - performed in input NROW = zeros(1,N); for I = 1 : N NROW(I) = I; end; % initialize row pointer for linear system NN = N-1; % STEP 3 Q = zeros(1, LM + 1); P = zeros(1, LN + 1); A = zeros(N,N+1); Q(1) = 1; P(1) = AA(1); % STEP 4 % Set up a linear system, but use A(i,j) instead of B(i,j). for I = 1 : N % STEP 5 for J = 1 : I- if J <= LN A(I,J) = 0; end; end;
XM = A(J1,I)/A(I1,I);
% STEP 17
for K = JJ : M A(J1,K) = A(J1,K)-XM * A(I1,K); end; % STEP 18 A(J1,I) = 0; end; end; I = I+1; end; if OK == TRUE % STEP 19 N1 = NROW(N); if abs(A(N1,N)) <= 1.0e- OK = FALSE; % system has no unique solution else % STEP 20 % Start backward substitution. if LM > 0 Q(LM+1) = A(N1,M)/A(N1,N); A(N1,M) = Q(LM+1); end; PP = 1; % STEP 21 for K = LN+1 : NN I = NN-K+LN+1; JJ = I+1; N2 = NROW(I); SUM = A(N2,N+1); for KK = JJ : N LL = NROW(KK); SUM = SUM-A(N2,KK)A(LL,M); end; A(N2,M) = SUM/A(N2,I); Q(LM-PP+1) = A(N2,M); PP = PP+1; end; % STEP 22 for K = 1 : LN I = LN-K+1; N2 = NROW(I); SUM = A(N2,N+1); for KK = LN+1 : N LL = NROW(KK); SUM = SUM-A(N2,KK)A(LL,M); end; A(N2,M) = SUM; P(LN-K+2) = A(N2,M); end; % STEP 23 % procedure completed successfully fprintf(1,'Choice of output method:\n');
fprintf(1,'1. Output to screen\n'); fprintf(1,'2. Output to text file\n'); fprintf(1,'Enter 1 or 2\n'); FLAG = input(' '); if FLAG == 2 fprintf(1,'Input the file name in the form - drive:\name.ext\n'); fprintf(1,'for example: A:\OUTPUT.DTA\n'); NAME = input(' ','s'); OUP = fopen(NAME,'wt'); else OUP = 1; end; fprintf(OUP, 'PADE RATIONAL APPROXIMATION\n\n'); fprintf(OUP, 'Denominator Coefficients Q(0), ..., Q(M) \n'); for I = 0 : LM fprintf(OUP, ' %11.8f', Q(I+1)); end; fprintf(OUP, '\n'); fprintf(OUP, 'Numerator Coefficients P(0), ..., P(N)\n'); for I = 0 : LN fprintf(OUP, ' %11.8f', P(I+1)); end; fprintf(OUP, '\n'); if OUP ~= 1 fclose(OUP); fprintf(1,'Output file %s created successfully \n',NAME); end; end; end; if OK == FALSE fprintf(1,'System has no unique solution\n'); end; end;
