Approximaths

Approximaths

  • Richardson extrapolation: example

    Richardson extrapolation acts as a powerful convergence accelerator in numerical methods. Instead of relying solely on extremely fine discretizations—which can be computationally expensive—it combines results from coarser simulations to estimate the zero-step (continuum) limit.

    This idea appears in several numerical techniques (e.g., Romberg integration, mesh convergence studies in engineering), where it can significantly improve accuracy at a modest additional cost.

    When the error behaves regularly, Richardson extrapolation can dramatically reduce the error—sometimes by an order of magnitude or more—without requiring prohibitively fine meshes.

    Consider the following (slowly) converging series (see this post Zeta summation: Basel problem):

    \displaystyle 1 + \frac{1}{4} + \frac{1}{9} + \frac{1}{16} + ... = \frac{\pi^2}{6}

    and calculate its corresponding Richardson extrapolation using up to 4-term partial sums:

    \displaystyle N = 2
    \displaystyle S_2 = 1 + \frac{1}{4} = \frac{5}{4}
    \displaystyle S_3 = 1 + \frac{1}{4} + \frac{1}{9} = \frac{49}{36}
    \displaystyle S_4 = 1 + \frac{1}{4} + \frac{1}{9} + \frac{1}{16} = \frac{205}{144}

    The Richardson extrapolation is (see this post Richardson extrapolation):

    \displaystyle \mathcal{R}_N = \frac{N^2S_N -2(N+1)^2S_{N+1} + (N+2)^2S_{N+2}}{2}

    in this case:

    \displaystyle \mathcal{R}_2 = \frac{4 \times (\frac{5}{4}) - 2 \times 9 \times (\frac{49}{36}) + 16 \times \frac{205}{144}}{2}
    \displaystyle = \frac{5 - \frac{49}{2} + \frac{205}{9}}{2}
    \displaystyle = 1.6389

    We have:

    \displaystyle \frac{\pi^2}{6} = 1.6449
    \displaystyle \mathcal{R}_{2} = 1.6389
    \displaystyle S_4 = 1.4236

    Using solely a partial series with 4 terms, we were able to speed up the convergence of the above-mentioned series considerably using Richardson extrapolation.

  • Richardson extrapolation

    Another method for accelerating convergence is Richardson extrapolation. This method is similar to the Shanks transformation. Let’s model the partial sum as:

    \displaystyle S_N = \sum_{n=0}^{N} a_n = S + \frac{a}{N} + \frac{b}{N^2} + \frac{c}{N^3} + \dots

    where S = \lim_{N \to \infty} S_N.

    \displaystyle S_N = S + \frac{a}{N} + \frac{b}{N^2} + \frac{c}{N^3} + \dots
    \displaystyle S_{N+1} = S + \frac{a}{N+1} + \frac{b}{(N+1)^2} + \frac{c}{(N+1)^3} + \dots
    \displaystyle S_{N+2} = S + \frac{a}{N+2} + \frac{b}{(N+2)^2} + \frac{c}{(N+2)^3} + \dots

    equivalently:

    \displaystyle N^2 S_N = N^2 S + N a + b + \frac{c}{N} + \dots
    \displaystyle (N+1)^2 S_{N+1} = (N+1)^2 S + (N+1) a + b + \frac{c}{N+1} + \dots
    \displaystyle (N+2)^2 S_{N+2} = (N+2)^2 S + (N+2) a + b + \frac{c}{N+2} + \dots

    and

    \displaystyle N^2 S_N = N^2 S + N a + b + \frac{c}{N} + \dots
    \displaystyle -2 (N+1)^2 S_{N+1} = -2 (N+1)^2 S - 2 (N+1) a - 2 b - 2 \frac{c}{N+1} + \dots
    \displaystyle (N+2)^2 S_{N+2} = (N+2)^2 S + (N+2) a + b + \frac{c}{N+2} + \dots

    Summing up the system above:

    \displaystyle N^2 S_N - 2 (N+1)^2 S_{N+1} + (N+2)^2 S_{N+2} = N^2 S + N a + b + \frac{c}{N} + \dots
    \displaystyle -2 (N+1)^2 S - 2 (N+1) a - 2 b - 2 \frac{c}{N+1} + \dots
    \displaystyle + (N+2)^2 S + (N+2) a + b + \frac{c}{N+2} + \dots

    Simplifying and taking the limit N \to \infty:

    \displaystyle N^2 S_N - 2 (N+1)^2 S_{N+1} + (N+2)^2 S_{N+2} = N^2 S + a N + b
    \displaystyle - 2 N^2 S - 4 N S - 2 S - 2 a N - 2 a - 2 b
    \displaystyle + N^2 S + 4 N S + 4 S + a N + 2 a + b
    \displaystyle = 2 S
    \displaystyle \implies S = \frac{N^2 S_N - 2 (N+1)^2 S_{N+1} + (N+2)^2 S_{N+2}}{2} \quad \text{(as } N \to \infty\text{)}

    On this basis, we define the Richardson extrapolation as:

    \boxed{\displaystyle \mathcal{R}_N := \frac{N^2 S_N - 2 (N+1)^2 S_{N+1} + (N+2)^2 S_{N+2}}{2}}

  • Shanks transformation: Example 2

    To illustrate the Shanks transformation we will use the following series:

    \displaystyle \ln(x) = (x-1)-\frac{1}{2}(x-1)^2+\frac{1}{3}(x-1)^3-\dots
    \displaystyle \ln(x) = \sum_{n=1}^{\infty} (-1)^{n+1} \frac{(x-1)^n}{n}

    This series is converging very slowly. To illustrate this we will compute ln(2):

    \displaystyle \ln(2) = 1 - \frac{1}{2} + \frac{1}{3} - \frac{1}{4} + \frac{1}{5} + \dots
    \displaystyle \ln(2) = \sum_{n=1}^{\infty} (-1)^{n+1} \frac{1}{n}

    The corresponding partial sum of ln(2) is:

    \displaystyle S_N = \sum_{n=1}^{N} (-1)^{n+1} \frac{1}{n}

    Let’s calculate some terms of this partial sum:

    \displaystyle S_1 = 1 \quad S_{10} = 0.6456
    \displaystyle S_2 = 0.5 \quad S_{400} = 0.6919
    \displaystyle S_3 = 0.8333 \quad \ln(2) = 0.6931

    Perform the Shanks transform (using SN = S4):

    \displaystyle S_4 = 1 - \frac{1}{2} + \frac{1}{3} - \frac{1}{4} = \frac{12}{12} - \frac{6}{12} + \frac{4}{12} - \frac{3}{12} = \frac{7}{12}
    \displaystyle S_4^2 = \frac{49}{144}
    \displaystyle S_5 = \frac{7}{12} + \frac{1}{5} = \frac{35}{60} + \frac{12}{60} = \frac{47}{60}
    \displaystyle S_3 = 1 -\frac{1}{2} + \frac{1}{3} = \frac{6}{6} - \frac{3}{6} + \frac{2}{6} = \frac{5}{6}
    \displaystyle \mathcal{S} = \frac{S_N^2 - S_{N+1}S_{N-1}}{2S_N -S_{N+1}- S_{N-1}}
    \displaystyle \mathcal{S} = \frac{\frac{49}{144} - (\frac{47}{60} \frac{5}{6})}{\frac{14}{12} - \frac{47}{60} -\frac{5}{6}}
    \displaystyle \mathcal{S} = \frac{\frac{49}{144} - \frac{235}{360}}{\frac{70}{60} - \frac{47}{60} -\frac{50}{60}}
    \displaystyle \mathcal{S} = \frac{-225/720}{-27/60}
    \displaystyle \mathcal{S} = 0.6944

    Since S_{358} = 0.6918 and \ln(2) = 0.6931 It is necessary to calculate up to 358 terms of the partial sum SN to obtain a relative error equal to that of the Shanks transformation \mathcal{S} based on the first 4 terms.

  • Shanks transformation: Example 1

    Let’s transform the geometric series:

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

    using the corresponding partial sum S2:

    \displaystyle 1+x+x^2
    \displaystyle S_2^2 = (1+x+x^2)^2 = 1 + 2x + 3x^2 + 2x^3 + x^4
    \displaystyle S_3 S_1 = (1+x+x^2+x^3)(1+x) = 1 + 2x + 2x^2 + 2x^3 + x^4
    \displaystyle 2S_2 = 2 + 2x + 2x^2
    \displaystyle S_3 = 1+x+x^2+x^3
    \displaystyle S_1 = 1 + x
    \displaystyle \mathcal{S} = \frac{(1+2x+3x^2+2x^3+x^4) - (1+2x+2x^2+2x^3+x^4)}{(2+2x +2x^2) - (1+x+x^2+x^3) - (1+x)}
    \displaystyle \mathcal{S} = \frac{x^2}{x^2-x^3}
    \displaystyle \mathcal{S} = \frac{1}{1-x}

    Therefore the Shanks transformation of the geometric series is correctly \frac{1}{1-x}. A similar result can be obtained for the series:

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

    The corresponding Shanks transformation is \frac{1}{1+x}.

    Remarkably, the Shanks transformation yields the exact sum of the geometric series using only a few terms, since the error model S_N = S + \alpha \beta^N perfectly describes such a series. The Gregory series for \pi is notoriously slow, requiring five hundred terms just to calculate two correct decimal places. However, by applying the Shanks transformation, we can accelerate this convergence dramatically, extracting high-precision results from only a few initial partial sums. This demonstration highlights how a simple mathematical shift can transform an inefficient infinite series into a powerful computational tool.

  • Accelerate convergence (Shanks Transformation)

    Imagine you’ve solved a problem using a Taylor series expansion of the solution, and the resulting series converges very slowly. The Shanks transformation is one way of improving the speed of convergence of a convergent series.
    Given the partial sum:

    \displaystyle S_N = a_0 + a_1 + a_2 + \dots + a_N

    Suppose that this partial sum converges to:

    \displaystyle S = a_0 + a_1 + a_2 + \dots
    \displaystyle S = \lim_{N \to \infty} \sum_{k=0}^{N} a_k

    We can model the difference between S and SN as follows:

    \displaystyle S_N = S + \alpha \beta^{N}

    Where |\beta| < 1:

    \displaystyle S_{N+1} = S + \alpha \beta^{N+1}
    \displaystyle S_N = S + \alpha \beta^{N}
    \displaystyle S_{N-1} = S + \alpha \beta^{N-1}

    and

    \displaystyle S_{N+1} - S = \alpha \beta^{N+1}
    \displaystyle S_N - S = \alpha \beta^{N}
    \displaystyle S_{N-1} - S = \alpha \beta^{N-1}

    Therefore:

    \displaystyle \frac{S_N - S}{S_{N-1}-S} = \frac{\alpha \beta^{N}}{\alpha \beta^{N-1}} = \beta
    \displaystyle \frac{S_{N+1} - S}{S_{N}-S} = \frac{\alpha \beta^{N+1}}{\alpha \beta^{N}} = \beta
    \displaystyle \frac{S_N - S}{S_{N-1}-S} = \frac{S_{N+1} - S}{S_{N}-S}
    \displaystyle (S_N - S)(S_{N}-S) = (S_{N+1} - S)(S_{N-1}-S)
    \displaystyle S_N^2 - 2S_NS + S^2 = S_{N+1}S_{N-1}- SS_{N+1} - SS_{N-1} + S^2
    \displaystyle S_N^2 - 2S_NS + S^2= S_{N+1}S_{N-1} + S^{2}- S(S_{N+1}+ S_{N-1})
    \displaystyle S_N^2 - 2S_NS = S_{N+1}S_{N-1}- S(S_{N+1}+ S_{N-1})
    \displaystyle S_N^2 - S_{N+1}S_{N-1} = - S(S_{N+1}+ S_{N-1}) - 2S_NS
    \displaystyle S_N^2 - S_{N+1}S_{N-1} = S(-S_{N+1}- S_{N-1})+ 2S_N
    \displaystyle S = \frac{S_N^2 - S_{N+1}S_{N-1}}{2S_N -S_{N+1}- S_{N-1}}

    According to this result, the Shanks transformation is defined as:

    \displaystyle \mathcal{S} := \frac{S_N^2 - S_{N+1}S_{N-1}}{2S_N -S_{N+1}- S_{N-1}}

  • Summation summary

    Series Standard Euler Borel Borel–Écalle Zeta
    \displaystyle 1+\frac{1}{4}+\frac{1}{9}+\dots
    \displaystyle \frac{\pi^2}{6}
    \displaystyle \frac{\pi^2}{6}
    \displaystyle \frac{\pi^2}{6}
    \displaystyle \frac{\pi^2}{6}
    \displaystyle \frac{\pi^2}{6}
    \displaystyle 1-1+1-1+\dots
    \displaystyle \frac12
    \displaystyle \frac12
    \displaystyle \frac12
    \displaystyle \frac12
    \displaystyle 1!-2!+3!-4!+\dots
    \displaystyle \approx 0.5963
    \displaystyle \approx 0.5963
    \displaystyle 1+2+3+4+\dots
    \displaystyle -\frac{1}{12}
    \displaystyle 1!+2!+3!+4!+\dots

    The various summation methods presented in this table are not arbitrary assignments of values to divergent series; they are governed by the fundamental principle of regularity and consistency. A method is “regular” if it preserves the sum of any convergent series, and “consistent” if different methods yield the same value when they both apply to a given series.

    In theoretical physics and asymptotic analysis, these techniques are indispensable. Borel and Euler summations are frequently used to handle “weakly” divergent perturbative expansions in quantum mechanics. Zeta-function regularization is a cornerstone of modern physics, particularly in the calculation of the Casimir effect and in string theory, where it provides a rigorous way to subtract infinities from physical observables. Meanwhile, the Borel–Écalle theory of resurgent functions offers the most sophisticated framework for “decoding” the non-perturbative information hidden within divergent power series. By bridging the gap between divergent expansions and their underlying analytic structures, these methods allow mathematicians and physicists to extract precise physical predictions from seemingly ill-defined mathematical objects.

  • Classic Summation Axioms

    So far, we’ve seen several ways of summing series (the usual method, Euler, Borel, generic, Borel-Écalle and Zeta summation). All of these methods fulfill all three properties — except Zeta summation, which fulfills none of them.

    \displaystyle \mathcal{S}(a_0+a_1+a_2+ ...) = \lim_{n \to \infty} \sum_{k =0}^{n}a_k \quad\text{(Regularity)}
    \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{(Linearity)}
    \displaystyle \mathcal{S}(a_0+a_1+a_2+ ...) = a_0 + \mathcal{S}(a_1+a_2+ ...) \quad\text{(Stability)}

    Concerning the Zeta summation we should restrain our enthusiasm by saying that a certain meromorphic complex function called \zeta which have the value \zeta(s)=\sum_{n=1}^{\infty} \frac{1}{n^s} for all Re(s)>1 can be defined on the whole complex plane except at 1, in such a way that, \zeta(0) = -\frac{1}{2}. This ‘summation’ does not rely on any of the properties cited above. So it’s important to be clear about the method we are using and what properties it fulfills.

  • Zeta summation VI

    If we want to calculate \zeta(-2) for example we need to have a proper value for \Gamma(s/2) in for Re(s) \leq 0. The problem is that the Bernoulli representation of the Gamma function presented above:

    \displaystyle \Gamma(s) = \int_{0}^{\infty} e^{-t} t^{s-1}\,dt

    is only valid for Re(s) > 0. We have to provide another representation of the Gamma function for negative values of s:

    \Gamma(s) = \int_{0}^{\infty} e^{-t}\, t^{s-1}\, dt
    = \int_{1}^{\infty} e^{-t}\, t^{s-1}\, dt + \int_{0}^{1} e^{-t}\, t^{s-1}\, dt
    = \int_{1}^{\infty} e^{-t}\, t^{s-1}\, dt + \int_{0}^{1} \sum_{n=0}^{\infty} \frac{(-1)^n}{n!} t^{n+s-1}\, dt
    = \int_{1}^{\infty} e^{-t}\, t^{s-1}\, dt + \sum_{n=0}^{\infty} \frac{(-1)^n}{n!} \int_{0}^{1} t^{n+s-1}\, dt
    = \int_{1}^{\infty} e^{-t}\, t^{s-1}\, dt + \sum_{n=0}^{\infty} \frac{(-1)^n}{n!} \left[ \frac{t^{n+s}}{n+s} \right]_{0}^{1}
    = \int_{1}^{\infty} e^{-t}\, t^{s-1}\, dt + \sum_{n=0}^{\infty} \frac{(-1)^n}{n!\,(n+s)}

    The sum \sum_{n=0}^{\infty} \frac{(-1)^n}{n!(n+s)} exists for negative s except for s= -1,-2,-3,... This implies that \frac{1}{\Gamma(s/2)} is entire with simple zeros at s= -2,-4,-6,.... The simple pole of \xi(s) at zero is cancelled by the corresponding zero of \frac{1}{\Gamma(s/2)}. As a consequence the only singularity of \zeta(s) is a single pole at s=1. Since

    \displaystyle \zeta(s) = \pi^{s/2}\frac{\xi(s)}{\Gamma(s/2)}

    and \frac{1}{\Gamma(s/2)} is entire with simple zeros at s= -2,-4,-6,..., this implies that for example:

    \displaystyle \zeta(-2) = 0

    We have in general

    \displaystyle \zeta(-2n)=0

    The so-called “trivial zeros” of the Zeta function. They occur at the negative even integers.

  • Zeta summation V

    In a previous post we have computed the meromorphic continuation into the entire complex plane of the Zeta function:

    \displaystyle \zeta(s) = \pi^{s/2}\frac{\xi(s)}{\Gamma(s/2)}

    This analytic continuation allows us to sum the Zeta function in the whole complex plane. For special (‘known’) values of the Xi function we can calculate the Zeta function. For example for

    s = 2 we have \xi(2) = \frac{\pi}{6} and \Gamma(2/2) = \Gamma(1) = 1 and therefore:

    \displaystyle \zeta(2) = \frac{\pi^2}{6}
    \displaystyle \frac{1}{1^2} + \frac{1}{2^2} + \frac{1}{3^2} + \frac{1}{4^2} + ... = \frac{\pi^2}{6}

    For s = 4 we have \xi(4) = \frac{\pi^2}{90} and \Gamma(4/2) = \Gamma(2) = 1:

    \displaystyle \zeta(4) = \frac{\pi^4}{90}
    \displaystyle \frac{1}{1^4} + \frac{1}{2^4} + \frac{1}{3^4} + \frac{1}{4^4} + ... = \frac{\pi^4}{90}

    For s = 6 we have \xi(6) = \frac{2\pi^3}{945} and \Gamma(6/2) = \Gamma(3) = 2:

    \displaystyle \zeta(6) = \frac{\pi^6}{945}
    \displaystyle \frac{1}{1^6} + \frac{1}{2^6} + \frac{1}{3^6} + \frac{1}{4^6} + ... = \frac{\pi^6}{945}

    So, Zeta of all the positive even integers has the form

    \displaystyle\zeta(2n) = (-1)^{n+1} \frac{B_{2n} (2\pi)^{2n}}{2(2n)!}

    Where B_{2n} are Bernoulli numbers.

  • 1 + 2 + 3 + 4 + …

    To calculate the value of the Riemann zeta function at s = -1, we use the functional equation (presented here):

    \displaystyle \zeta(s) = \pi^{s/2} \frac{\xi(s)}{\Gamma(s/2)}

    The Dirichlet series for the Riemann zeta function is defined for \Re(s) > 1 as:

    \displaystyle \zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}

    For s = -1, the formal series (which is divergent in the classical sense) is:

    \displaystyle \zeta(-1) = \sum_{n=1}^{\infty} \frac{1}{n^{-1}}= \sum_{n=1}^{\infty} n = 1 + 2 + 3 + 4 + \dots

    Substituting s = -1 into the formula, we obtain:

    \displaystyle \zeta(-1) = \pi^{-1/2} \frac{\xi(-1)}{\Gamma(-1/2)}

    Using the recurrence property of the Gamma function, \Gamma(z) = \frac{\Gamma(z+1)}{z}, and the fact that \Gamma(1/2) = \sqrt{\pi}:

    \displaystyle \Gamma(-1/2) = \frac{\Gamma(-1/2+1)}{-1/2} = \frac{\Gamma(1/2)}{-1/2} = -2\sqrt{\pi}

    We exploit the functional equation \xi(s) = \xi(1-s), which implies:

    \displaystyle \xi(-1) = \xi(1 - (-1)) = \xi(2)

    From the definition of \xi(s), we have:

    \displaystyle \xi(2) = \pi^{-2/2} \Gamma(2/2) \zeta(2)

    Using \Gamma(1) = 1 and the solution of the Basel problem presented in the previous post \zeta(2) = \frac{\pi^2}{6}:

    \displaystyle \xi(-1) = \pi^{-1} \cdot 1 \cdot \frac{\pi^2}{6} = \frac{\pi}{6}

    Combining these results into our original expression:

    \displaystyle \zeta(-1) =\frac{1}{\sqrt{\pi}} \cdot \frac{\frac{\pi}{6}}{-2\sqrt{\pi}}
    \displaystyle = \frac{\pi}{-12 \cdot (\sqrt{\pi} \cdot \sqrt{\pi})}
    \displaystyle = \frac{\pi}{-12\pi}
    \displaystyle \boxed{\zeta(-1) = -\frac{1}{12}}

    Though surprising at first, turns out to be extremely useful in physics — most notably in string theory and the calculation of the Casimir effect, where the infinite sum 1 + 2 + 3 + 4 + … = −1/12 naturally appears when regularizing the vacuum energy between two conducting plates.