How would I generate a Pade Approximant using the coefficients of a Taylor Series?

Question

I would like to find an effective way to make a Pade Approximant using the coefficients of a Taylor Series. I've heard of Wynn's epsilon algorithm and using the Extended Euclidean Algorithm, but what is the most efficient way and how is it done step by step? Any help or guidance would be appreciated.

Answer 1

Wynn's epsilon method is an iterative algorithm which accelerates a slowly converging series. It can be used to produce Padé approximants from a sequence of Taylor series polynomials.

The core of the algorithm is quite simple. Let $D_0$ be a sequence of length $n$ that we wish to accelerate. $n$ must be odd. Let $D_{-1}$ be a sequence of zeros (of length $n-1$). The outer loop of the algorithm produces sequence $D_{i+1}$ from sequence $D_{i-1}$ and sequence $D_i$.

The inner loop of the algorithm steps through the items of $D_{i-1}$ and $D_i$. It gets the difference of adjacent pairs of items of $D_i$, adding the reciprocal of that difference to the corresponding term of $D_{i-1}$ to produce items of sequence $D_{i+1}$. The length of $D_{i+1}$ is one less than the length of $D_i$, so at the end of the algorithm, $D_{n-1}$ consists of a single item. When $i$ is even, $D_i$ is an accelerated version of the sequence.

I'll illustrate the algorithm in Python, using the classic Leibniz series for $\pi$, which is essentially the Taylor series for $4\arctan(1)$: $$\pi = 4 - \frac43 + \frac45 - \frac47 + \cdots$$ This series converges very slowly: after two million terms we barely have six decimals.

The inner loop of the algorithm consists of this list comprehension:

[p + 1 / (r - q) for p, q, r in zip(pd, d, d1)]   

pd is the previous $D$, d and d1 give us the adjacent pairs of items from the current $D$.

Code

""" Wynn's epsilon algorithm
    Leibniz pi. Plain Python
    Written by PM 2Ring 2022.06.22
"""

# Number of terms. Must be odd.
# ZeroDivisionError when num > 21
num = 9

# Create base sequence
a = p = 4.0
d = [a]
for i in range(3, 2 * num, 2):
    p = -p
    a += p / i
    d.append(a)

print("Original", *d, "", sep="\n")

# Accelerate using Wynn's epsilon algorithm
pd = [0] * num
for i in range(1, num):
    d1 = d[1:]
    pd, d = d1, [p + 1 / (r - q) for p, q, r in zip(pd, d, d1)]
    print(i, *d, "", sep="\n")

Output

Original
4.0
2.666666666666667
3.466666666666667
2.8952380952380956
3.3396825396825403
2.9760461760461765
3.2837384837384844
3.017071817071818
3.2523659347188767

1
-0.7500000000000002
1.2500000000000002
-1.7500000000000009
2.249999999999999
-2.749999999999999
3.2499999999999973
-3.750000000000001
4.249999999999999

2
3.166666666666667
3.1333333333333337
3.1452380952380956
3.13968253968254
3.1427128427128435
3.1408813408813416
3.142071817071818

3
-28.750000000000107
82.2500000000003
-177.75000000000065
327.24999999996163
-542.7500000000019
836.250000000003

4
3.1423423423423427
3.1413919413919418
3.141662737702342
3.141563417425487
3.141606504043053

5
-969.9375000000699
3515.0624999998927
-9741.187500039121
22666.312500039647

6
3.1416149068322987
3.141587301587302
3.1415942744801812

7
-32709.93749997997
133671.31249657096

8
3.1415933118799284

The final value, $3.1415933118799284$, is correct to six decimals. We can get more decimals by increasing num, but if we make it too large the difference r - q becomes zero, leading to a ZeroDivisionError. So if you want more precision you need to use the Python decimal library, or some other library supporting high precision arithmetic, eg mpmath.

---

But how can we use this algorithm to produce Padé approximants? The procedure is identical, we just need to transform a list of Taylor series sums. To implement this from scratch in plain Python, you will need to write code that can do arithmetic with polynomials of a single variable. Then you will need to implement rational functions. Your code will need to implement the Euclidean algorithm, so you can calculate the GCD of two polynomials so that you can reduce the rational functions to their lowest terms. You can store polynomial coefficients in a list, but I strongly recommend creating classes for your polynomials and rational functions.

I'm not going to implement such classes here. It's a fair bit of work, and it's really not on-topic for math.SE. But I can demonstrate the procedure using SageMath, a free open-source mathematics software system built on top of Python. SageMath makes it easy to manipulate polynomials and other symbolic expressions. For convenience, Sage lets us use ^ for exponentiation.

The following code produces Padé approximants from the Taylor series for $\arctan(x)$.

Code

""" Wynn's epsilon algorithm
    Arctan Padé
    Written by PM 2Ring 2023.12.10
"""

var('x')

# Number of terms. Must be odd.
num = 7

# Create base sequence
a = p = x
d = [a]
for i in range(3, 2 * num, 2):
    p *= -x^2
    a += p / i
    d.append(a)

print("Original", *d, "", sep="\n")

# Accelerate using Wynn's epsilon algorithm
pd = [0] * num
for i in range(1, num):
    d1 = d[1:]
    pd, d = d1, [(p + 1 / (r - q)).simplify_full() for p, q, r in zip(pd, d, d1)]
    print(i, *d, "", sep="\n")

pade(x) = d[0]
print("Padé", pade)

g(x) = pade(x) - arctan(x)
u = sqrt(2.0) - 1
print("pi ~=", 8 * pade(u))

dom = -u, u
P = plot(g, dom, frame=True, gridlines=True)
P.show()

Output

Original
x
-1/3*x^3 + x
1/5*x^5 - 1/3*x^3 + x
-1/7*x^7 + 1/5*x^5 - 1/3*x^3 + x
1/9*x^9 - 1/7*x^7 + 1/5*x^5 - 1/3*x^3 + x
-1/11*x^11 + 1/9*x^9 - 1/7*x^7 + 1/5*x^5 - 1/3*x^3 + x
1/13*x^13 - 1/11*x^11 + 1/9*x^9 - 1/7*x^7 + 1/5*x^5 - 1/3*x^3 + x

1
-3/x^3
5/x^5
-7/x^7
9/x^9
-11/x^11
13/x^13

2
1/3*(4*x^3 + 15*x)/(3*x^2 + 5)
-1/15*(4*x^5 - 40*x^3 - 105*x)/(5*x^2 + 7)
1/105*(12*x^7 - 56*x^5 + 420*x^3 + 945*x)/(7*x^2 + 9)
-1/315*(20*x^9 - 72*x^7 + 252*x^5 - 1680*x^3 - 3465*x)/(9*x^2 + 11)
1/3465*(140*x^11 - 440*x^9 + 1188*x^7 - 3696*x^5 + 23100*x^3 + 45045*x)/(11*x^2 + 13)

3
-5/4*(15*x^4 + 42*x^2 + 35)/x^7
7/4*(35*x^4 + 90*x^2 + 63)/x^9
-9/4*(63*x^4 + 154*x^2 + 99)/x^11
11/4*(99*x^4 + 234*x^2 + 143)/x^13

4
1/15*(64*x^5 + 735*x^3 + 945*x)/(15*x^4 + 70*x^2 + 63)
-1/105*(64*x^7 - 1344*x^5 - 9765*x^3 - 10395*x)/(35*x^4 + 126*x^2 + 99)
1/315*(64*x^9 - 576*x^7 + 8064*x^5 + 47355*x^3 + 45045*x)/(63*x^4 + 198*x^2 + 143)

5
-21/64*(175*x^8 + 1260*x^6 + 3690*x^4 + 4620*x^2 + 2079)/x^11
27/64*(735*x^8 + 4620*x^6 + 11242*x^4 + 12012*x^2 + 4719)/x^13

6
1/35*(256*x^7 + 5943*x^5 + 19250*x^3 + 15015*x)/(35*x^6 + 315*x^4 + 693*x^2 + 429)

Padé x |--> 1/35*(256*x^7 + 5943*x^5 + 19250*x^3 + 15015*x)/(35*x^6 + 315*x^4 + 693*x^2 + 429)
pi ~= 3.14159265432570

Here's a plot showing the error of our Padé:

You can run these scripts on the SageMathCell server. Leibniz pi. Arctan Padé.