5

I try to find out for quite some time now, how Matlab implements the calculation of Padé Approximants using its symbolic pade function. (the code of is buried in a compiled mex-file)

I compared its output with the "direct calculation" via coefficient comparison and solving the resulting linear systems of equations. I used floating point, variable precision arithmetic and pure symbolic variables and even in the last case Matlab's implementation achieves far better results for orders $> 7$.

So there must be a more sophisticated algorithm involved, but the documentation remains silent. The only tiny bit of information I found was the "Euclidian Algorithm" - but in my understanding that does not directly lead to the Padé-Approximant, just to the greatest common divisor of two numbers, am I right?

I contacted the Matlab Support and received the answer, that if no source is given the algorithm would be their own development. I actually can't believe that they developed an algorithm from scratch, if there are books full of alternatives (e.g. Baker, see below), or do they? At least their algorithm should be based on something?

Baker, George A.jun.; Graves-Morris, Peter, Padé approximants., Encyclopedia of Mathematics and Its Applications. 59. Cambridge: Cambridge Univ. Press. xiv, 746 p. (1996). ZBL0923.41001.

Can you give me some hints on how they could have implemented the algorithm? What approach for calculating Padé approximants could be based on the mentioned "Euclidian Algorithm"? How could I proceed from there?

I actually could just use their implementation, but I feel for any publication one should know how it is done.

Robert Seifert
  • 285
  • You can check out the documentation of expm. Matlab's code for matrix exponential is based on padé-approximation. There are also expmdemo-files that should show how padé is used. Maybe there are some hints in there? – mdot Dec 31 '17 at 10:40
  • Did you check out https://en.wikipedia.org/wiki/Pad%C3%A9_approximant? There I explained once how to use the extended Euclidean algorithm to get a diagonal of the Padé table. I think I wrote it on the talk page and someone else transferred it with proper references to the article page. – Lutz Lehmann Dec 31 '17 at 10:47
  • @max I already read about the implementation of expm, I will have a deeper look in the algortihm. But I guess the case of $e^A$ is much simpler than the generic case, probably the reason why Matlab does not use a real Padé approximation within the code (it does not call its own pade function), but rather hard coded (!!!) Padé coefficients. – Robert Seifert Dec 31 '17 at 11:13

1 Answers1

4

Citing myself on https://en.wikipedia.org/wiki/Talk:Pad%C3%A9_approximant, one of the classical ways to compute Padé approximants uses the extended Euclidean algorithm. To compute the $[m/n]$ approximant, set \begin{align} &&r_0&=x^{m+n+1},&v_0&=0,\\ &&r_1&=f(x), &v_1&=1,\\ q_1&=r_0\text{ div }r_1,& r_2&=r_0-q_1r_1,& v_2&=v_0-q_1v_1\\ &\ \ \vdots\\ q_k&=r_{k-1}\text{ div }r_k,&r_{k+1}&=r_{k-1}-q_kr_k,&v_{k+1}&=v_{k-1}-q_kv_k, \\ &\ \ \vdots \end{align} where $f(x)$ is the Taylor polynomial of degree $m+n$. Then in every step $$ v_kf(x)=r_k\mod x^{m+n+1} $$ which means that $\frac{r_k}{v_k}$ is the $[\deg r_k/\deg v_k]$ Padé approximant. For most $f$ the quotients $q_k$ will all have degree $1$ so that $\deg v_k=k-1$ and $\deg r_k=m+n+1-k$ so that in step $k=n+1$ one obtains the $[m/n]$ approximant.

Lutz Lehmann
  • 131,652
  • That looks promising! I already read the according wikipedia article on that, but couldn't understand it. But your answer seems more straightforward. I will try to implement it after the holidays and come back to you. Happy New Year! – Robert Seifert Dec 31 '17 at 11:18
  • Sorry for my late reaction, my boss kept me doing other things. Unfortunately the results I get do not match with what I expect for a given example. So maybe you want to have a look into my follow-up question. Your answer here, however, matches the approach in a paper I found. – Robert Seifert Jan 25 '18 at 15:40