This cell adds /home/docs/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/checkouts/latest/src to your path, and contains some definitions for equations and some CSS for styling the notebook. If things look a bit strange, please try the following:
- Choose "Trust Notebook" from the "File" menu.
- Re-execute this cell.
- Reload the notebook.
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*}\]
Do It!
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?
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)\).
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
np.float64(-0.0010005001666666402)
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*}\]
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)
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/mul.py:1980, in Mul._eval_nseries(self, x, n, logx, cdir)
1979 for t in self.args:
-> 1980 coeff, exp = t.leadterm(x)
1981 if not coeff.has(x):
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:3558, in Expr.leadterm(self, x, logx, cdir)
3557 from sympy.functions.elementary.exponential import log
-> 3558 l = self.as_leading_term(x, logx=logx, cdir=cdir)
3559 d = Dummy('logx')
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
71 try:
---> 72 retval = cfunc(*args, **kwargs)
73 except TypeError as e:
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:3521, in Expr.as_leading_term(self, logx, cdir, *symbols)
3520 return self
-> 3521 obj = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
3522 if obj is not None:
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/add.py:1076, in Add._eval_as_leading_term(self, x, logx, cdir)
1075 while res.is_Order:
-> 1076 res = old._eval_nseries(x, n=n0+incr, logx=logx, cdir=cdir).cancel().powsimp().trigsimp()
1077 incr *= 2
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/add.py:507, in Add._eval_nseries(self, x, n, logx, cdir)
506 def _eval_nseries(self, x, n, logx, cdir=0):
--> 507 terms = [t.nseries(x, n=n, logx=logx, cdir=cdir) for t in self.args]
508 return self.func(*terms)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:3466, in Expr.nseries(self, x, x0, n, dir, logx, cdir)
3465 else:
-> 3466 return self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/mul.py:2046, in Mul._eval_nseries(self, x, n, logx, cdir)
2045 if res != self:
-> 2046 if (self - res).subs(x, 0) == S.Zero and n > 0:
2047 lt = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/decorators.py:248, in _SympifyWrapper.make_wrapped.<locals>._func(self, other)
247 return retval
--> 248 return func(self, other)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/decorators.py:118, in call_highest_priority.<locals>.priority_decorator.<locals>.binary_op_wrapper(self, other)
117 return f(self)
--> 118 return func(self, other)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:244, in Expr.__sub__(self, other)
241 @sympify_return([('other', 'Expr')], NotImplemented)
242 @call_highest_priority('__rsub__')
243 def __sub__(self, other) -> Expr:
--> 244 return Add(self, -other)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
71 try:
---> 72 retval = cfunc(*args, **kwargs)
73 except TypeError as e:
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/operations.py:108, in AssocOp.__new__(cls, evaluate, _sympify, *args)
106 return args[0]
--> 108 c_part, nc_part, order_symbols = cls.flatten(args)
109 is_commutative = not nc_part
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/add.py:382, in Add.flatten(cls, seq)
380 for t in newseq:
381 # x + O(x) -> O(x)
--> 382 if not any(o.contains(t) for o in order_factors):
383 newseq2.append(t)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/add.py:382, in <genexpr>(.0)
380 for t in newseq:
381 # x + O(x) -> O(x)
--> 382 if not any(o.contains(t) for o in order_factors):
383 newseq2.append(t)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
71 try:
---> 72 retval = cfunc(*args, **kwargs)
73 except TypeError as e:
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/series/order.py:453, in Order.contains(self, expr)
452 obj = self.func(expr, *self.args[1:])
--> 453 return self.contains(obj)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
71 try:
---> 72 retval = cfunc(*args, **kwargs)
73 except TypeError as e:
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/series/order.py:428, in Order.contains(self, expr)
427 from sympy.series.limits import Limit
--> 428 l = Limit(ratio, s, point).doit(heuristics=False)
429 if not isinstance(l, Limit):
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/series/limits.py:372, in Limit.doit(self, **hints)
371 else:
--> 372 raise NotImplementedError("Not sure of sign of %s" % ex)
374 # gruntz fails on factorials but works with the gamma function
375 # If no factorial term is present, e should remain unchanged.
376 # factorial is defined to be zero for negative inputs (which
377 # differs from gamma) so only rewrite for non-negative z0.
NotImplementedError: Not sure of sign of 2 - q
During handling of the above exception, another exception occurred:
NotImplementedError Traceback (most recent call last)
Cell In[3], line 6
4 n = 1/x
5 a_n = 1/n**s + 1/(1-s)*(1/n**(s-1) - 1/(n+1)**(s-1))
----> 6 a_n.series(x, 0, 2)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:3092, in Expr.series(self, x, x0, n, dir, logx, cdir)
3090 from sympy.series.order import Order
3091 if n is not None: # nseries handling
-> 3092 s1 = self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
3093 o = s1.getO() or S.Zero
3094 if o:
3095 # make sure the requested order is returned
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/add.py:507, in Add._eval_nseries(self, x, n, logx, cdir)
506 def _eval_nseries(self, x, n, logx, cdir=0):
--> 507 terms = [t.nseries(x, n=n, logx=logx, cdir=cdir) for t in self.args]
508 return self.func(*terms)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:3466, in Expr.nseries(self, x, x0, n, dir, logx, cdir)
3464 return self.series(x, x0, n, dir, cdir=cdir)
3465 else:
-> 3466 return self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/mul.py:2004, in Mul._eval_nseries(self, x, n, logx, cdir)
2002 if n0.is_nonnegative:
2003 n0 = S.Zero
-> 2004 facs = [t.nseries(x, n=ceiling(n-n0), logx=logx, cdir=cdir) for t in self.args]
2005 from sympy.simplify.powsimp import powsimp
2006 res = powsimp(self.func(*facs).expand(), combine='exp', deep=True)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:3466, in Expr.nseries(self, x, x0, n, dir, logx, cdir)
3464 return self.series(x, x0, n, dir, cdir=cdir)
3465 else:
-> 3466 return self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/add.py:507, in Add._eval_nseries(self, x, n, logx, cdir)
506 def _eval_nseries(self, x, n, logx, cdir=0):
--> 507 terms = [t.nseries(x, n=n, logx=logx, cdir=cdir) for t in self.args]
508 return self.func(*terms)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:3466, in Expr.nseries(self, x, x0, n, dir, logx, cdir)
3464 return self.series(x, x0, n, dir, cdir=cdir)
3465 else:
-> 3466 return self._eval_nseries(x, n=n, logx=logx, cdir=cdir)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/mul.py:2046, in Mul._eval_nseries(self, x, n, logx, cdir)
2043 return res
2045 if res != self:
-> 2046 if (self - res).subs(x, 0) == S.Zero and n > 0:
2047 lt = self._eval_as_leading_term(x, logx=logx, cdir=cdir)
2048 if lt == S.Zero:
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/decorators.py:248, in _SympifyWrapper.make_wrapped.<locals>._func(self, other)
246 if not isinstance(other, expectedcls):
247 return retval
--> 248 return func(self, other)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/decorators.py:118, in call_highest_priority.<locals>.priority_decorator.<locals>.binary_op_wrapper(self, other)
116 if f is not None:
117 return f(self)
--> 118 return func(self, other)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/expr.py:244, in Expr.__sub__(self, other)
241 @sympify_return([('other', 'Expr')], NotImplemented)
242 @call_highest_priority('__rsub__')
243 def __sub__(self, other) -> Expr:
--> 244 return Add(self, -other)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
69 @wraps(func)
70 def wrapper(*args, **kwargs):
71 try:
---> 72 retval = cfunc(*args, **kwargs)
73 except TypeError as e:
74 if not e.args or not e.args[0].startswith('unhashable type:'):
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/operations.py:108, in AssocOp.__new__(cls, evaluate, _sympify, *args)
105 if len(args) == 1:
106 return args[0]
--> 108 c_part, nc_part, order_symbols = cls.flatten(args)
109 is_commutative = not nc_part
110 obj = cls._from_args(c_part + nc_part, is_commutative)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/add.py:382, in Add.flatten(cls, seq)
379 newseq2 = []
380 for t in newseq:
381 # x + O(x) -> O(x)
--> 382 if not any(o.contains(t) for o in order_factors):
383 newseq2.append(t)
384 newseq = newseq2 + order_factors # type: ignore
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/add.py:382, in <genexpr>(.0)
379 newseq2 = []
380 for t in newseq:
381 # x + O(x) -> O(x)
--> 382 if not any(o.contains(t) for o in order_factors):
383 newseq2.append(t)
384 newseq = newseq2 + order_factors # type: ignore
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
69 @wraps(func)
70 def wrapper(*args, **kwargs):
71 try:
---> 72 retval = cfunc(*args, **kwargs)
73 except TypeError as e:
74 if not e.args or not e.args[0].startswith('unhashable type:'):
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/series/order.py:453, in Order.contains(self, expr)
450 return rv
452 obj = self.func(expr, *self.args[1:])
--> 453 return self.contains(obj)
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/core/cache.py:72, in __cacheit.<locals>.func_wrapper.<locals>.wrapper(*args, **kwargs)
69 @wraps(func)
70 def wrapper(*args, **kwargs):
71 try:
---> 72 retval = cfunc(*args, **kwargs)
73 except TypeError as e:
74 if not e.args or not e.args[0].startswith('unhashable type:'):
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/series/order.py:428, in Order.contains(self, expr)
426 for s in common_symbols:
427 from sympy.series.limits import Limit
--> 428 l = Limit(ratio, s, point).doit(heuristics=False)
429 if not isinstance(l, Limit):
430 l = l != 0
File ~/checkouts/readthedocs.org/user_builds/wsu-phys-581-computation/conda/latest/lib/python3.13/site-packages/sympy/series/limits.py:372, in Limit.doit(self, **hints)
370 return S.ComplexInfinity
371 else:
--> 372 raise NotImplementedError("Not sure of sign of %s" % ex)
374 # gruntz fails on factorials but works with the gamma function
375 # If no factorial term is present, e should remain unchanged.
376 # factorial is defined to be zero for negative inputs (which
377 # differs from gamma) so only rewrite for non-negative z0.
378 if z0.is_extended_nonnegative:
NotImplementedError: Not sure of sign of 2 - q
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)))
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[4], line 99
96 def get_a(n, s=2.0):
97 return
---> 99 l = np.array(list(Levin(get_a, eps=1e-12)))
Cell In[4], line 59, in Levin.__iter__(self)
57 n = self.n
58 a_n = self.get_a(n)
---> 59 self.S_n += a_n
60 if 'u' == self.type:
61 w_n = (self.beta + n) * a_n
TypeError: unsupported operand type(s) for +=: 'int' and 'NoneType'
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.
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:])
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[5], line 29
27 ax.loglog(Ns, abs(S-S2s), ':+', label='$\zeta(2)$ Improved')
28 ps = np.arange(11)
---> 29 ax.loglog(2**ps, [abs(S-s(p, 0)) for p in ps],
30 ':.', label="Richardson extrapolation")
31 ax.legend()
32 ax.set(xlabel="$N$", ylabel="$S$")
TypeError: 'float' object is not callable
(f)
\[\begin{gather*} \sum_{n=1}^{\infty}\frac{(-1)^{n}}{n}. \end{gather*}\]