Page 142 - Applied Numerical Methods Using MATLAB

PADE APPROXIMATION BY RATIONAL FUNCTION 131

function [num,den] = padeap(f,xo,M,N,x0,xf)
%Input:f= function to be approximated around in [xo, xf]
%Output: num = numerator coeffs of Pade approximation of degree M
% den = denominator coeffs of Pade approximation of degree N
a(1) = feval(f,xo);
h = .01; tmp = 1;
fori=1:M+N
tmp = tmp*i*h; %i!h^i
dix = difapx(i,[-i i])*feval(f,xo+[-i:i]*h)’; %derivative(Section 5.3)
a(i + 1) = dix/tmp; %Taylor series coefficient
end
for m = 1:N
n = 1:N; A(m,n) = a(M +1+m-n);
b(m)=-a(M+1+m);
end
d=A\b’; %Eq.(3.4.4b)
form=1:M+1
mm = min(m - 1,N);
q(m) = a(m:-1:m - mm)*[1; d(1:mm)]; %Eq.(3.4.4a)
end
num = q(M + 1:-1:1)/d(N); den = [d(N:-1:1)’ 1]/d(N); %descending order
if nargout == 0 % plot the true ftn, Pade ftn and Taylor expansion
if nargin < 6, x0 = xo - 1; xf = xo + 1; end
x = x0+[xf-x0]/100*[0:100]; yt = feval(f,x);
x1 = x-xo; yp = polyval(num,x1)./polyval(den,x1);
yT = polyval(a(M + N + 1:-1:1),x1);
clf, plot(x,yt,’k’, x,yp,’r’, x,yT,’b’)
end

x
Example 3.2. Pade Approximation for f(x) = e . Let’s find the Pade approx-
o
x
imation p 3,2 (x) = Q 3 (x)/D 2 (x) for f(x) = e around x = 0. We make the
MATLAB program “do_pade.m”, which uses the routine “padeap()” for this
job and uses it again with no output argument to see the graphic results as
depicted in Fig. 3.6.
>>do_pade %Pade approximation
n = 0.3333 2.9996 11.9994 19.9988
d = 1.0000 -7.9997 19.9988
%do_pade.m to get the Pade approximation for f(x) = e^x
f1 = inline(’exp(x)’,’x’);
M=3;N=2; %the degrees of Numerator Q(x) and Denominator D(x)
xo = 0; %the center of Taylor series expansion
[n,d] = padeap(f1,xo,M,N) %to get the coefficients of Q(x)/P(x)
x0 = -3.5; xf = 0.5; %left/right boundary of the interval
padeap(f1,xo,M,N,x0,xf) %to see the graphic results

To confirm and support this result from the analytical point of view and to help
the readers understand the internal mechanism, we perform the hand-calculation