Padé approximant for a numerical data array

20 Ansichten (letzte 30 Tage)
Saeid
Saeid am 16 Jul. 2020
Kommentiert: John D'Errico am 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?

Akzeptierte Antwort

KSSV
KSSV am 16 Jul. 2020
syms x ;
pade(atanh(x))
  2 Kommentare
Saeid
Saeid am 16 Jul. 2020
Hi, thanks for the reply. I guess my question has not been formulated properly. How about when y is not symbolically defined but just numerically in the form of an array?
KSSV
KSSV am 16 Jul. 2020
Pade supports symbolic variable.

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (1)

John D'Errico
John D'Errico am 16 Jul. 2020
Bearbeitet: John D'Errico am 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 Kommentare
Saeid
Saeid am 16 Jul. 2020
Thanks John, this will still help me a lot, since fortunately I have the curve fitting toolbox! And now I need to get to work on the real (x,y) dataset, which is just a bit more complicated than this, but whose behavior qualitatively resembles that of tanh(x).
John D'Errico
John D'Errico am 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.

Melden Sie sich an, um zu kommentieren.

Community Treasure Hunt

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

Start Hunting!

Translated by