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