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.
How to generate a Pade Approximant using the coefficients of a Taylor Series
computational mathematicspade approximation
Best Answer
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:
pd
is the previous $D$,d
andd1
give us the adjacent pairs of items from the current $D$.Code
Output
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 differencer - q
becomes zero, leading to a ZeroDivisionError. So if you want more precision you need to use the Pythondecimal
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
Output
Here's a plot showing the error of our Padé:
You can run these scripts on the SageMathCell server. Leibniz pi. Arctan Padé.