Approximaths

Author: approximaths

  • ‘Generic’ summation I

    Let imagine that we have engineered a summation ‘machine’ called \mathcal{S}() and consider the geometric series:

    \displaystyle 1+x+x^2+x^3+x^4+\dots

    Apply our summation ‘machine’ on this series and call the result s:

    \displaystyle s = \mathcal{S}(1+x+x^2+x^3+x^4+\dots)

    We assign the following two properties for the ‘machine’ \mathcal{S}:

    \displaystyle \mathcal{S}(a_0+a_1+a_2+\dots) = a_0 + \mathcal{S}(a_1+a_2+\dots) \quad \text{(first property)}
    \displaystyle \mathcal{S}(\sum(\alpha a_n) + \sum(\beta b_n)) = \alpha \mathcal{S}(\sum a_n) + \beta \mathcal{S}(\sum b_n)) \quad \text{(second property)}

    where \alpha and \beta are constants. Equipped with the machine \mathcal{S}() and its two properties consider again the geometric series:

    \displaystyle \begin{aligned} s &= \mathcal{S}(1+x+x^2+x^3+x^4+\dots) &\quad\text{(definition)}\\ s &= 1+ \mathcal{S}(x+x^2+x^3+x^4+\dots) &\quad\text{(first property)}\\ s &= 1+ x\mathcal{S}(1+x+x^2+x^3+\dots) &\quad\text{(second property)}\\ s &= 1+ xs &\quad\text{(definition)}\\ s &= \frac{1}{1-x} \end{aligned}

    Consider the series:

    \displaystyle 1-1+1-1+1-1+\dots

    If we use traditional, ‘rigorous’, summation techniques we will conclude that this alternating series is not converging. Now apply the summation “machine” as described above:

    \displaystyle \begin{aligned} s &= \mathcal{S}(1-1+1-1+1-1+\dots ) \\ s &= 1+ \mathcal{S}(-1+1-1+1-1+\dots ) \\ s &= 1 -\mathcal{S}(1-1+1-1+1-\dots ) \\ s &= 1-s \\ s &= \frac{1}{2} \end{aligned}

  • Borel summation

    Let introduce the Borel summation by first recalling the following property:

    \displaystyle n!= \int_{0}^{\infty} e^{-t}t^{n}dt
    \displaystyle 1 = \frac{\int_{0}^{\infty} e^{-t}t^{n}dt}{n!}

    We would like to sum the following series:

    \displaystyle \sum_{n=0}^{\infty}a_n = a_0 + a_1 + a_2 + a_3 + a_4 +...
    \displaystyle \phantom{\sum_{n=0}^{\infty}a_n} = 1a_0 + 1a_1 + 1a_2 + 1a_3 + 1a_4 +...
    \displaystyle \phantom{\sum_{n=0}^{\infty}a_n} = \sum_{n=0}^{\infty} \frac{\int_{0}^{\infty} e^{-t}t^{n}dt}{n!} a_n

    The Borel sum B is defined as follow:

    \displaystyle B := \sum_{n=0}^{\infty} \frac{\int_{0}^{\infty} e^{-t}t^{n}dt}{n!} a_n
    \displaystyle \phantom{B :=} = \int_{0}^{\infty}e^{-t} \sum_{n=0}^{\infty} \frac{t^n}{n!}a_n\,dt

    The factorial term makes this sum with better chance to converge. In the section concerning the Euler summation we have seen that:

    \displaystyle E(1 - 1 + 1 - 1 + 1 - 1 + ...) = \frac{1}{2}

    Let calculate the corresponding Borel sum:

    \displaystyle B(1 - 1 + 1 - 1 + 1 - 1 + ...) = \int_{0}^{\infty}e^{-t} \sum_{n=0}^{\infty} \frac{t^n}{n!}(-1)^n \,dt
    \displaystyle \phantom{B(1 - 1 + 1 - 1 + 1 - 1 + ...)} = \int_{0}^{\infty}e^{-t} \sum_{n=0}^{\infty} \frac{(-t)^n}{n!} \,dt
    \displaystyle \phantom{B(1 - 1 + 1 - 1 + 1 - 1 + ...)} = \int_{0}^{\infty} e^{-t}e^{-t}\,dt
    \displaystyle \phantom{B(1 - 1 + 1 - 1 + 1 - 1 + ...)} = \int_{0}^{\infty} e^{-2t}\,dt
    \displaystyle \phantom{B(1 - 1 + 1 - 1 + 1 - 1 + ...)} = -\frac{1}{2} e^{-2t}\Big|_{0}^{\infty}
    \displaystyle \phantom{B(1 - 1 + 1 - 1 + 1 - 1 + ...)} = 0 - (-\frac{1}{2})
    \displaystyle \phantom{B(1 - 1 + 1 - 1 + 1 - 1 + ...)} = \frac{1}{2}

    and therefore:

    \displaystyle B(1-1+1-1+1-1+...) = E(1-1+1-1+1-1+...)

    Consider the series:

    \displaystyle -1+2!-3!+4!-5!+6!-...
    \displaystyle B(-1+2!-3!+4!-...) = \int_{0}^{\infty}e^{-t} \sum_{n=1}^{\infty} \frac{t^n}{n!}(-1)^n n! \,dt
    \displaystyle \phantom{B(-1+2!-3!+4!-...)} = \int_{0}^{\infty}e^{-t} \sum_{n=1}^{\infty} (-t)^n \,dt

    \sum_{n=0}^{\infty} (-t)^n converges for t < 1 to -\frac{1}{1+t} therefore:

    \displaystyle B = -\int_{0}^{\infty} \frac{te^{-t}}{(1+t)} \,dt \approx -0.40365

    If a series is summable in the sense of Euler, then it is also summable in the sense of Borel, and both summation methods yield the same value. The converse is false: there exist series that are summable in the sense of Borel but not in the sense of Euler. In other words, Borel summation is more powerful, as it applies to more strongly divergent series.

  • Euler summation II

    In the previous post, we introduced Euler summation. The following are two examples where it fails to produce a finite result.

    Consider the divergent series:

    \displaystyle 0 + 1 + 2 + 3 + 4 + \dots

    Define:

    \displaystyle \begin{array}{rcl} f(x) &=& 0x^{0} + 1x^{1} + 2x^{2} + 3x^{3} + 4x^{4} + \dots \\ &=& \sum_{n=0}^{\infty} nx^{n} = \frac{x}{(1-x)^{2}} \\ \end{array}

    The Euler sum is:

    \displaystyle E = \lim_{x \to 1_{-}} \frac{x}{(1-x)^{2}} = \infty
    \displaystyle E(0+ 1 + 2 + 3 + 4 + \dots) = \infty

    Now consider the divergent series:

    \displaystyle 1 + 4 + 9 + 16 + 25 + 36 + \dots

    Define:

    \displaystyle \begin{array}{rcl} f(x) &=& 1^{2}x^{1}+2^{2}x^{2}+3^{2}x^{3}+4^{2}x^{4}+5^{2}x^{5}+6^{2}x^{6}+\dots \\ &=& \sum_{n=1}^{\infty} n^{2} x^{n} = \frac{x(1+x)}{(1-x)^{3}} \\ \end{array}

    The Euler sum is:

    \displaystyle E = \lim_{x \to 1_{-}}\frac{x(1+x)}{(1-x)^{3}} = \infty
    \displaystyle E(1^{2} + 2^{2} + 3^{2} + 4^{2} + \dots) = \infty

    Advantages

    Regularization of slowly divergent series:
    Euler summation can assign a finite value to some divergent series that oscillate or diverge slowly, such as

    \displaystyle 1 - 1 + 1 - 1 + \dots

    where E(series) = \tfrac{1}{2}.

    Improved convergence:
    For many convergent series, Euler transformation accelerates convergence, making it useful for numerical computations.

    Analytic continuation link:
    It provides a bridge between ordinary summation and more advanced summation methods (e.g. Borel or zeta regularization).

    Disadvantages

    Limited domain of applicability:
    Euler summation fails for series that diverge too rapidly, such as

    \displaystyle 1 + 2 + 3 + 4 + \dots

    where E(series) = \infty.

    Not uniquely defined for all divergent series:
    Some series cannot be assigned a finite Euler sum, or the method may yield inconsistent results depending on the transformation order.

    Weaker than analytic regularization:
    Compared to zeta or Borel summation, Euler’s method handles fewer classes of divergent series and lacks a rigorous analytic continuation framework.

  • Euler summation I

    So far we’ve seen Padé’s approximants. These have enabled us to approximate a function from its corresponding Taylor series and then transform this (potentially non-convergent) series into a convergent rational fraction.

    We’d like to introduce other techniques that could be used to sum non-convergent series. First, we’ll take a look at Euler summation.

    If the series \sum_{n=0}^{\infty} a_n is algebraically divergent (the terms blow up like some power of n), then the series:

    \displaystyle f(x) = \sum_{n=0}^{\infty}a_n x^n

    converges for x \in (-1,1). If the limit

    \displaystyle E := \lim_{x \to 1_{-}}f(x)

    exists and is finite then it is called the Euler sum E of the original series.

    For example, consider the divergent series:

    \displaystyle 1 - 1 + 1 - 1 + 1 - 1 + ...

    and multiply each n-term with x^n (starting at n = 0):

    \displaystyle f(x) = x^0 - x^1 + x^2 - x^3 + x^4 - x^5 + ...
    \displaystyle = \frac{1}{1-(-x)} = \frac{1}{1+x}
    \displaystyle E = \lim_{x \to 1_{-}}{\frac{1}{1+x}} = \frac{1}{2}

    Therefore (according to Euler summation):

    \displaystyle 1 - 1 + 1 - 1 + 1 - 1 + ... = \frac{1}{2}

  • Padé approximants: Convergence III

    For a function f(z) analytic at z=0, with Taylor series

    \displaystyle f(z) = \sum_{k=0}^{\infty} a_k z^k

    valid within its radius of convergence R, the Padé approximant of order [m/n], denoted \frac{P_m(z)}{Q_n(z)} with \deg P_m \leq m and \deg Q_n \leq n, satisfies

    \displaystyle f(z) - \frac{P_m(z)}{Q_n(z)} = O(z^{m+n+1})

    near z=0. The rational structure of \frac{P_m(z)}{Q_n(z)} allows it to approximate f(z) beyond the disk |z| < R by modeling singularities (e.g., poles or branch points) through the zeros of Q_n(z). This enables analytic continuation into regions where the Taylor series diverges.

    Formally, for a meromorphic function f(z) in a domain D, the diagonal Padé approximants [n/n] often converge to f(z) in D \setminus S, where S is the set of poles of f:

    Let f(z) be meromorphic in a domain D \subseteq \mathbb{C}, with a set of poles S of finite total multiplicity. The diagonal Padé approximants [n/n], defined as rational functions \frac{P_n(z)}{Q_n(z)} satisfying

    \displaystyle f(z) - \frac{P_n(z)}{Q_n(z)} = O(z^{2n+1})

    near z=0, converge uniformly to f(z) on compact subsets of D \setminus S as n \to \infty.

    The zeros of Q_n(z) approximate the poles in S, enabling analytic continuation of f(z) beyond the radius of convergence of its Taylor series.

  • Nuttall’s Padé approximant

    Let f(z) = \sum_{k=0}^\infty C_k z^k be a power series. The denominator Q_{n-1}(z) of the Padé approximant P(n, n-1) is given by Nuttall’s compact form:

    \displaystyle Q_{n-1}(z) = \frac{ \begin{vmatrix} C_0 & C_1 & C_2 & \cdots & C_{n-1} & C_n \\ C_1 & C_2 & C_3 & \cdots & C_n & C_{n+1} \\ C_2 & C_3 & C_4 & \cdots & C_{n+1} & C_{n+2} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ C_{n-2} & C_{n-1} & C_n & \cdots & C_{2n-3} & C_{2n-2} \\ 1 & z & z^2 & \cdots & z^{n-2} & z^{n-1} \end{vmatrix} }{ \begin{vmatrix} C_0 & C_1 & C_2 & \cdots & C_{n-1} \\ C_1 & C_2 & C_3 & \cdots & C_n \\ C_2 & C_3 & C_4 & \cdots & C_{n+1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ C_{n-1} & C_n & C_{n+1} & \cdots & C_{2n-1} \end{vmatrix} }

    The numerator P_n(z) is obtained by satisfying the Padé approximation condition: f(z) Q_{n-1}(z) - P_n(z) = O(z^{2n}).

    The compact form of Nuttall’s Padé approximant P(n, n-1) is particularly valuable in numerical analysis and theoretical physics for its efficiency in computing Padé approximants without explicitly solving large linear systems.

    By expressing the denominator Q_{n-1}(z) as a ratio of determinants, it provides a direct and elegant method to capture the approximant’s poles, which is crucial for analyzing singularities of functions, especially in Stieltjes series or meromorphic functions.

    This formulation simplifies calculations, facilitates the study of convergence properties, and connects Padé approximants to orthogonal polynomials, enabling applications in areas like quantum field theory and asymptotic analysis where rapid computation and singularity detection are essential.

  • Continued fractions II

    A well-known continued fraction representation of exp(x) is:

    \displaystyle exp(x) = 1 + \frac{x}{1 - \frac{x}{2 + \frac{x}{3 - \frac{x}{4 + \cdots}}}}

    Using the procedure described in the previous post, we can show that the near-diagonal Padé coefficients presented in this post can be converted to a sequence of truncated fractions corresponding to the continued fraction above (see figures below).

    …

    P(m,n) 0 1 2 3
    0
    \displaystyle 1
    \displaystyle 1 + x
    		 …
    1
    \displaystyle 1 + \cfrac{x}{1 - \cfrac{x}{2}}
    \displaystyle 1 + \cfrac{x}{1 - \cfrac{x}{2 + \cfrac{x}{3}}}
    		 …
    2
    \displaystyle 1 + \cfrac{x}{1 - \cfrac{x}{2 + \cfrac{x}{3 - \cfrac{x}{4}}}}
    \displaystyle 1 + \cfrac{x}{1 - \cfrac{x}{2 + \cfrac{x}{3 - \cfrac{x}{4 + \cfrac{x}{5}}}}}
    		 …
    3 …
    … … … … … …

    Caption: Padé Approximants P(m,n) for exp(x) expressed as continued fractions.

    The relationship between Padé approximants and continued fractions is a profound connection in mathematical analysis, particularly for approximating functions like exp(x). Padé approximants, which are rational functions that match the Taylor series of a function up to a specified order, can often be expressed as continued fractions. This representation is advantageous because continued fractions can provide better convergence properties for certain functions, especially near singularities. For instance, the near-diagonal Padé approximants for exp(x), as shown above, can be systematically converted into a sequence of truncated continued fractions, revealing a structured pattern known as the “main Padé sequence.”

    A function with a convergent Taylor series around x = a can be approximated by a sequence of diagonal Padé approximants P(n,n)(x) = \frac{P_n(x)}{Q_n(x)}, provided the associated Hankel matrices H_n, built from the Taylor coefficients, have nonzero determinants.

    When this “Padé table” is normal (i.e., \det(H_n) \ne 0 for all n), each Padé approximant is uniquely defined.

    Under these conditions, one can systematically derive a continued fraction whose successive convergents exactly match the Padé approximants.

    This establishes a rigorous connection between the Taylor series, Padé approximation, and continued fraction representation of the function.

  • Continued fractions I

    We can use Padé approximants to build a sequence of truncated continued fractions representing a given function. For example using P(1,1):

    P(1,1) = \frac{A_0 + A_1 x}{1 + B_1 x} = C_0 + C_1 x + C_2 x^2

    Solving the linear systems presented in post Computing Padé approximants sequentially leads to:

    C_1 B_1 = - C_2 \quad \Rightarrow \quad B_1 = -\frac{C_2}{C_1}
    \begin{pmatrix} C_0 & 0 \\ C_1 & C_0 \end{pmatrix} \begin{pmatrix} 1 \\ B_1 \end{pmatrix} = \begin{pmatrix} A_0 \\ A_1 \end{pmatrix}
    A_0 = C_0, \quad A_1 = C_1 + C_0 B_1

    Setting C_0 = 1:

    P(1,1) = \frac{1 + (C_1 + B_1) x}{1 + B_1 x}
    P(1,1) = \frac{1 + (C_1 - \frac{C_2}{C_1}) x}{1 - \frac{C_2}{C_1} x}
    P(1,1) = \frac{1 + (C_1 - \frac{C_2}{C_1}) x}{1 + (C_1 - \frac{C_2}{C_1}) x - C_1 x}
    P(1,1) = \frac{1}{1 - \frac{C_1 x}{1 + (C_1 - \frac{C_2}{C_1}) x}}
    P(1,1) = \frac{1}{1 - \frac{C_1 x}{1 - (\frac{C_2}{C_1} - C_1) x}}

    Now define:

    \frac{b_0}{1 - \frac{b_1 x}{1 - b_2 x}} = \frac{1}{1 - \frac{C_1 x}{1 - \left(\frac{C_2}{C_1} - C_1\right) x}}

    We have:

    b_0 = 1, \quad b_1 = C_1, \quad b_2 = \frac{C_2}{C_1} - C_1, \quad \ldots

  • Python code for Padé approximants

    In this post, we provide readers with Python code that enables the derivation of Padé approximants from series coefficients given as input to the algorithm in the form of a coefficient vector. It is also necessary to predefine the degree of the numerator and denominator.

    import sympy as sp

def pade_approximation_function(series, m, n):
    """
    Returns the Padé approximation [m/n] in the form of a rational function
    with coefficients simplified as fractions.
    
    Parameters:
    series : list - Coefficients of the Taylor series (e.g., [1, 0, -1/2, 0, 1/24, ...] for cos(x))
    m : int - Degree of the numerator
    n : int - Degree of the denominator
    
    Returns:
    sympy.Expr - Rational function P(x)/Q(x) with simplified fractions
    """
    # Check that there are enough terms in the series
    if len(series) < m + n + 1:
        raise ValueError(f"The series must contain at least {m + n + 1} terms for an [m/n] approximation")

    # Convert the series coefficients to simplified fractions
    series = [sp.simplify(sp.Rational(str(c))) for c in series]
    
    # Symbolic variable
    x = sp.Symbol('x')
    
    # Coefficients of the numerator P(x) = a0 + a1*x + ... + am*x^m
    a = [sp.Symbol(f'a{i}') for i in range(m + 1)]
    # Coefficients of the denominator Q(x) = 1 + b1*x + ... + bn*x^n (b0 = 1)
    b = [1] + [sp.Symbol(f'b{i}') for i in range(1, n + 1)]
    
    # Polynomials P(x) and Q(x)
    P = sum(ai * x**i for i, ai in enumerate(a))
    Q = sum(bi * x**i for i, bi in enumerate(b))
    
    # The truncated series in polynomial form
    S = sum(c * x**i for i, c in enumerate(series))
    
    # Equation to solve: P(x) - Q(x)*S(x) = 0 up to order m+n
    expr = P - Q * S
    
    # Extract coefficients of x^0 to x^(m+n) and set the equations to 0
    equations = [sp.expand(expr).coeff(x, k) for k in range(m + n + 1)]
    
    # Variables to solve for (a0, a1, ..., am, b1, b2, ..., bn)
    unknowns = a + b[1:]
    
    # Solve the system
    solution = sp.solve(equations, unknowns)
    
    if not solution:
        raise ValueError("No solution found for this approximation")
    
    # Simplified coefficients for P(x) and Q(x)
    num_coeffs = [sp.simplify(sp.Rational(str(solution[ai]))) if ai in solution else 0 for ai in a]
    den_coeffs = [1] + [sp.simplify(sp.Rational(str(solution[bi]))) if bi in solution else 0 for bi in b[1:]]
    
    # Construct the final polynomials
    P_final = sum(coef * x**i for i, coef in enumerate(num_coeffs))
    Q_final = sum(coef * x**i for i, coef in enumerate(den_coeffs))
    
    # Return the simplified rational function
    return sp.simplify(P_final / Q_final)

    Here is the code to compute the Padé(2,2) for \cos(x)

        # Compute Padé(2,2) for cos(x)
    series = [1, 0, -sp.Rational(1,2), 0, sp.Rational(1,24)]
    m, n = 2, 2
    pade_approx = pade_approximation_function(series, m, n)
    print("Padé(2,2) approximant for cos(x):")
    sp.pprint(pade_approx)

  • Anharmonic oscillator I

    We would like to illustrate the use of Padé approximants in the context of the anharmonic oscillator. A harmonic oscillator is an oscillating system that experiences a restoring force. Anharmonicity is the deviation of a system from being a harmonic oscillator. The anharmonic oscillator is described by the following differential equation:

    \displaystyle \left(-\frac{d^2}{dx^2} + x^2 + x^4\right) \phi(x) = E_n \phi(x)

    This differential equation is very hard to solve exactly. In order to obtain an approximate solution, we can use perturbation theory. By inserting a small dimensionless parameter \epsilon, the equation becomes:

    \displaystyle \left(-\frac{d^2}{dx^2} + x^2 + \epsilon x^4\right) \phi(x) = E_n(\epsilon) \phi(x)

    If we are interested in the energy levels, the basic idea is to write E_n(\epsilon) as a geometric series. For the ground state E_0 we write:

    \displaystyle E_0(\epsilon) = \sum_{n=0}^{\infty} E_{0,n} \epsilon^n

    The solution is a divergent perturbation series (C.M. Bender and T.T. Wu, Phys. Rev. 184, 1969):

    \displaystyle E_0(\epsilon) = \frac{1}{2} + \frac{3}{4}\epsilon - \frac{21}{8}\epsilon^2 + \frac{333}{16}\epsilon^3 + \mathcal{O}(\epsilon^4)

    Let’s calculate the P(1,1) approximant of E_0 using the techniques presented in post Computing Padé approximants:

    \displaystyle P(1,1) = \frac{A_0 + A_1 \epsilon}{1 + B_1 \epsilon} = C_0 + C_1 \epsilon + C_2 \epsilon^2

    Solving the linear systems for computing Padé approximants sequentially leads to:

    \displaystyle C_1 B_1 = - C_2 \implies B_1 = -\frac{C_2}{C_1}
    \displaystyle \frac{3}{4} B_1 = \frac{21}{8} \implies B_1 = \frac{84}{24} = \frac{7}{2}
    \displaystyle \begin{pmatrix} C_0 & 0 \\ C_1 & C_0 \\ \end{pmatrix} \begin{pmatrix} 1 \\ B_1 \end{pmatrix} = \begin{pmatrix} A_0 \\ A_1 \end{pmatrix}
    \displaystyle \begin{pmatrix} \frac{1}{2} & 0 \\ \frac{3}{4} & \frac{1}{2} \\ \end{pmatrix} \begin{pmatrix} 1 \\ \frac{7}{2} \end{pmatrix} = \begin{pmatrix} \frac{1}{2} \\ \frac{3}{4} + \frac{7}{4} \end{pmatrix}
    \displaystyle P(1,1)_{E_0} = \frac{\frac{1}{2} + \left(\frac{3}{4} + \frac{7}{4}\right)\epsilon}{1 + \frac{7}{2}\epsilon}

    Setting \epsilon = 1 (to recover the initial differential equation) we have:

    \displaystyle P(1,1)_{E_0} = \frac{\frac{1}{2} + \frac{3}{4} + \frac{7}{4}}{1 + \frac{7}{2}} = 0.6666

    The P(4,4)_{E_0} = 1.3838 and the exact solution is 1.3924.

    The relative error (definition here) of the approximant P(4,4)_{E_0} is:

    \displaystyle \text{Relative error} = \frac{1.3838 - 1.3924}{1.3924} = -0.0062

    It’s therefore interesting to note that in order to solve a very difficult differential equation, we can use perturbation theory. In the case of the anharmonic oscillator, the perturbative series is divergent. Using Padé approximants, we can make use of a divergent perturbative series and approximate the exact answer arbitrarily.