MATLAB Answers

Padé approximant for a numerical data array

10 views (last 30 days)
Saeid
Saeid on 16 Jul 2020
Commented: John D'Errico on 16 Jul 2020
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?

Accepted Answer

KSSV
KSSV on 16 Jul 2020
syms x ;
pade(atanh(x))
  2 Comments
KSSV
KSSV on 16 Jul 2020
Pade supports symbolic variable.

Sign in to comment.

More Answers (1)

John D'Errico
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
John D'Errico on 16 Jul 2020
Things are never as simple as we want them to be. :) It is good you have the curve fitting TB. I have found it to be tremendously valuable over the years.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!