--- jupytext: encoding: '# -*- coding: utf-8 -*-' formats: md:myst,ipynb text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.16.4 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- ```{code-cell} :tags: [hide-input] import mmf_setup; mmf_setup.nbinit() import os, sys from pathlib import Path FIG_DIR = Path(mmf_setup.ROOT) / '../Docs/_build/figures/' os.makedirs(FIG_DIR, exist_ok=True) import logging; logging.getLogger("matplotlib").setLevel(logging.CRITICAL) %matplotlib inline import numpy as np, matplotlib.pyplot as plt try: from myst_nb import glue except: glue = None from matplotlib.animation import FuncAnimation from phys_581.contexts import FPS ``` (sec:SequenceAcceleration)= # Sequence Acceleration Techniques As a sample problem, consider accurately computing the Riemann zeta function $\zeta(s)$ for values of $s$ close to $s\approx 1$: \begin{gather*} \zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^s}. \end{gather*} :::{doit} Spend a bit of time trying to implement this so you get an idea of how slowly the series converges before continuing. Try to estimate how many terms you need to achieve a specified tolerance. *You did remember to sum from smallest to largest, didn't you?* ::: :::{admonition} Estimate for convergence. For real $s$, the terms are monotonically decreasing, so we can bound the truncation error with integrals. *(Draw a picture if needed to convince yourself.)* \begin{gather*} \zeta_{N}(s) = \sum_{n=1}^{N-1} \frac{1}{n^s}, \\ \frac{1}{(s-1)N^{s-1}} = \int_{N}^{\infty} \frac{1}{n^s}\d{n} < \underbrace{\sum_{n=N}^{\infty} \frac{1}{n^s}}_{\zeta(s) - \zeta_N(s)} < \int_{N}^{\infty} \frac{1}{n^s}\d{n} + \frac{1}{N^s} = \frac{1}{(s-1)N^{s-1}} + \frac{1}{N^s},\\ \zeta(s) = \zeta_{N}(s) + \frac{1}{(s-1)N^{s-1}} + \frac{1}{2N^{s}} \pm \frac{1}{2N^{s}}. \end{gather*} Simply put, the error in the truncated sum can be estimated by the integral of the tail: \begin{gather*} \zeta(s) = \zeta_{N}(s) \pm \int_{N}^{\infty}\frac{1}{n^s}\d{n} = \zeta(s) \pm \frac{1}{(s-1)N^{s-1}}, \end{gather*} but we can improve the convergence by including this tail. To achieve a tolerance of $\epsilon$, one thus needs about \begin{gather*} N \gtrapprox \left(\frac{s-1}{\epsilon}\right)^{\frac{1}{s-1}}. \end{gather*} I.e., for 12 digits we need about $N\gtrapprox 10^{12}$ terms for $\zeta(2)$ (doable, but slow) and $N \gtrapprox 10^{110}$ terms for $\zeta(1.1)$ (impossible). Using the integral-corrected formula, the one needs \begin{gather*} N \gtrapprox \left(\frac{1}{2\epsilon}\right)^{\frac{1}{s}}. \end{gather*} I.e., for 12 digits we need about $N\gtrapprox 7\times 10^{5}$ terms for $\zeta(2)$, and $N \gtrapprox 4\times 10^{10}$ terms for $\zeta(1.1)$. ::: ```{code-cell} def zeta1(N, s=2): n = np.arange(1, N)[::-1] return (1/n**s).sum() def zeta2(N, s=2): return zeta1(N, s=s) + 1/(s-1)/N**(s-1) + 1/2/N**s zeta1(1000) - np.pi**2/6 ``` ## Simple Acceleration An obvious thing to do is to use the integral approximation to try to accelerate the convergence. The idea is that \begin{gather*} \zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^{s}} \approx \int_1^{\infty}\frac{1}{n^{s}}\d{n} = \left.\frac{-1}{(s-1)n^{s-1}}\right|_{n=1}^{\infty} = \frac{1}{s-1}. \end{gather*} Subtracting this should cancel the leading order terms, and we can do this integral term by term to get a correction factor: \begin{gather*} \zeta(s) = \sum_{n=1}^{\infty} \frac{1}{n^{s}} - \underbrace{\int_1^{\infty}\frac{1}{n^s}\d{n}} _{\sum_{n=1}^{\infty}\int_{n}^{n+1}\frac{1}{n^s}\d{n}} + \int_1^{\infty}\frac{1}{n^s}\d{n}\\ = \sum_{n=1}^{\infty} \underbrace{\left( \frac{1}{n^{s}} + \frac{1}{1-s}\left(\frac{1}{n^{s-1}} - \frac{1}{(n+1)^{s-1}}\right) \right)}_{a_n} + \frac{1}{s-1}. \end{gather*} For large $n$, we have \begin{gather*} a_n = n^{-s} + \frac{n^{1-s}}{1-s}\Bigl(1 - (1+1/n)^{1-s}\Bigr)\\ = n^{-s} + \frac{n^{1-s}}{1-s}\Biggl( 1 - 1 - \frac{(1-s)}{n} - \frac{(1-s)(1-s-1)}{2!n^2} + \order(n^{-3}) \Biggr)\\ = \frac{s}{2n^{s+1}} + \order(n^{-s-2}). \end{gather*} Integrating, we can estimate the error to be $\sim 1/2n^{s}$. \begin{gather*} \end{gather*} ```{code-cell} from sympy import var, oo x, q = var('x, q', positive=True) s = 1+q n = 1/x a_n = 1/n**s + 1/(1-s)*(1/n**(s-1) - 1/(n+1)**(s-1)) a_n.series(x, 0, 2) ``` ```{code-cell} from functools import cache @cache def zeta_N(N, s=2): return (1/np.arange(1, N)[::-1]**s).sum() def zeta1_N(N, s=2): ns = np.arange(1, N)[::-1] return 1/(s-1) + (1/ns**s + 1/(1-s)*(1/ns**(s-1) - 1/(ns+1)**(s-1))).sum() def S(n): return zeta_N(n, s=2) def S1(n, S=S): return S(n+1) - (S(n+1) - S(n))**2/(S(n+1) - 2*S(n) + S(n-1)) def S2(n, S=S1): return S1(n, S=S1) s = 1.1 Ns = np.linspace(1, 1000) plt.semilogy(Ns, 1/(Ns)**(s-1) - 1/(Ns+1)**(s-1)) plt.semilogy(Ns, (s-1)/Ns**s) class Levin: """Levin trsnformation. This is a fairly direct translation of the algorithm given in Numerical Recipes. It provide a generator that computes successive terms. """ def __init__(self, get_a, n0=0, nmax=100, type='u', beta=1.0, eps=np.finfo(float).eps): """ Arguments --------- get_a : function Return the nth term of the sum starting at n=0. """ self.get_a = get_a self.nmax = nmax self.eps = eps self.S_n = 0 # Partial sums - unaccelerated self.lastval = 0 self.lasteps = 0.0 self.n = 0 self.ncv = 0 self.cnvgd = False self.small = np.finfo(float).tiny * 10.0 self.big = np.finfo(float).max self.numer = [] self.denom = [] self.type = type self.beta = beta def __iter__(self): while not self.cnvgd and self.n < self.nmax: n = self.n a_n = self.get_a(n) self.S_n += a_n if 'u' == self.type: w_n = (self.beta + n) * a_n elif 't' == self.type: w_n = a_n elif 'd' == self.type: w_n = self.get_a(n+1) elif 'v' == self.type: a_n1 = self.get_a(n+1) w_n = a_n * a_n1/(a_n - a_n1) term = 1.0/(self.beta + n) self.denom.append(term/w_n) self.numer.append(self.S_n*self.denom[n]) if n > 0: ratio = (self.beta + n - 1) * term for j in range(1, n+1): fact = (n - j - self.beta) * term self.numer[n - j] = self.numer[n - j + 1] - fact * self.numer[n - j] self.denom[n - j] = self.denom[n - j + 1] - fact * self.denom[n - j] term *= ratio self.n += 1 val = self.lastval if abs(self.denom[0]) < self.small else self.numer[0]/self.denom[0] self.lasteps = abs(val - self.lastval) if (self.lasteps <= self.eps): self.ncv += 1 if (self.ncv >= 2): self.cnvgd = True self.lastval = val yield val def get_a(n, s=2.0): return zeta_N(2**(n+1), s=s) - zeta_N(2**n, s=s) def get_a(n, s=2.0): return l = np.array(list(Levin(get_a, eps=1e-12))) ``` ## Richardson Extrapolation Consider the partial sums $S_{N}$ and remainders $E_{N} = S-S_{N}$. As discussed above, the error can be expressed as a polynomial in $1/N$. More generally, we can write: \begin{gather*} S = S_{N} + \sum_{n=0}^{\infty}\frac{a_{n}}{N^{k_n}}. \end{gather*} In our case, $k_{n} = n+1$. The idea is to use $S_{2N}$ to remove the first error term, improving the convergence: \begin{align*} S &= S_{N} + \frac{a_0}{N^{k_0}} + \frac{a_1}{N^{k_1}} + \cdots,\\ &= S_{2N} + \frac{a_0}{2^{k_0}N^{k_0}} + \frac{a_1}{2^{k_1}N^{k_1}} + \cdots,\\ (1 - 2^{k_0})S &= S_{N} - 2^{k_0}S_{2N} + \frac{a_1}{2^{k_1-k_0}N^{k_1}} + \cdots,\\ S &= \underbrace{\frac{S_{N} - 2^{k_0}S_{2N}}{(1 - 2^{k_0})}}_{S^{(1)}_{N}} + \frac{b_1}{N^{k_1}} + \cdots. \end{align*} We now have a new approximation $S^{(1)}_{N}$ that converges more quickly (as $1/N^{k_1}$). The idea is to repeat the process on $S^{(1)}_{N}$ to obtain a new series $S^{(2)}_{N}$ the converges like $1/N^{k_2}$, and so on. This gives us a recursive formula: \begin{gather*} S^{(n+1)}_{N} = \frac{S^{(n)}_{N} - 2^{k_n}S^{(n)}_{2N}}{(1 - 2^{k_n})}. \end{gather*} We can arrange these results in a triangular table: $s_{n,p} = S^{(n)}_{2^p}$. This defines the following recurrence: \begin{align*} s_{0, p} &= S_{2^{p}}, \\ s_{n, p} &= \frac{s_{n-1, p} - 2^{k_{n-1}}s_{n-1, p+1}}{1 - 2^{k_{n-1}}}. \end{align*} The key points are that, to form the new series, we only need to know $k_{n}$: we do **not** need to know the values of $a_{n}$. We also do not need to know any special values like $\zeta(2)$ etc. We use the series itself to perform the corrections. The downside is that we need to double the number of terms at each step, but we can start with $N=1$, and often this converges very quickly. Another potential problem is that you need to properly characterize $k_{n}$ for this to work. ::: ```{code-cell} :tags: [hide-input, margin] import numpy as np from functools import cache def get_S0(N): """Naive sum.""" # To prevent round-off error, we sum from small to large. n = np.arange(0, N+1)[::-1] return sum(1/(2*n**2 + 1)) @cache # We use the cache decorator to memoize def s(n, p): if n == 0: return get_S0(2**p) k = (n-1)+1 return (s(n-1, p) - 2**k * s(n-1, p+1))/(1 - 2**k) # To decide when to terminate, we look at successive terms in the table. # This can fail, so be careful. n = 1 while abs(s(n, 0) - s(n-1, 1)) > 1e-16: n += 1 S = s(n, 0) print(f"S = s({n}, 0) = {S}") # In this case, the exact answer is known... S_ = (2 + np.sqrt(2)*np.pi / np.tanh(np.pi/np.sqrt(2)))/4 assert abs(S-S_) < 1e-15 # Print out the table nmax = 11 ns = np.arange(nmax) print(" n: p=" + " ".join([f"{p:2}" for p in range(nmax)])) for n in ns: print(f"{n:2}: " + ",".join([f"{int(np.floor(-np.log10(abs(s(n, p)-S)+1e-16))):2}" for p in range(nmax-n)])) ``` :::{margin} **Richardson Extrapolation:** Number of correct digits in $s_{n, p}$. Note that 16 digits of accuracy are obtained for $s_{10, 0}$ which requires evaluating only $N = 2^{10} = 1024$ terms. For a termination criterion, we set $\abs{s_{n, 0} - s_{n-1, 1}} < 10^{-16}$. ::: ```{code-cell} def get_S0(N): """Naive sum.""" # To prevent round-off error, we sum from small to large. n = np.arange(0, N+1)[::-1] return sum(1/(2*n**2 + 1)) def get_S1(N): """Improved convergence.""" n = np.arange(1, N+1)[::-1] return 1 + 0.5 + sum(1/(2*n**2 + 1) - 1/(2*n*(n+1))) def get_S2(N): """Improved with zeta(2).""" n = np.arange(1, N+1)[::-1] zeta2 = np.pi**2/6 return 1 + zeta2/2 + sum(1/(2*n**2 + 1) - 1/2/n**2) Ns = 2**np.arange(22) S0s = np.array(list(map(get_S0, Ns))) S1s = np.array(list(map(get_S1, Ns))) S2s = np.array(list(map(get_S2, Ns))) S = (2 + np.sqrt(2)*np.pi / np.tanh(np.pi/np.sqrt(2)))/4 fig, ax = plt.subplots() ax.loglog(Ns, abs(S-S0s), '-+', label='Naive') ax.loglog(Ns, abs(S-S1s), '-o', label='Improved') ax.loglog(Ns, abs(S-S2s), ':+', label='$\zeta(2)$ Improved') ps = np.arange(11) ax.loglog(2**ps, [abs(S-s(p, 0)) for p in ps], ':.', label="Richardson extrapolation") ax.legend() ax.set(xlabel="$N$", ylabel="$S$") print(S1s[-3:]) ``` * **(f)** \begin{gather*} \sum_{n=1}^{\infty}\frac{(-1)^{n}}{n}. \end{gather*}