pade.m

用MATLAB语言编写的pade近似程序

function [a,b,p,z] = pade(f,N);

% Compute the coefficients for an O(N)/O(N) Pade' 
% approximant to the symbolic function f(x) about x=0
% 
% Parameters:
%
% f - orginal function (using symbolic variable)
% N - highest order term in approximant
% y - point to expand about (number)
% a - coefficients of the numerator (decending)
% b - coefficients of the denominator (decending)
% p - poles of the approximant
% z - zeros of the approximant
%
% Note: Recommend restricting Pade' approximant to 
% bi-quintic or less
%
% Present version only expands around x=0
%
% Written by: Greg von Winckel 01/18/2004
% 

clear a b p z lhs rhs c;

N1=N+1; Nc=2*N+1;

% Symbolic computation of first 2N+1 terms of
% the taylor series
fser=taylor(f,0,Nc);

% Extract the expansion coefficients
c=sym2poly(fser)';

% Reverse the order (ascending powers)
c=c(Nc:-1:1);

lhs=eye(N+1);
lhs(2:N1,2:N1)=toeplitz(c(N1:Nc-1),c(N1:-1:2));
rhs=[1;-c(N1+1:Nc)];

% Find denominator coefficients
b=qmr(lhs,rhs,1e-6,20);

% Compute the numerator coefficients
lhs=toeplitz(c(1:N1),[c(1) zeros(1,N)]);
a=lhs*b;

% reverse order to descending powers
a=a(N1:-1:1);
b=b(N1:-1:1);

% Find poles and zeros
z=roots(a);
p=roots(b);

% Convert coefficients to fractions
a=rats(a); 
b=rats(b);