Direct algorithm

We can calculate \([p/q](x)\,\) using a direct method. Defining \(N(x)\) and \(D(x)\) as \begin{equation*} N(x)=\sum_{n=0}^{p} c_n\,x^n\,, \qquad D(x)=\sum_{n=0}^{q} b_n\,x^n\,, \end{equation*} we have the following equivalence \begin{equation*} f(x) - \frac{\sum_{n=0}^{p} c_n\,x^n}{\sum_{n=0}^{q} b_n\,x^n}=\mathcal{O}(x^{p+q+1}) \quad \Leftrightarrow \quad \left(\sum_{n=0}^{+\infty}a_n\,x^n\right) \left(\sum_{n=0}^{q} b_n\,x^n\right) - \sum_{n=0}^{p} c_n\,x^n = \mathcal{O}(x^{p+q+1})\,. \end{equation*}

Matching the coefficients of \(N(x)\) and \(f(x)\,D(x)\,\) leads to a system of linear equations

\[\begin{split}\left\{\begin{array}{c c} \begin{aligned} c_0&=a_0\, b_0\\ c_1&=a_1 \,b_0 + a_0\, b_1\\ & \ \ \ \ \ \ \ \ \ \ \ \vdots \ \ \ \ \ \ \ \ \vdots \\ c_p&=a_p\, b_0+ a_{p-1}\,b_1+\dots +a_{p-q}\,b_q\\ 0&=a_{p+1}\,b_0 + a_p\,b_1+\dots +a_{p-q+1}\,b_q\\ 0&=a_{p+2}\,b_0 + a_{p+1}\,b_1+\dots +a_{p-q+2}\,b_q\\ & \ \ \ \ \ \ \ \ \ \ \ \vdots \ \ \ \ \ \ \ \ \vdots \\ 0&=a_{p+q}\, b_0+ a_{p+q-1}\,b_1+\dots+ a_p\,b_q \,\,\\ \end{aligned} \end{array}\right.\end{split}\]

with \(a_n=0\,,\) \(\forall n <0\,.\) Selecting the last \(q\) equations from the above system, we obtain the following matrix system

\[\begin{split}\left(\begin{array}{c c c c} a_{p-q+1}&a_{p-q+2}&\dots \dots&a_p\\ a_{p-q+2}&a_{p-q+3}&\dots \dots&a_{p+1}\\ \vdots&\vdots&\ddots&\vdots\\ a_{p}&a_{p+1}&\dots \dots&a_{p+q-1}\\ \end{array}\right) \left(\begin{array}{c} b_q\\ b_{q-1}\\ \vdots\\ b_1\\ \end{array}\right) =\left(\begin{array}{c} -a_{p+1}\\ -a_{p+2}\\ \vdots\\ -a_{p+q} \end{array}\right)\,,\end{split}\]

which can be expressed in a compact form \begin{equation*}\ {\bf A}\,b=a \,. \end{equation*} When possible, there are several ways to solve the above system [1]. A classical approach is to use LU factorization, consisting of expressing \({\bf A}={\bf LU}\), it requires the order of \(q^3\) operations to solve the system. Without loss of generality, assume that is possible to solve

\[{\bf A}\,b=a\]

and the denominator coefficients \(b_0\,,b_1\,,b_2\,, \dots\,,b_q\,,\) are well defined. The remaining \(p+1\) equations

\[\begin{split}\begin{array}{c c} \begin{aligned} c_0&=a_0\,b_0\,\\ c_1&=a_1\,b_0 + a_0\, b_1\\ c_2&=a_2\,b_0 + a_1\, b_1 + a_0\,b_2\\ \ \ \ \ \ \ \ \ & \ \ \vdots\\ c_p&=a_p\, b_0+ a_{p-1}\,b_1+\dots +a_{p-q}\,b_q\\ \end{aligned} \end{array}\end{split}\]

gives the numerator coeficients \(c_0\,,c_1\,,c_2\,, \dots\,,c_p\,\) values, and the direct algorithm ends. It requires the order of \(q^3+pq\,\) operations to derive an approximant.

[1] Trefethen, L. N. and Bau, I. I. I. D. (1997). Numerical Linear Algebra. Society for Industrial and Applied Mathematics, Philadelphia.

Examples

import sympy as sp
from padepy import direct_algorithm as da

var = sp.Symbol("x")
func = sp.exp(var)
p = 3
q = 1
func

$\displaystyle e^{x}$

pade = da.pade(p, q, var, func)
pade

$\displaystyle \frac{\frac{x^{3}}{24} + \frac{x^{2}}{4} + \frac{3 x}{4} + 1}{1 - \frac{x}{4}}$

var = sp.Symbol("z")
func = (var + 1) / sp.sqrt(var**2 + 1)
p = 1
q = 7
func

$\displaystyle \frac{z + 1}{\sqrt{z^{2} + 1}}$

pade = da.pade(p, q, var, func)
pade

$\displaystyle \frac{\frac{179 z}{184} + 1}{- \frac{5 z^{7}}{128} + \frac{147 z^{6}}{1472} - \frac{55 z^{5}}{1472} - \frac{31 z^{4}}{368} - \frac{15 z^{3}}{368} + \frac{97 z^{2}}{184} - \frac{5 z}{184} + 1}$