sympy/series/acceleration.py

   1 """
   2 Convergence acceleration / extrapolation methods for series and
   3 sequences.
   4
   5 References:
   6 Carl M. Bender & Steven A. Orszag, "Advanced Mathematical Methods for
   7 Scientists and Engineers: Asymptotic Methods and Perturbation Theory",
   8 Springer 1999. (Shanks transformation: pp. 368-375, Richardson
   9 extrapolation: pp. 375-377.)
  10 """
  11
  12 from sympy import Symbol, factorial, pi, Integer, E, Sum2
  13
  14
  15 def richardson(A, k, n, N):
  16     """
  17     Calculate an approximation for lim k->oo A(k) using Richardson
  18     extrapolation with the terms A(n), A(n+1), ..., A(n+N+1).
  19     Choosing N ~= 2*n often gives good results.
  20
  21     A simple example is to calculate exp(1) using the limit definition.
  22     This limit converges slowly; n = 100 only produces two accurate
  23     digits:
  24
  25         >>> n = Symbol('n')
  26         >>> e = (1 + 1/n)**n
  27         >>> print round(e.subs(n, 100).evalf(), 10)
  28         2.7048138294
  29
  30     Richardson extrapolation with 11 appropriately chosen terms gives
  31     a value that is accurate to the indicated precision:
  32
  33         >>> print round(richardson(e, n, 10, 20).evalf(), 10)
  34         2.7182818285
  35         >>> print round(E.evalf(), 10)
  36         2.7182818285
  37
  38     Another useful application is to speed up convergence of series.
  39     Computing 100 terms of the zeta(2) series 1/k**2 yields only
  40     two accurate digits:
  41
  42         >>> k = Symbol('k'); n = Symbol('n')
  43         >>> A = Sum2(k**-2, (k, 1, n))
  44         >>> print round(A.subs(n, 100).evalf(), 10)
  45         1.6349839002
  46
  47     Richardson extrapolation performs much better:
  48
  49         >>> print round(richardson(A, n, 10, 20).evalf(), 10)
  50         1.6449340668
  51         >>> print round(((pi**2)/6).evalf(), 10)     # Exact value
  52         1.6449340668
  53
  54     """
  55     s = Integer(0)
  56     for j in range(0, N+1):
  57         s += A.subs(k, Integer(n+j)) * (n+j)**N * (-1)**(j+N) / \
  58             (factorial(j) * factorial(N-j))
  59     return s
  60
  61
  62 def shanks(A, k, n, m=1):
  63     """
  64     Calculate an approximation for lim k->oo A(k) using the n-term Shanks
  65     transformation S(A)(n). With m > 1, calculate the m-fold recursive
  66     Shanks transformation S(S(...S(A)...))(n).
  67
  68     The Shanks transformation is useful for summing Taylor series that
  69     converge slowly near a pole or singularity, e.g. for log(2):
  70
  71         >>> n = Symbol('n')
  72         >>> k = Symbol('k')
  73         >>> A = Sum2(Integer(-1)**(k+1) / k, (k, 1, n))
  74         >>> print round(A.subs(n, 100).evalf(), 10)
  75         0.6881721793
  76         >>> print round(shanks(A, n, 25).evalf(), 10)
  77         0.6931396564
  78         >>> print round(shanks(A, n, 25, 5).evalf(), 10)
  79         0.6931471806
  80
  81     The correct value is 0.6931471805599453094172321215.
  82     """
  83     table = [A.subs(k, Integer(j)) for j in range(n+m+2)]
  84     table2 = table[:]
  85
  86     for i in range(1, m+1):
  87         for j in range(i, n+m+1):
  88             x, y, z = table[j-1], table[j], table[j+1]
  89             table2[j] = (z*x - y**2) / (z + x - 2*y)
  90         table = table2[:]
  91     return table[n]
  92