^{Department of Physics and Astronomy}_{The Forbes Group
Pade Approximants}

Padé Approximants¶

A Padé approximant is a rational functions $A(x) = P(x)/Q(x)$ where $P(x)$ and $Q(x)$ are polinomials. The provide an alternative to Taylor series, often with better behaviour. Here we consider using them to model functions with well-defined behaviour both at small and large parameter values.

Our running example will be the equation of state for a Fermi gas with S-wave scattering length $a$. We shall consider the region from weakly interacting – the BCS region – to the strongly interacting unitary Fermi gas (UFG). The equation of state will be expressed in terms of a dimensionless parameter $x = k_F a$ where $k_F$ is the Fermi momentum and $a$ is the (negative) scattering length. The equation of state will be normalized by the energy density of a non-interacting gas. We then have the following limits:

\begin{align} \frac{\mathcal{E}}{\mathcal{E}_{FG}} &= 1 + \frac{10}{9\pi}x + \frac{4(11-2\ln 2)}{21\pi^2} x^2 + \order(x^3), \tag{when $x$ is small}\\ &= \xi - \frac{\zeta}{x} + \order(x^{-2}).\tag{when $x$ is large} \end{align}

Our goal is to find a Padé approximant reproducing both of these behaviours. Let's simplify this to:

\begin{align} \frac{\mathcal{E}}{\mathcal{E}_{FG}} &= \sum_{n=0}^{N_a} a_n x^n,\tag{when $x$ is small}\\ &= \sum_{n=0}^{N_b} b_n x^{-n},\tag{when $x$ is large} \end{align}

In [4]:

from IPython.display import display
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
a = [1.0, 10./9./np.pi, 4*(11 - 2*np.log(2))/12/np.pi**2]
b = [0.370, -0.9, 0.05]
x = np.linspace(-5, -0.01, 100)
plt.plot(x, np.polyval(a[::-1], x), label='BCS')
plt.plot(x, np.polyval(b[::-1], 1./x), label='UFG')
plt.ylim(0.37, 1)
plt.xlabel(r"$k_F a$")
plt.ylabel(r"$E/E_{FG}$")
plt.legend(loc='best')

Out[4]:

<matplotlib.legend.Legend at 0x128f5cf50>

Our goal is to find a nice smooth interpolating function that works over this range of $k_F a$. We have 6 parameters to fit here, and constants at both asymptotes which means that $P(x)$ and $Q(x)$ must have the same order $N$. Removing the overall normalization which will cancell, we have $2N - 1$ degrees of freedom, so we need $N=3$ and will have one free parameter to play with.

In [5]:

%install_ext http://bitbucket.org/mforbes/notes_ipython/raw/tip/notebook_utils/mmf_ext.py
%load_ext mmf_ext
import sympy
from sympy import S, var
sympy.init_printing(use_latex=False)

var('x')
var('a_0, a_1, a_2, a_3')
var('b_0, b_1, b_2')
var('p_0, p_1, p_2, p_3')
var('q_1, q_2, q_3')

BCS = a_0 + a_1*x + a_2*x**2 + a_3*x**3
UFG = b_0 + b_1/x + b_2/x**2
P = p_0 + p_1*x + p_2*x**2 +  p_3* x**3
Q = 1 + q_1*x + q_2*x**2 +  q_3* x**3
A = P/Q
%math A(x) = A

Installed mmf_ext.py. To use it, type:
  %load_ext mmf_ext
The mmf_ext extension is already loaded. To reload it, use:
  %reload_ext mmf_ext

$$A(x) = \frac{p_{0} + p_{1} x + p_{2} x^{2} + p_{3} x^{3}}{q_{1} x + q_{2} x^{2} + q_{3} x^{3} + 1}$$

In [6]:

x, y = sympy.symbols('x, y')
N = 4
M = N
_as = sympy.symbols(map("a_{0}".format, range(M)))
_bs = sympy.symbols(map("b_{0}".format, range(2*N-M-1)))
_ps = sympy.symbols(map("p_{0}".format, range(N)))
_qs = sympy.symbols(map("q_{0}".format, range(N)))
_qs[0] = 1
A = np.polyval(_as[::-1], x)
B = np.polyval(_bs[::-1], 1/x)
P = np.polyval(_ps[::-1], x)
Q = np.polyval(_qs[::-1], x)
P/Q

Out[6]:

p₀ + x⋅(p₁ + x⋅(p₂ + p₃⋅x))
───────────────────────────
 x⋅(q₁ + x⋅(q₂ + q₃⋅x)) + 1

Here is a direct solution using series and letting sympy do some hard work (the equations for $p_i$ and $q_i$ are non-linear in this approach.)

In [7]:

eqsA = sympy.poly((P/Q-A).series(x, n=len(_as)).removeO(), x).coeffs()
sol = sympy.solve(eqsA, _ps, dict=True)[0]; display(sol)
(P/Q).subs(sol).series(x, n=len(_as)).removeO().simplify()

{p₀: a₀, p₁: a₀⋅q₁ + a₁, p₂: a₀⋅q₂ + a₁⋅q₁ + a₂, p₃: a₀⋅q₃ + a₁⋅q₂ + a₂⋅q₁ + a
₃}

Out[7]:

                2       3
a₀ + a₁⋅x + a₂⋅x  + a₃⋅x

Now we do it explicitly by collecting the coefficients of the polynomial to obtain a set of linear equations. Note that the results are the same.

In [8]:

eqsA = sympy.poly((P - A*Q), x).coeffs()[::-1][:len(_as)]
sol = sympy.solve(eqsA, _ps, dict=True)[0]; display(sol)
(P/Q).subs(sol).series(x, n=len(_as)).removeO().simplify()

{p₀: a₀, p₁: a₀⋅q₁ + a₁, p₂: a₀⋅q₂ + a₁⋅q₁ + a₂, p₃: a₀⋅q₃ + a₁⋅q₂ + a₂⋅q₁ + a
₃}

Out[8]:

                2       3
a₀ + a₁⋅x + a₂⋅x  + a₃⋅x

Now for the denominator:

In [9]:

eqsB = sympy.poly((P/Q-B).subs(x, 1/y).series(y, n=len(_bs)).removeO(), y).coeffs()
sol = sympy.solve(eqsB, _qs[1:], dict=True)[0]; display(sol)
(P/Q).subs(sol).subs(x, 1/y).series(y, n=len(_bs)).removeO().simplify()

⎧                    2                   2                           ⎫
⎪    b₀⋅p₁⋅p₃ - b₀⋅p₂  + b₁⋅p₂⋅p₃ - b₂⋅p₃       b₀⋅p₂ - b₁⋅p₃      p₃⎪
⎨q₁: ─────────────────────────────────────, q₂: ─────────────, q₃: ──⎬
⎪                      2                               2           b₀⎪
⎩                    b₀ ⋅p₃                          b₀              ⎭

Out[9]:

                2
b₀ + b₁⋅y + b₂⋅y

And the smart way:

In [51]:

eqsB = sympy.poly(((P-B*Q)*x**(len(_bs)-1)).expand(), x).coeffs()[:len(_bs)]
sol = sympy.solve(eqsB, _qs[1:], dict=True)[0]; display(sol)
(P/Q).subs(sol).subs(x, 1/y).series(y, n=len(_bs)).simplify()

⎧      2                    ⎛          2⎞                           ⎫
⎪    b₀ ⋅p₁ - b₀⋅b₁⋅p₂ - p₃⋅⎝b₀⋅b₂ - b₁ ⎠      b₀⋅p₂ - b₁⋅p₃      p₃⎪
⎨q₁: ────────────────────────────────────, q₂: ─────────────, q₃: ──⎬
⎪                      3                              2           b₀⎪
⎩                    b₀                             b₀              ⎭

Out[51]:

    2  2                           2             2  2        
  b₁ ⋅y        2          2⋅b₁⋅p₂⋅y         b₀⋅p₂ ⋅y     ⎛ 3⎞
- ────── + b₂⋅y  + b₁⋅y + ────────── + b₀ - ───────── + O⎝y ⎠
    b₀                        p₃                 2           
                                               p₃

This SEEMS TO FAIL due to Issue #9173. Be sure to simplify first before calling `series or limit!

In [52]:

(P/Q).subs(sol).subs(x, 1/y).simplify().series(y, n=len(_bs)).simplify()

Out[52]:

    2                ⎛ 3⎞
b₂⋅y  + b₁⋅y + b₀ + O⎝y ⎠

Here is the complete solution:

In [36]:

eqsA = sympy.poly((P - A*Q), x).coeffs()[::-1][:len(_as)]
eqsB = sympy.poly(((P-B*Q).subs(x, 1/y)*y**(N-1)).expand(), y).coeffs()[::-1][:len(_bs)]
sols = sympy.solve(eqsA + eqsB, _ps + _qs[1:])

In [49]:

as_ = np.random.rand(len(_as))
bs_ = np.random.rand(len(_bs))
vals = dict(zip(_as, as_)+zip(_bs, bs_))
vals.update({_s: sols[_s].subs(vals) for _s in sols})
display((P/Q - A).subs(vals).simplify().series(x, n=len(_as)).simplify(),
 (P/Q-B).subs(vals).subs(x, 1/y).simplify().series(y, n=2+0*len(_bs)).simplify())

                                               2                         3    
1.11022302462516e-16⋅x + 2.10229541421588e-16⋅x  + 3.45642189044409e-16⋅x  + O

⎛ 4⎞
⎝x ⎠

                                                 ⎛ 2⎞
8.78967335516104e-17 + 5.06570180747383e-16⋅y + O⎝y ⎠

The equations for the coefficients are easily found by multiplying the approrpiate polynomials and collecting terms. For the BCS limit this is simply $P(x) = Q(x)f(x)$ and then we collect the lowest powers in $x$. For the UFG limit, we do the same, but multiply by the highest inverse power and then collect the highest powers. Thus:

In [13]:

eqs = (sympy.Poly(P - Q*BCS, x).coeffs()[::-1][:3] 
       + sympy.Poly((P - UFG*Q)*x**2, x).coeffs()[:3])
sols = sympy.solve(eqs, [p_0, p_1, p_2, q_1, q_2, q_3])
A_p3 = A.subs(sols).simplify()

In [25]:

display(A.series(x, n=3), 
        A.subs(x, 1/x).series(x, n=3).subs(x, 1/x))
sols_q = sympy.solve(sympy.Poly((P - UFG*Q)*x**2, x).coeffs()[:3], [q_1, q_2, q_3])
A.subs(sols_q).subs(x, 1/x).series(x, n=3).subs(x, 1/x).simplify().zeries

$$x \left(- p_{0} q_{1} + p_{1}\right) + x^{2} \left(p_{0} q_{1}^{2} - p_{0} q_{2} - p_{1} q_{1} + p_{2}\right) + p_{0} + \mathcal{O}\left(x^{3}\right)$$

$$\frac{1}{x^{2}} \left(\frac{p_{1}}{q_{3}} - \frac{p_{2} q_{2}}{q_{3}^{2}} - \frac{p_{3} q_{1}}{q_{3}^{2}}\right) + \frac{1}{x} \left(\frac{p_{2}}{q_{3}} - \frac{p_{3} q_{2}}{q_{3}^{2}}\right) + \frac{p_{3}}{q_{3}} + \mathcal{O}\left(\frac{1}{x^{3}}; x\rightarrow\infty\right)$$

Out[25]:

$$- \frac{b_{1}^{2}}{b_{0} x^{2}} + \frac{b_{2}}{x^{2}} + \frac{b_{1}}{x} + \frac{2 b_{1} p_{2}}{p_{3} x^{2}} + b_{0} - \frac{b_{0} p_{2}^{2}}{p_{3}^{2} x^{2}} + \mathcal{O}\left(\frac{1}{x^{3}}; x\rightarrow\infty\right)$$

Now let's try some values. The parameter $p_3$ is still undetermined, but we must be careful to ensure that there are no poles in the region of interest $x<0$. This is ensured by choosing $p_3\lessapprox -0.04$:

In [166]:

a = [1.0, 10./9./np.pi, 4*(11 - 2*np.log(2))/12/np.pi**2]
b = [0.370, -0.9, 0.05]
subs = dict(zip([a_0, a_1, a_2], a) + zip([b_0, b_1, b_2], b))
subs[p_3] = -1
np.roots(sympy.Poly(Q.subs(sols).subs(subs)).coeffs())

Out[166]:

array([ 5.00709122,  2.02189374,  0.03654752])

In [169]:

interp = sympy.lambdify(
    [x, p_3], A_p3.subs(subs), np)

_x = np.linspace(-10, -0.001, 1000)
plt.plot(_x, np.polyval(a[::-1], _x), label='BCS')
plt.plot(_x, np.polyval(b[::-1], 1./_x), label='UFG')
for p3 in np.linspace(-1000, -0.1,20):
    plt.plot(_x, interp(_x, p3))
plt.ylim(0.37, 1)
plt.xlabel(r"$k_F a$")
plt.ylabel(r"$E/E_{FG}$")
plt.legend(loc='best')

Out[169]:

<matplotlib.legend.Legend at 0x1215d8a90>

Note that the interpolation is largely insensitive to the value of $p_3$ since the poles are all at positive values of $x$.

In [ ]: