CODEWITHSUNDEEP
pythonmath
Simd
Simd

Reputation: 21223

How to perform infinite sum numerically in python

I am moving from Maple to python for my mathematical programming. As part of this I am trying to work out what the right tools are to perform infinite sums numerically.

I would like to compute numerically for example:

sum(exp(-x^2), x = -infinity..infinity)

In Maple this would just be

evalf(sum(exp(-x^2), x = -infinity..infinity));

                             1.772637205

Can you do something similar in python, perhaps using the scientific python ecosystem (scipy,numpy, sympy etc)?

Upvotes: 2

Views: 22473

Answers (3)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

For terms that fall fast enough to zero

def infinisum(f):
    n, res = 0, f(0)
    while True:
        term = sum( f(k) for k in range(2**n,2**(n+1)) )
        if (res+term)-res == 0:
             break;
        n,res = n+1, res+term
    return res

Then your particular series is

2*infinisum(lambda x: exp(-x**2)) - 1

giving 1.772637204826652 in 3 steps, or in general, forgetting that this function is symmetric,

f = lambda x:exp(-x**2)
infinisum(f)+infinisum(lambda x : f(-1-x))

Upvotes: 2

DSM
DSM

Reputation: 353059

I've used mpmath's nsum for this in the past with good results:

>>> from mpmath import nsum, exp, inf
>>> nsum(lambda x: exp(-x**2), [-inf, inf])
mpf('1.7726372048266521')

which is slightly different from

>>> from mpmath import quad
>>> quad(lambda x: exp(-x**2), [-inf, inf])
mpf('1.7724538509055161')

We can get it at higher precision, and compare it against the analytic value:

>>> import mpmath
>>> mpmath.mp.dps = 50
>>> nsum(lambda x: exp(-x**2), [-inf, inf])
mpf('1.7726372048266521530312505511578584813433860453722459')
>>> mpmath.jtheta(3, 0, 1/exp(1))
mpf('1.7726372048266521530312505511578584813433860453722465')

Upvotes: 6

blazs
blazs

Reputation: 4845

I expect you could do some of the sums by summing the range where the value of the expression you are summing is sufficiently large (i.e. bigger than some prespecified eps). But I don't know how one could solve this more generally.

Whenever possible you should of course use the closed-form. For example, 1+x+x^2+x^3+... (up to infinity) sums to 1/(1-x) when |x|<1.

Upvotes: 0

Related Questions