3

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)?

Community
  • 1
  • 1
Simd
  • 19,447
  • 42
  • 136
  • 271
  • 1
    What do you expect to happen if the sum doesn't converge, or if its convergence is [unknown](http://math.stackexchange.com/questions/20555/are-there-any-series-whose-convergence-is-unknown)? (I'm curious how Maple handles this.) In general Python doesn't handle operations like this that involve symbolic/theoretical constructs like summing "to infinity". You can compute the sum between whatever finite limits you want, but you have to set the finite limits yourself (and decide for yourself if you consider the result to converge). – BrenBarn Apr 03 '16 at 07:01
  • @BrenBarn If it doesn't converge or is generally very hard to compute I don't expect it to succeed. – Simd Apr 03 '16 at 07:12
  • Nothing infinite can be done on a computer in a finite period of time. You can calculate by adding terms until a desired level of accuracy is achieved. You also need to know that IEEE double precision floating point numbers only have a limited accuracy: 17 significant digits. You write "sum". Do you really mean "integral"? – duffymo Sep 25 '19 at 17:08

3 Answers3

7

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')
DSM
  • 342,061
  • 65
  • 592
  • 494
  • 1
    I had no idea nsum existed! Thank you. I believe 1.77263720482665215303125055115785848134338604537224605 is the correct answer up to 54 significant figures. – Simd Apr 03 '16 at 16:01
  • Using Mathematica, you will get the result, `EllipticTheta(3,0,1/e)`. – Jens Munk Apr 03 '16 at 19:28
2

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))
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
0

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.

blazs
  • 4,705
  • 24
  • 38