Padé approximant for a numerical data array
11 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
Plus de réponses (1)
John D'Errico
le 16 Juil 2020
Modifié(e) : John D'Errico
le 16 Juil 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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/332708/image.jpeg)
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 commentaires
John D'Errico
le 16 Juil 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.
Voir également
Catégories
En savoir plus sur Linear Algebra dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!