# Padé approximant for a numerical data array

25 views (last 30 days)

Show older comments

I am looking for the Padé approximant of a data array such as y. y cannot be given explicitly and is generated by another routine (e.g., solution to a differential equation)

x=[some data array];

y=[some data array];

pade(y,x,0,'Order',[3 3])

I receive an error message. I assume this is because Padé only works with symbolic math, but is there an alternative solution within matlab for this?

##### 0 Comments

### Accepted Answer

### More Answers (1)

John D'Errico
on 16 Jul 2020

Edited: John D'Errico
on 16 Jul 2020

Pade is only defined for a symbolic functions. If all you have is a list of (x,y) pairs, then you cannot use pade. And that appears to be your problem. Instead, you can use the curve fitting toolbox, which does offer a rational polynomial fit. That is all a pade approximant really is.

x = linspace(.01,0.99,99);

y = atanh(x);

Note that I excluded 0 and 1 as problematic points for a rational polynomial fit here.

ft = fittype('rat22')

ft =

General model Rat22:

ft(p1,p2,p3,q1,q2,x) = (p1*x^2 + p2*x + p3) / (x^2 + q1*x + q2)

[mdl,stuff] = fit(x',y',ft)

Warning: Start point not provided, choosing random start point.

> In curvefit.attention/Warning/throw (line 30)

In fit>iFit (line 307)

In fit (line 116)

mdl =

General model Rat22:

mdl(x) = (p1*x^2 + p2*x + p3) / (x^2 + q1*x + q2)

Coefficients (with 95% confidence bounds):

p1 = -5480 (-3.295e+06, 3.284e+06)

p2 = 6451 (-3.866e+06, 3.879e+06)

p3 = -70.87 (-4.265e+04, 4.251e+04)

q1 = -5772 (-3.469e+06, 3.458e+06)

q2 = 6080 (-3.643e+06, 3.655e+06)

stuff =

struct with fields:

sse: 0.0183300436668179

rsquare: 0.999383646351906

dfe: 94

adjrsquare: 0.999357418537093

rmse: 0.0139642566769813

plot(mdl),hold on,plot(x,y)

That is not unreasonable. An intelligent choice of starting values is possible. Be careful with high order rational polynomial fits ove a large domain. But a better result can come just from understanding the model you want to fit, and the character of the singularity, as well as where it lives.

The atanh function has a pole at x == 1. So build a model that reflects this essential behavior. And this time, I can even choose more intelligent starting points. The only one that matters is the one that locates the singularity.

ft = fittype('rat41')

ft =

General model Rat41:

ft(p1,p2,p3,p4,p5,q1,x) = (p1*x^4 + p2*x^3 + p3*x^2 + p4*x + p5) /

(x + q1)

[mdl,stuff] = fit(x',y',ft,'start',[1 1 1 1 1 -1])

mdl =

General model Rat41:

mdl(x) = (p1*x^4 + p2*x^3 + p3*x^2 + p4*x + p5) /

(x + q1)

Coefficients (with 95% confidence bounds):

p1 = 0.8268 (0.7445, 0.9091)

p2 = -1.377 (-1.556, -1.198)

p3 = 1.619 (1.486, 1.753)

p4 = -1.156 (-1.194, -1.118)

p5 = 0.006645 (0.003171, 0.01012)

q1 = -1.025 (-1.026, -1.024)

stuff =

struct with fields:

sse: 0.00129250606291084

rsquare: 0.999956539065507

dfe: 93

adjrsquare: 0.999954202456126

rmse: 0.00372799069941909

plot(mdl),hold on,plot(x,y)

So the errors are cut considerably, just by recognizing we might be able to get away with a rational polynomial fit that is first order in the denominator. And since we know where the singularity lives, we should do pretty well.

The rational polynomial approximation has a singularity at x = -mdl.q1, so it found x = 1.025. That seems reasonable.

Is this a reasonable approximation? We might get a hint if the numerator polynomial also has a root also near x==1.

roots([mdl.p1,mdl.p2,mdl.p3,mdl.p4,mdl.p5])

ans =

0.301423835666943 + 1.10495965334969i

0.301423835666943 - 1.10495965334969i

1.05726192810239 + 0i

0.00579519514280635 + 0i

Some further understanding of the nature of the singularity at x == 1 would useful. We could gain that by looking at the function tanh itself. But then, I've already used a lot of space here, and while chases that wander down rabbit holes seem to be my specialty, I will only go so far.

##### 2 Comments

John D'Errico
on 16 Jul 2020

### See Also

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!