G Rational Padé Approximations

  The (p,q) Padé approximant to the matrix exponential e^X is, by definition, the unique (p,q) rational function [100]
$$\begin{gather}R_{pq}(\boldsymbol{X})=\frac{N_{pq}(\boldsymbol{X})}{D_{pq}(\boldsymbol{X})}, \end{gather}$$
which matches the Taylor series expansion of e^X through terms to the power p+q. Its coefficients are therefore determined by solving the algebraic equations
$$\begin{gather}e^{\boldsymbol{X}}=\sum_{j=0}^{\infty}\frac{\boldsymbol{X}^j}{j!}=... ...{X})}{D_{pq}(\boldsymbol{X})}+O\left(\boldsymbol{X}^{p+q+1}\right). \end{gather}$$
The result is
$$\begin{gather}N_{pq}(\boldsymbol{X})=\sum_{j=0}^p\frac{(p+q-j)!q!}{(p+q)!j!(p-j)... ...\sum_{j=0}^q\frac{(p+q-j)!p!}{(p+q)!j!(q-j)!} (\boldsymbol{-X})^j. \end{gather}$$
Choosing p=q one obtains the diagonal Padé approximation. This choice is to prefer, because it yields a higher order approximation with the same amount of computation.
$$\begin{gather}R_{pp}=\frac{N_{pp}(\boldsymbol{X})}{N_{pp}(\boldsymbol{-X})}\qqua... ...j!(p-j)!}\boldsymbol{X}^j \equiv \sum_{j=0}^p c_j\boldsymbol{X}^j, \end{gather}$$
where the coefficients c_j can be conveniently constructed by means of the recursion
$$\begin{gather}c_0=1,\qquad c_j=c_{j-1}\frac{p+1-j}{j(2p+1-j)}. \end{gather}$$
The computation of the polynomials are best done with a Horner scheme. C. Moler and C. Van Loan [87] showed that if $\Vert\boldsymbol{X}\Vert _2/2^m\leq 1/2$, then
$$\begin{gather}\left[R_{pp}\left(\frac{\boldsymbol{X}}{2^m}\right)\right]^{2^m}=e... ...^{2p-3}\frac{(p!)^2}{(2p)!(2p+1)!}\approx 0.34\cdot 10^{-15} (p=6) \end{gather}$$