On the Jacobi formula for Bivariate Pade Approximants of Rectangular Type

In this paper a recursive algorithm is presented for evaluating multivariate Padé approximants (of the rectangular type described in the work of Lutterodt) which is analogous to the Jacobi formula for univariate Padé approximants. This algorithm is then applied to a (singular) Riccati differential equation to generate fast and accurate approximate solutions.

💡 Research Summary

The paper introduces a recursive algorithm for constructing rectangular‑type multivariate Padé approximants, extending the classic Jacobi recurrence known for univariate Padé approximants. Starting from a bivariate power series f(x,y)=∑{n,m≥0}c{n,m}x^n y^m with c_{0,0}=0, the authors define left‑hand approximants f_L^{n,m}(x,y)=A_{n,m}(x,y)/B_{n,m}(x,y) and right‑hand approximants f_R^{n,m}(x,y) by swapping the variables. The numerator and denominator are expressed as collections of polynomials A_{n,m}^{(p)}(x) and B_{n,m}^{(p)}(x) of degree n for each y‑degree p (0≤p≤m). The index sets follow Lutterodt’s rectangular scheme: the same set N_{n,m}=D_{n,m}={ (i,j): i≤n, j≤m } is used for both numerator and denominator, while the equations to be satisfied are drawn from Q_{n,m}={ (i,j): 0≤i≤2n, 0≤j≤m }∪{ (0,j): m+1≤j≤2m }.
To compute the coefficients efficiently, the authors develop a double induction. First, for a fixed n, they treat the scalar parameters b_{0,1},…,b_{0,m} as unknowns (with b_{0,0}=1) and, assuming the polynomials B_{n,p} for p<m are already known, solve a linear system derived from the requirement that the error term be O(x^{2n+1}). This yields A_{n,m} and B_{n,m}. Second, they extend the classic Jacobi recurrence A_n=(1+β_n x)A_{n-1}+α_n x^2 A_{n-2}, B_n=(1+β_n x)B_{n-1}+α_n x^2 B_{n-2} by adding correction terms \hat A_{n,m} and \hat B_{n,m}. These corrections are linear combinations of previously computed polynomials for lower y‑degrees, with coefficients \hat α_{n,m,p} and \hat β_{n,m,p}. Crucially, the base case p=0 coefficients \hat α_{n,m,0} and \hat β_{n,m,0} are expressed in closed form using the univariate Jacobi coefficients F(n) (see equations 20‑21). Proposition 2.1 proves that the combined recurrence guarantees the error polynomial E_{n,m}=A_{n,m}−∑{k=0}^{m}B{n,k}C_{2n,k} remains O(x^{2n+1}), by expanding E_{n,m}, defining auxiliary σ and τ coefficients, and performing a term‑by‑term cancellation.
Implementation details include storing all A_{n,p} and B_{n,p} in matrices, which allows simultaneous updates for all p and reduces the overall computational complexity to roughly O(n·m). The correction matrices \hat α, \hat β are also pre‑computed, making the algorithm amenable to vectorised or parallel execution. For n=1 the correction coefficients are linear in the scalar parameters b_{0,1}…b_{0,m}, while for n>1 non‑linear dependence appears; nevertheless, in many practical cases this dependence can be eliminated.
The authors then apply the framework to a singular Riccati differential equation x w’ – β w + β w² + α x = 0, w(0)=w(1)=0, with positive constants α, β. By the substitution u=w’/w the equation reduces to a second‑order linear ODE whose exact solution is a ratio of Bessel functions: w(x)=z^{2β}