Padé approximant
Summary
Padé approximant is a rationale function approximation expanded from a function value and the derivative function values similar to Taylor Series. Often, Padé approximant’s can converge when Taylor Series does not even though both use the same information at a point. The key is Taylor Series is a polynomial while Padé approximant is a rationale function. Rationale functions can have singularities/vertical-asymptotes and horizontal asymptotes. Polynomials do not have asymptotes or singularities.
High-order Padé Approximants can be hard to compute because of how fast the coefficients grow. Total orders (denominator plus numerator order) of less than 20 fail at times for MATLAB’s symbolic pade function in the code below. padeapprox from chebfun also struggles. It takes a lot of effort to compute a high-order Padé approximate.
The most useful way to understand Padé approximant is to watch the animations below. For details on the construction, see the resources section below.
Resources
Wikipedia Padé approximate - A good overview of Pade’s approximate.
Matlab Symbolic Pade Approximant Page
Chebfun has a function called padeapprox that calculates a Padé approximant numerically. padeapprox is based on this paper, Robust Padé Approximation via SVD.
Animations
\(exp(x)\) about 0 and 1
\(tanh(x)\) about 0 and 0.5
\(arctan(x)\) about 0 and 0.5
\(1/x\) about 0.5
\(1/(x-1)\) about 0.75
\(sin(x)\) about 0 and 0.5
\(x/(x^2 + 0.01^2)^\frac{1}{4}\) about 0.5
Runge function, \(1/(1+25x^2)\) about 0 and 0.5
Code to Create Animations
MATLAB 2020a was used to write and execute the code below.
clc; clear; close all;
%% exp(x) about 0
titleName = "exp(x) about 0";
func = @(x) exp(x);
expansionPoint = 0;
domain = [-5 5];
maxOrder = 20;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% exp(x) about 1
titleName = "exp(x) about 1";
func = @(x) exp(x);
expansionPoint = 1;
domain = [-5 5];
maxOrder = 20;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% tanh(x) about 0
titleName = "tanh(x) about 0";
func = @(x) tanh(x);
expansionPoint = 0;
domain = [-5 5];
maxOrder = 25;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% tanh(x) about 0.5
titleName = "tanh(x) about 0.5";
func = @(x) tanh(x);
expansionPoint = 0.5;
domain = [-5 5];
maxOrder = 30;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% arctan(x) about 0
titleName = "arctan(x) about 0";
func = @(x) atan(x);
expansionPoint = 0;
domain = [-5 5];
maxOrder = 30;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% arctan(x) about 0.5
titleName = "arctan(x) about 0.5";
func = @(x) atan(x);
expansionPoint = 0.5;
domain = [-5 5];
maxOrder = 30;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% 1/x about 0.5
titleName = "1/x about 0.5";
func = @(x) 1./x;
expansionPoint = 0.5;
domain = [-5 5];
maxOrder = 30;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% 1/(x-1) about 0.75
titleName = "1/(x-1) about 0.75";
func = @(x) 1./(x-1);
expansionPoint = 0.75;
domain = [-5 5];
maxOrder = 30;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% sin(x) about 0
titleName = "sin(x) about 0";
func = @(x) sin(x);
expansionPoint = 0;
domain = [-5 5];
maxOrder = 30;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% sin(x) about 0.5
titleName = "sin(x) about 0.5";
func = @(x) sin(x);
expansionPoint = 0.5;
domain = [-5 5];
maxOrder = 30;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% x/(x^2 + 0.01^2)^(1/4)
titleName = "x/(x^2 + 0.01^2)^(1/4) about 0.5";
func = @(x) x./(x.^2 + 0.01^2).^(1/4);
expansionPoint = 0.5;
domain = [-5 5];
maxOrder = 30;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% Runge function
titleName = "Runge function, 1/(1+25*x^2) about 0";
func = @(x) 1./(1+25*x.^2);
expansionPoint = 0;
domain = [-5 5];
maxOrder = 20;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%% Runge function
titleName = "Runge function, 1/(1+25*x^2) about 0.5";
func = @(x) 1./(1+25*x.^2);
expansionPoint = 0.5;
domain = [-5 5];
maxOrder = 20;
createVideo(titleName, func, expansionPoint, domain, maxOrder);
%%
function createVideo(titleName,func,expansionPoint,domain,maxOrder,approach)
arguments
titleName (1,1) string;
func (1,1) function_handle;
expansionPoint (1,1) double;
domain (1,2) double;
maxOrder (1,1);
approach (1,1) string {mustBeMember(approach,["chebfun","sym"])} = "sym";
end
% setup figure
width = 1920;
height = 1080;
figureObject = figure("Position",[10 10 width height]);
% setup video file
videoName = regexprep(titleName,["/","\*"],"_") + ".mp4";
video = VideoWriter(videoName,"MPEG-4");
video.FrameRate = 2;
try
video.open();
% setup function to approximate
fplot(func,domain,"linewidth",3,"DisplayName",func2str(func))
hold all;
xlabel("x"); ylabel("y")
title(titleName)
ylim(ylim);
xlim(xlim);
% find the permutations of the numerator and denominator
[numeratorOrder, denominatorOrder] = ndgrid(0:maxOrder,0:maxOrder);
totalOrder = numeratorOrder + denominatorOrder;
isOrderToHigh = totalOrder > maxOrder;
numeratorOrder(isOrderToHigh) = [];
denominatorOrder(isOrderToHigh) = [];
[~,sortedIndex] = sortrows([denominatorOrder(:)+numeratorOrder(:), numeratorOrder(:), denominatorOrder(:)]);
numeratorOrder = numeratorOrder(sortedIndex);
denominatorOrder = denominatorOrder(sortedIndex);
% loop
for i=1:length(numeratorOrder)
% calculate next pade approximant
order = [numeratorOrder(i), denominatorOrder(i)];
switch approach
case "chebfun"
% padeapprox is used from the Chebfun toolbox by Trefethen and others. The advantage is computation is
% faster then symbolic computation. The disadvantage is has numerical issues at total order
% (numeratorOrder + denominatorOrder) of 16 or so.
%
% http://www.chebfun.org/
% P. Gonnet, S. Guettel, and L. N. Trefethen, "Robust Pade approximation via SVD", SIAM Rev., 55:101-117, 2013.
%
% Commentary: numerical issues are more of a statement of ill-conditioning of monomial basis combined
% with constructing a pade approximate. Finding a pade approximate may be possible via an orthogonal
% polynomial basis but I am not aware of any software that does this.
%
padeApproximant = padeapprox(@(x) func(x+expansionPoint), order(1),order(2));
padeApproximant = @(x) padeApproximant(x-expansionPoint);
case "sym"
% The symbolic approach is accurate but slow. For high enough total order (numeratorOrder +
% denominatorOrder), large integer errors occur stopping further computation. To avoid this as much as
% possible, we try to tell the symbolic engine that we want variable precision arithmetic, vpa().
if i==1
syms x real
digits(20);
func = vpa(func(vpa(1.0)*x));
end
padeApproximant = pade(func, x, expansionPoint, 'order', order);
end
label = sprintf("Padé approximant, Pade Order = [%3d,%3d],Total Order = %3d",order,sum(order));
if i ==1
% setup plot
fplotObject = fplot(padeApproximant,"--","LineWidth",3,"DisplayName",label);
legend("location","best");
else
% plot
fplotObject.Function = padeApproximant;
fplotObject.DisplayName = label;
end
% write frame
video.writeVideo(getframe(figureObject));
end
video.close();
catch e %#ok
video.close();
fprintf("%s around %g failed at pade order = [%d, %d]\n",string(func),expansionPoint,order(1),order(2));
% if exist(videoName,"file")
% delete(videoName);
% end
rethrow(e);
end
end