Padé approximant for a numerical data array

20 次查看（过去 30 天）
Saeid
Saeid 2020-7-16
评论： John D'Errico 2020-7-16
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?

请先登录，再回答此问题。

采纳的回答

KSSV
KSSV 2020-7-16
syms x ;
pade(atanh(x))
  2 个评论
Saeid
Saeid 2020-7-16
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?

请先登录，再进行评论。

更多回答（1 个）

John D'Errico
John D'Errico 2020-7-16
编辑：John D'Errico 2020-7-16
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 个评论
Saeid
Saeid 2020-7-16
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 2020-7-16
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Curve Fitting Toolbox 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by