I have two lines in the code:
print(expr.subs(t,2*l))
print(simplify(expr.subs(t,2*l)))
Which yield:
2*Sum(2**(-2*j + 2*l + 1)*2**(2*j + 2*l - 3)*(l + 1)**2*factorial(-j + 2*l + 1)*factorial(j + 2*l - 1)/((-j + l + 1)*(j + l - 1)*factorial(j)*factorial(2 - j)*factorial(-2*j + 2*l + 1)*factorial(2*j + 2*l - 3)), (j, 0, 0))
0
The first command works correctly and the expression is what it should be, but somehow the output of the second command is 0. When I calculate this doubled sum manually, it is not 0. There is only one term in the sum (as j goes from 0 to 0). Can somebody explain why the last command returns 0?
You have an inner expression that looks like:
from sympy import Sum, factorial, sin , gamma, pi
from sympy.abc import j, l
inner = (2 ** (-2 * j + 2 * l + 1)
* 2 ** (2 * j + 2 * l - 3)
* (l + 1) ** 2
* factorial(-j + 2 * l + 1)
* factorial(j + 2 * l - 1)
/ ((-j + l + 1)
* (j + l - 1)
* factorial(j)
* factorial(2 - j)
* factorial(-2 * j + 2 * l + 1)
* factorial(2 * j + 2 * l - 3)))
Calling inner.simplify() substitutes factorial() with the gamma() function. This substitution has some conditions, which sometimes may make it invalid. The simplified version looks like:
(16 ** l * (l + 1) ** 2
* sin(pi * j)
* gamma(j + 2 * l)
* gamma(-j + 2 * l + 2)
/ (pi * j * (j - 2) * (j - 1)
* gamma(-2 * j + 2 * l + 3)
* gamma(2 * j + 2 * l - 1)))
Then, substituting j with 0 gives an expression of zero (due to sin(pi * j)) divided by zero (due to pi * j in the denominator). Further evaluation lets sympy conclude that the limit value must be zero.
All this can be avoided by first calling .doit() (this tries to evaluate the unevaluated sum).
expr.subs(t,2*l).doit().simplify()
# gives 2**(4*l - 1)*(l + 1)*(2*l - 1)