Approximaths

Computing Padé approximants

Written by

approximaths

in

We have seen that Padé approximants offer a more efficient and flexible method for approximating functions than the Taylor expansion. They have been used in many areas of mathematics and physics.

We would like to find a more convenient way to calculate the Padé coefficients of any Taylor or Maclaurin series.

In the following post, we admit the existence of Padé approximants P(m,n). It is generally possible to find the coefficients of the Padé approximants except in the case of degeneracy due to particular values of the coefficients of the corresponding series. We will not consider these cases in the rest of our study. Remember that we defined Padé approximants as rational functions:

\displaystyle P(m,n) := \frac{A_0 + A_1 x + A_2 x^2 + \dots + A_m x^m}{1 + B_1 x + B_2 x^2 + \dots + B_n x^n}

If we want to solve a very hard or even impossible-to-solve-exactly problem using perturbation theory we will end up with a power series like \sum_{k=0}^{l} C_k x^k. To convert this (potentially not converging) series to its corresponding Padé approximant we write (as described in the previous post):

\displaystyle \frac{A_0 + A_1 x + \dots + A_m x^m}{1 + B_1 x + \dots + B_n x^n} = C_0 + C_1 x + \dots + C_{n+m}x^{n+m} + \mathcal{O}(x^{m+n+1})

From a purely algorithmic point of view and in order to calculate the Padé coefficients we drop \mathcal{O}(x^{m+n+1}):

\displaystyle (C_0 + C_1 x + \dots + C_{n+m}x^{n+m})(1 + B_1 x + \dots + B_n x^n) - (A_0 + A_1 x + \dots + A_m x^m) = 0
\displaystyle \sum_{i=0}^{n+m}\sum_{j=0}^{n}C_iB_jx^{i+j} - (A_0 + A_1 x + \dots + A_m x^m)= 0

Comparing each term of the geometric series on the right and the left of this equation (Uniqueness of power series coeffficients) we obtain the following set of equations:

\displaystyle (C_0 B_0 - A_0) = 0 \quad (i+j = 0)
\displaystyle (C_1 B_0 + C_0 B_1 - A_1)x = 0 \quad (i+j = 1)
\displaystyle (C_2 B_0 + C_1 B_1 + C_0 B_2 - A_2)x^2 = 0 \quad (i+j = 2)
\dots
\displaystyle (C_m B_0 + C_{m-1} B_1 + \dots + C_{m-n} B_n - A_m)x^m = 0 \quad (i+j = m)
\displaystyle (C_{m+1} B_0 + C_m B_1 + \dots + C_{m-n+1} B_n)x^{m+1} = 0 \quad (i+j = m+1)
\dots
\displaystyle (C_{m+n} B_0 + C_{m+n-1} B_1 + \dots + C_m B_n)x^{m+n} = 0 \quad (i+j = m+n)

Using A_k = 0 for k > m. We can split this system of equations in two parts and get rid of the x-terms. One part containing the A_i coefficients (up to i+j = m) and the part without A_i coefficients. The first system of equation is:

\displaystyle C_0 B_0 - A_0 = 0
\displaystyle C_1 B_0 + C_0 B_1 - A_1 = 0
\displaystyle C_2 B_0 + C_1 B_1 + C_0 B_2 - A_2 = 0
\dots
\displaystyle C_m B_0 + C_{m-1} B_1 + \dots + C_{m-n} B_n - A_m = 0

The second system of equations:

\displaystyle C_{m+1} B_0 + C_m B_1 + \dots + C_{m-n+1} B_n = 0
\dots
\displaystyle C_{m+n} B_0 + C_{m+n-1} B_1 + \dots + C_m B_n = 0

Setting B_0 = 1 without loss of generality the second system of equations becomes:

\displaystyle C_m B_1 + \dots + C_{m-n+1} B_n = -C_{m+1}
\displaystyle C_{m+1} B_1 + \dots + C_{m-n+2} B_n = -C_{m+2}
\dots
\displaystyle C_{m+n-1} B_1 + \dots + C_m B_n = -C_{m+n}

Changing the order of the coefficients gives:

\displaystyle C_{m-n+1} B_n + \dots + C_m B_1 = -C_{m+1}
\displaystyle C_{m-n+2} B_n + \dots + C_{m+1} B_1 = -C_{m+2}
\dots
\displaystyle C_m B_n + \dots + C_{m+n-1} B_1 = -C_{m+n}

We can write both systems using matrix notation:

\displaystyle \left( \begin{array}{cccccc} C_0 & 0 & \dots & & & 0 \\ C_1 & C_0 & 0 & \dots & & 0 \\ C_2 & C_1 & C_0 & 0 & \dots & 0 \\ \vdots & & & & & \vdots \\ C_m & \dots & & & \dots & C_{m-n} \end{array} \right) \left( \begin{array}{c} B_0 \\ B_1 \\ B_2 \\ \vdots \\ B_n \end{array} \right) = \left( \begin{array}{c} A_0 \\ A_1 \\ A_2 \\ \vdots \\ A_m \end{array} \right)
\displaystyle \left( \begin{array}{ccccc} C_{m-n+1} & C_{m-n+2} & \dots & & C_m \\ C_{m-n+2} & C_{m-n+3} & \dots & & C_{m+1} \\ \vdots & & & & \vdots \\ C_m & C_{m+1} & & \dots & C_{m+n-1} \end{array} \right) \left( \begin{array}{c} B_n \\ B_{n-1} \\ \vdots \\ B_1 \end{array} \right) = - \left( \begin{array}{c} C_{m+1} \\ C_{m+2} \\ \vdots \\ C_{m+n} \end{array} \right)

Practically, we start calculations where m-n+k \geq 0. The resolution of the second system gives the values of the B_j coefficients which are injected into the first system to obtain the A_i coefficient (using B_0 = 1). So if the following determinant (called a ‘Hankel determinant’):

\displaystyle H_{n,m}(f) := \left| \begin{array}{ccccc} C_{m-n+1} & C_{m-n+2} & \dots & & C_m \\ C_{m-n+2} & C_{m-n+3} & \dots & & C_{m+1} \\ \vdots & & & & \vdots \\ C_m & C_{m+1} & & \dots & C_{m+n-1} \end{array} \right|

Is not null. i.e.:

\displaystyle H_{n,m}(f) \neq 0

the Padé approximant will be unique. The calculations for the Padé approximant coefficients yield a unique solution for any Taylor or Maclaurin series, provided the Hankel determinant is not null, excluding cases of degeneracy due to specific coefficient values.

Comments

Leave a comment

More posts