Padé approximant for a numerical data array

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?

Answer 1

[KSSV]

syms x ;
pade(atanh(x))

[Saeid]

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]

Pade supports symbolic variable.

Answer 2

[John D'Errico]

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.

[Saeid]

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]

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.