How to index a Python list in an Sympy Sum?

Viewed 458

I want to compute following sum (last line of the code) in Python using Sympy:

# Theta Functions, n = 0, 1, ... , x
theta_n = [0]*(x+1)
sig = sy.symbols('sigma^2', real=True)
j = sy.symbols('j', integer=True)
theta_n[0] = 1
theta_n[1] = 1 
for n in range(2,x+1):
   theta_n[n] = sy.Sum( sig**j * theta_n[n-2*j], (j,1,int(n/2))).doit()

but I get the following error

TypeError: list indices must be integers or slices, not Add

I am a complete Sympy newbie and would like to know how to do this in the right way. One workaround would be to compute the sum with a for loop which works nice but i think isnt the right way to do it. Maybe i need a other representation for thetan or I can somehow convert j to an integer in the sum to get access to the elements of the python list.

2 Answers

Instead of using Sum use a loop over j

for j in range(1, n//2 + 1):
    theta_n[n] += theta_n[n-2*j]*sig**j

These two lines replace your single line with Sum.

Once you look at the results that you are getting you might notice the pattern and rewrite the nth term as

>>> theta_n = lambda x: 2**max(0, x//2-1)*sig**(x//2)
>>> theta_n(100)
562949953421312*sigma^2**50

The main problem is that you try to use a sympy symbol (j) to index into a list. List indexing only works with Python integers. In your case, you can use j to a normal Python integer, and use a regular Python sum. This regular sum will get converted first to addition and then to sympy Add when needed. Note that the Python range command needs one number higher than the last index.

import sympy as sy

x = 12

theta_n = [0]*(x+1)
sig = sy.symbols('sigma^2', real=True)
theta_n[0] = 1
theta_n[1] = 1
for n in range(2, x+1):
    theta_n[n] = sum(sig ** j * theta_n[n - 2 * j] for j in range(1, n//2+1) )

Note that in the sympy expression, sy.Sum( sig**j * theta_n[n-2*j], (j,1,int(n/2))), j is a complete symbolic variable, and j doesn't really take all the values one by one such as in a Python iteration.

Consider for example:

print(sy.Sum(j*j, (j, 1, 1_000_000_000)).doit())

Sympy immediately gives the correct result (333333333833333333500000000) by converting the expression to n**3/3 + n**2/2 + n/6. This would take a long while to calculate with a regular Python summation.

Related