sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.8k stars 4.39k forks source link

summation over poles of rational functions #8746

Open akritas opened 9 years ago

akritas commented 9 years ago

summation( 1 / (k * (k + 1)), (k, -3, n) ) returns -(n + 4)/(3*(n + 1)) whereas it should return nothing as Wolframalpha and Maxima.

gxyd commented 9 years ago

I have modified the 'summation' method in the Concrete part. It seems to be working fine for me.. Can anyone tell me about it... ?

def summation(f, _symbols, *_kwargs): r""" Compute the summation of f with respect to symbols.

The notation for symbols is similar to the notation used in Integral.
summation(f, (i, a, b)) computes the sum of f with respect to i from a to b,
i.e.,

::

                                b
                              ____
                              \   `
    summation(f, (i, a, b)) =  )    f
                              /___,
                              i = a

If it cannot compute the sum, it returns an unevaluated Sum object.
Repeated sums can be computed by introducing additional symbols tuples::

>>> from sympy import summation, oo, symbols, log
>>> i, n, m = symbols('i n m', integer=True)

>>> summation(2*i - 1, (i, 1, n))
n**2
>>> summation(1/2**i, (i, 0, oo))
2
>>> summation(1/log(n)**n, (n, 2, oo))
Sum(log(n)**(-n), (n, 2, oo))
>>> summation(i, (i, 0, n), (n, 0, m))
m**3/6 + m**2/2 + m/3

>>> from sympy.abc import x
>>> from sympy import factorial
>>> summation(x**n/factorial(n), (n, 0, oo))
exp(x)

See Also
========

Sum
Product, product

"""

''' def L(expr):
e1 = simplify(expr)
denom = e1.as_numer_denom()[1]
r = solve(denom)
r_int = []
for i in r:
    if isinstance(i, Integer)==True:
        r_int.append(i)
return r_int'''

e1 = simplify(f)
denom = e1.as_numer_denom()[1]
r = solve(denom)
r_int = []
for i in r:
    if isinstance(i, Integer):
        r_int.append(i)

if r_int!=[]:
    if isinstance(symbols[0][1], int) and isinstance(symbols[0][2], Symbol):
        for val in r_int:
            if val>symbols[0][1]:
                return "Sum(%s, %s)" %(f, symbols[0])
        for val in r_int:
            if val==symbols[0][1]:
                return "Sum(%s, %s)" %(f, symbols[0])
    elif isinstance(symbols[0][1], Symbol) and isinstance(symbols[0][2], int):
        gr = 0
        lt = 0
        for val in r_int:
            if val>symbols[0][2]:
                gr=gr+1
            elif val<symbols[0][2]:
                lt=lt+1
            if gr!=0 and lt!=0:
                return "Sum(%s, %s)" %(f, symbols[0])
        if gr==0 or lt==0:
            for val in r_int:
                if val==symbols[0][2]:
                    return "Sum(%s, %s)" %(f, symbols[0])

return Sum(f, *symbols, **kwargs).doit(deep=False)

///

The tests are as:

print(summation(1/(k(k+3)), (k, -5, n))) Sum(1/(k(k + 3)), (k, -5, n))

print(summation(1/(k(k+3)), (k, -3, n))) Sum(1/(k(k + 3)), (k, -3, n))

print(summation(1/(k(k+3)), (k, -2, n))) Sum(1/(k(k + 3)), (k, -2, n))

print(summation(1/(k(k+3)), (k, 0, n))) Sum(1/(k(k + 3)), (k, 0, n))

print(summation(1/(k(k+3)), (k, 2, n))) (n - 1)(13_n__2 + 55n + 54)/(36(n + 1)_(n + 2)*(n + 3))

print(summation(1/(k(k+3)), (k, n, -5))) (n + 4)(13_n__2 + 23_n + 6)/(36n(n + 1)*(n + 2))

print(summation(1/(k(k+3)), (k, n, -2))) Sum(1/(k(k + 3)), (k, n, -2))

print(summation(1/(k(k+3)), (k, n, 0))) Sum(1/(k(k + 3)), (k, n, 0))

print(summation(1/(k(k+3)), (k, n, 2))) -(n - 3)(47_n__2 + 102_n + 40)/(180n(n + 1)*(n + 2))

skirpichev commented 9 years ago

@gxyd, please read this

asmeurer commented 9 years ago

@gxyd it would be easier to review the code if you open a pull request, even if the code is not ready to be merged yet.

asmeurer commented 9 years ago

Wouldn't a better answer be a Piecewise (c.f. what integrate() does in similar situations). In the example above, if n <= -2 then the summation does not go across any poles.

gxyd commented 9 years ago

@asmeurer the pull request is at #8839 . i think it will be difficult to represent the summation in that case. For ex. summation(1/(k*(k+1)), (k, 2, n)) . If n== pi, in this case the documentation says should include all the integer value from 2 to pi. for that 'greatest integer function should be used'.

asmeurer commented 9 years ago

I'm not sure what our policy is regarding summation where the bounds are noninteger, but note that summations can have reversed bounds (using the Karr convention, see the docstring of Sum).

asmeurer commented 9 years ago

@raoulb thoughts on this?

gxyd commented 9 years ago

@asmeurer According to the current definition of the summation in the docstring, the issue i am fixing does not impact the answer.