Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Spring 2013
Lecture topic: Pade Approximants

Initializations

$TextStyle={FontFamily->"Times",FontSize->16};
Off[General::"spell1"]

Taylor Series Approximation

Here is a known function:

Clear[g]
g[x_]=Exp[x]/(1+x^2)

plot1=Plot[g[x],{x,-3,3},PlotStyle->{Red,Thick},ImageSize->{500,Automatic},AxesLabel->{"x","y"}]

Here is the Taylor series approximation to a function (using our own code to generate the Taylor series):

PT[f_,a_,n_,x_] :=(∑_(j=0)^n)((D[f,{x,j}]/.x->a)/(j!))(((x-a))^j)

taylor[x_]=PT[g[x],1/2,12,x]

plot2=Plot[taylor[x],{x,-3,3},PlotStyle->{Blue,Thick},ImageSize->{500,Automatic},AxesLabel->{"x","y"}]

Show[plot2,plot1,AxesLabel->{"x","y"},PlotRange->{All,{0,2}}]

We can show radius of convergence is ρ=((Sqrt[5])/2)by examining the complex poles of the original function.

Pade Approximant

We can use the Taylor series to construct another approminant based on rational functions.

Clear[p,q]
f=g[x];
a=1/2;

m=2;
n=3;
P=(∑_(i=0)^m)(Subscript[p, i])(((x-a))^i);
Q=1+(∑_(i=1)^n)(Subscript[q, i])(((x-a))^i);

expression1=P+(PT[f,a,m+n+1,x])*Q;

unknowns={};
Do[unknowns=Append[unknowns,(Subscript[p, i])],{i,0,m}];
Do[unknowns=Append[unknowns,(Subscript[q, i])],{i,1,n}];
unknowns

Collect powers of z=x-a

expression2=Collect[expression1/.x->z+a,z]

Set coefficients of powers of z=x-a to zero (we need m+n+1 equations):

equations={};
Do[equations=Append[equations,Coefficient[expression2,z,i]==0],{i,0,m+n}];
equations

Solve the system of equations:

sol=Solve[equations,unknowns]

Now we form the rational function approximation:

pade[x_]=-(P/Q)/.Flatten[sol]

plot3=Plot[pade[x],{x,-3,3},PlotStyle->{Green,Thick},ImageSize->{500,Automatic},AxesLabel->{"x","y"}]

Show[plot3,plot2,plot1,AxesLabel->{"x","y"},PlotRange->{All,{0,2}}]

The Power of Pade Approximants

Pade Approximants are extremely useful, since there are many applications (energy of the quantum anharmonic oscillator is one example) where you have access to the Taylor series, but not the underlying function it represents. The Taylor series may converge slowly, or not at all. Yet Pade Approximants allow you to pull information out of the Taylor series and reconstruct the underlying function. If the underlying function is a rational function, the Pade Approximants will get the exact function for you!

Pade Approximants approximates the function well

g[x_]=Exp[x]/(1+x^2)

The Taylor series had a radius of convergence of ρ= (Sqrt[5])/2. Let's include that in the plot

plot1=Plot[g[x],{x,-3,3},PlotStyle->{Red,Thick},ImageSize->{500,Automatic},AxesLabel->{"x","y"},GridLines->{{1/2-(Sqrt[5])/2,1/2+(Sqrt[5])/2},None}]

taylor[x_]=PT[g[x],1/2,100,x];

plot2=Plot[taylor[x],{x,-3,3},PlotStyle->{Blue,Thick},ImageSize->{500,Automatic},AxesLabel->{"x","y"},PlotRange->{{-3,3},{0,2}}];

Show[plot1,plot2,AxesLabel->{"x","y"},PlotRange->{All,{0,2}}]

Here is the error:

ploterr1=Plot[Abs[taylor[x]-g[x]],{x,-6,6},GridLines->{{1/2-(Sqrt[5])/2,1/2+(Sqrt[5])/2},None},PlotStyle->{Orange,Thick},ImageSize->{500,Automatic},AxesLabel->{"x","y"},PlotRange->{{-3,3},{0,2}}]

No matter how many terms we keep in the Taylor series, we are not going to approximate the function well outside the interval of convergence.

Here is the Pade code, cleaned up to surpress output:

f=g[x];

(*define the point we are expanding about and the degree of the polynomials*)
a=1/2;

m=8;
n=8;

(*define the rational function numerator P and denominator Q*)
P=(∑_(i=0)^m)(Subscript[p, i])(((x-a))^i);
Q=1+(∑_(i=1)^n)(Subscript[q, i])(((x-a))^i);

(*define the espression that we will expand to collect powers of z+a in*)
expression1=P+(PT[f,a,m+n+1,x])*Q;

(*set up a list of all the unknown coefficients*)
unknowns={};
Do[unknowns=Append[unknowns,(Subscript[p, i])],{i,0,m}];
Do[unknowns=Append[unknowns,(Subscript[q, i])],{i,1,n}];

(*collect powers of z+a*)
expression2=Collect[expression1/.x->z+a,z];

(*pull of the coefficients to get the set of equations to solve*)
equations={};
Do[equations=Append[equations,Coefficient[expression2,z,i]==0],{i,0,m+n}];

(*solve the system of equations for the coefficients*)
sol=Solve[equations,unknowns];

(*put everything back into the Pade approximant*)
pade[x_]=-(P/Q)/.Flatten[sol];

Here is a plot of the Pade approximant [8,8] to the function

plot3=Plot[pade[x],{x,-3,3},PlotStyle->{Green,Thick},ImageSize->{500,Automatic},AxesLabel->{"x","y"}];
Show[plot3,plot1,AxesLabel->{"x","y"},PlotRange->{All,{0,2}}]

Looks pretty good! You really can't tell there are two functions being plotted.

Let's look at the error:

ploterr2=Plot[Abs[pade[x]-g[x]],{x,-6,6},GridLines->{{1/2-(Sqrt[5])/2,1/2+(Sqrt[5])/2},None},PlotStyle->{Blue,Thick},PlotRange->{All,All},ImageSize->{500,Automatic},AxesLabel->{"x","y"}]

This approximates the function extremely well far outside the interval of convergence of the Taylor series!

Show[ploterr2,ploterr1,PlotRange->{All,{0,0.0002}}]

Pade Approximants can assign a sum to a divergent series

Let's look at the Taylor series about x=0 for Log(x+1). We'll use Mathematica commands for the Taylor series to speed things up a bit.

g[x_]=Log[x+1]

Notice the Taylor series converges nicely at x=1/2:

g[1/2]
Normal[Series[g[x],{x,0,12}]]
Sum[(-1)^(n+1)(1/2)^n/(n),{n,1,Infinity}]

Notice the Taylor series diverges at x=2:

g[2]
Normal[Series[g[x],{x,0,12}]]
Sum[(-1)^(n+1)(2)^n/(n),{n,1,Infinity}]

Pade can recover the Log[3] result for this divergent series! Notice m and n are relatively small.

f=g[x];

a=0;

m=8;
n=8;
P=(∑_(i=0)^m)(Subscript[p, i])(((x-a))^i);
Q=1+(∑_(i=1)^n)(Subscript[q, i])(((x-a))^i);

expression1=P+(PT[f,a,m+n+1,x])*Q;

unknowns={};
Do[unknowns=Append[unknowns,(Subscript[p, i])],{i,0,m}];
Do[unknowns=Append[unknowns,(Subscript[q, i])],{i,1,n}];

expression2=Collect[expression1/.x->z+a,z];

equations={};
Do[equations=Append[equations,Coefficient[expression2,z,i]==0],{i,0,m+n}];

sol=Solve[equations,unknowns];

pade[x_]=-(P/Q)/.Flatten[sol];

SetPrecision[pade[2],12]

SetPrecision[Log[3.0],12]

And that is impressive indeed.

Pade Approximants can "approximate" a rational function exactly!

Consider the rational function:

g[x_]=((1-x)/(((1+x))^2))

Notice the Taylor series is an infinite series:

Normal[Series[g[x],{x,0,20}]]

Let's see what the Pade result is!

f=g[x];

a=0;

m=4;
n=6;
P=(∑_(i=0)^m)(Subscript[p, i])(((x-a))^i);
Q=1+(∑_(i=1)^n)(Subscript[q, i])(((x-a))^i);

expression1=P+(PT[f,a,m+n+1,x])*Q;

unknowns={};
Do[unknowns=Append[unknowns,(Subscript[p, i])],{i,0,m}];
Do[unknowns=Append[unknowns,(Subscript[q, i])],{i,1,n}];

expression2=Collect[expression1/.x->z+a,z];

equations={};
Do[equations=Append[equations,Coefficient[expression2,z,i]==0],{i,0,m+n}];

sol=Solve[equations,unknowns];

pade[x_]=-(P/Q)/.Flatten[sol]

OK, this looks bad, but it really isn't. What is happening is that the system of equations has some slack in it, and some variables are left unspecified. That means these variables are arbitrary, and we can set them to be anything we want. So let's algebraically simplify the result (remember, those variables in the solution can be anything we want!) and see what happens.

Simplify[%]

Impressive again. We've recovered the original function. The fact this happened is not at all surprising if you think about it.

Mathematica Built-in Command: Pade

Mathematica has an intrinsic command to implement Pade Approximants. You must load a package first:

<<"FunctionApproximations`"

g[x_]=((Cos[x])/(Sqrt[1+x]));
pade[x_]=PadeApproximant[g[x],{x,0,{15,18}}]

ploterr1=Plot[{g[x],pade[x]},{x,-1,12},PlotStyle->{{Blue,Thick},{Red,Thick}},PlotRange->{All,{-1,1}},ImageSize->{500,Automatic},AxesLabel->{"x","y"}]

Using Wynn's Identity to Construct a Pade Table to Accelerate Convergence of a Sequence

We can set up Wynn's Identity as a simple set of recursion relations, which will fill out a Pade table and improve the convergence of a sequence. Let's get a sequence to work with using Newton's method.

Clear[f]

f[x_] =Exp[6x]+3Log[2]^2Exp[2x]-Log[8]Exp[4x]-Log[2]^3
Plot[f[x],{x,-0.3,-0.1},PlotStyle->Blue,ImageSize->{500,Automatic}]

Note that the zero here does not appear to be a simple zero, so we might not acheive quadratic convergence with Newton's method.

Newton's Method

Here is Newton’s method to generate a sequence whose limit is the root of f(x).

Clear[p0,p1,tol, max, i,p,q]

f[x_] =Exp[6x]+3Log[2]^2Exp[2x]-Log[8]Exp[4x]-Log[2]^3
p0=0.0;
tol=10.0^(-12);
max=300;

data={};
i=1;
Print["Newton's Method"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[i≤max ,
p=p0-f[p0]/f'[p0];
data=Append[data,{i,p,Abs[p-p0]}];
If[Abs[f[p]]≤tol,{Print["Result: ","p= ",p," Error ",Abs[p-p0], " Steps= ",i],Break[]}];
p0=p;
i=i+1]
If[i>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"n","p","Abs[p-p0]"}}]

Error in Newton's Method

error=Table[{n,Abs[data[[n,3]]]},{n,1,Length[data]}];
nmerrplot=ListLogPlot[error,PlotStyle->{Orange,PointSize[0.02]},AxesLabel->{"n","Logarithm of error"}]

Using Wynn's Identity to Create a Pade Table

Here we use Wynn's identity to create a Pade Table, where the initial column is the original sequence. Wynn's identity was solved for [L+1,M], and then we relabelled L+1->L.

The sequence is stored in data:

Length[data]

Table[data[[l,2]],{l,1,Length[data]}]

Wynn's identity was solved for [L+1,M], and then we relabelled L+1->L, which makes it look a bit messy.

Clear[MyPade,a];

(*Initializations*)
MyPade[l_,0]:=MyPade[l,0]=data[[l+1,2]]
MyPade[l_,-1]:=MyPade[l,-1]=a;
MyPade[-1,0]:=MyPade[-1,0]=0;
MyPade[-1,m_]:=MyPade[-1,m]=0;

(*Wynn'sIdentity - note I use the Limit to deal with the bit that is infinite*)
MyPade[l_,m_]:=MyPade[l,m]=Simplify[Limit[(MyPade[-1+l,-1+m] ((MyPade[l,-1+m])^2)-MyPade[l,-2+m] ((MyPade[l,-1+m])^2)+MyPade[-1+l,-1+m] MyPade[l,-2+m] MyPade[1+l,-1+m]-2 MyPade[-1+l,-1+m] MyPade[l,-1+m] MyPade[1+l,-1+m]+((MyPade[l,-1+m])^2) MyPade[1+l,-1+m])/(MyPade[-1+l,-1+m] MyPade[l,-2+m]-2 MyPade[l,-2+m] MyPade[l,-1+m]+((MyPade[l,-1+m])^2)-MyPade[-1+l,-1+m] MyPade[1+l,-1+m]+MyPade[l,-2+m] MyPade[1+l,-1+m]),a->Infinity]];

Here is the first column (the initial data we got from Newton’s Method):

Table[{l,MyPade[l,0]},{l,0,Length[data]-1}]//TableForm

Here is the second column (should be the same as Aitken's method, so am going to implement Aitken’s method to compare):

TableForm[Table[{l,MyPade[l,1],MyPade[l-1,0]-((((MyPade[l,0]-MyPade[l-1,0]))^2)/(MyPade[l+1,0]-2MyPade[l,0]+MyPade[l-1,0]))},{l,0,Length[data]-1-1}],TableHeadings->{None,{"n","Pade[L,1]","Aitken's (Δ^2)"}}]

Notice that for the second column I had to stop at 19 rather than 20 (also notice we atarted at 0 rather than 1). This is due to how we are implementing Wynn’s identity (see Figure 3.4.1 Compass Points). Let's create the Pade Table!

final=7;
TableForm[Table[SetPrecision[MyPade[l,m],7],{l,0,final},{m,0,final}],TableSpacing->{1,1},TableHeadings->{{"L"},{"M"}}]

Although there are convergence results only for specific types of functions, it is common for the Pade Approximant to converge as L,M-> ∞ along the diagonal of the Pade Table. Let's examine the error along the diagonal.

Error in Diagonal of Pade Table

error=Table[{n,Abs[Abs[MyPade[n,n]-MyPade[n-1,n-1]]]},{n,1,Floor[((Length[data])/2)]}];
pterrplot=ListLogPlot[error,PlotStyle->{Red,PointSize[0.02`]},AxesLabel->{"n","Logarithm of error"}]

Comparison of Errors (Newton and Pade)

Show[nmerrplot,pterrplot,PlotRange->{All,All},ImageSize->{500,Automatic}]

So the idea we have exploited here is to use Newton's method to get the first 21 elements in a sequence, and the from the 21 elements a Pade Table is constructed that refines the accuracy of the root.