Problem with partial_fraction

partial_fraction

asked 2014-12-16 11:27:21 +0100

Peter Luschny
548 ●14 ●37 ●58

updated 2014-12-17 11:38:52 +0100

With Maple I can compute the partial fractions of the sums given below

M := [1, x, x^2 + 1, x^3 + 3*x, x^4 + 6*x^2 + 2, x^5 + 10*x^3 + 10*x, x^6 + 15*x^4 + 30*x^2 + 5]:
for k from 1 to 7 do
p := expand(sum(x^j*subs(x=j,M[k]),j=0..infinity)):
lprint([k], convert(p,parfrac)) od:

and get:

[1], -1/(x-1)
[2], 1/(x-1)+1/(x-1)^2
[3], -2/(x-1)-3/(x-1)^2-2/(x-1)^3
[4], 6/(x-1)^4+4/(x-1)+10/(x-1)^2+12/(x-1)^3
[5], -60/(x-1)^4-9/(x-1)-24/(x-1)^5-33/(x-1)^2-62/(x-1)^3
[6], 120/(x-1)^6+450/(x-1)^4+21/(x-1)+360/(x-1)^5+111/(x-1)^2+300/(x-1)^3
[7], -2520/(x-1)^6-3000/(x-1)^4-51/(x-1)-720/(x-1)^7-3720/(x-1)^5-378/(x-1)^2-1412/(x-1)^3

With Sage I can do

from sage.calculus.calculus import symbolic_sum
x, j = SR.var('x, j')
M = [1, x, x^2 + 1, x^3 + 3*x, x^4 + 6*x^2 + 2, x^5 + 10*x^3 + 10*x, x^6 + 15*x^4 + 30*x^2 + 5]
for k in (0..6):
    p = symbolic_sum(x^j*M[k](x=j), j, 0, oo)
    print [k], p.partial_fraction()

and get:

[0] TypeError: 'sage.rings.integer.Integer' object is not callable
[1] 1/(x - 1) + 1/(x - 1)^2
[2] -2/(x - 1) - 3/(x - 1)^2 - 2/(x - 1)^3
[3] 4/(x - 1) + 10/(x - 1)^2 + 12/(x - 1)^3 + 6/(x - 1)^4
[4] 2*hypergeometric((1, -1/2*sqrt(-4*sqrt(7) - 12) + 1, 1/2*sqrt(-4*sqrt(7) - 12) + 1, 
    -1/2*sqrt(4*sqrt(7) - 12) + 1, 1/2*sqrt(4*sqrt(7) - 12) + 1), (-sqrt(-sqrt(7) - 3), 
    sqrt(-sqrt(7) - 3), -sqrt(sqrt(7) - 3), sqrt(sqrt(7) - 3)), x)
[5] 21/(x - 1) + 111/(x - 1)^2 + 300/(x - 1)^3 + 450/(x - 1)^4 + 360/(x - 1)^5 + 120/(x - 1)^6
[6] TypeError: ECL says: Error executing code in Maxima: expt: undefined: 0 to a negative exponent.

Is it true that something basic like the partial fraction decomposition is completely unusable in Sage?

Edit: By request of rws the input p explicit:

[0], -1/(x-1)
[1], x/(x-1)^2
[2], -(2*x^2-x+1)/(x-1)^3
[3], 2*x*(-x+2*x^2+2)/(x-1)^4
[4], -(2-x+17*x^2-3*x^3+9*x^4)/(x-1)^5
[5], 3*x*(7+2*x+22*x^2+2*x^3+7*x^4)/(x-1)^6
[6], -(5+16*x+177*x^2+112*x^3+287*x^4+72*x^5+51*x^6)/(x-1)^7

edit retag flag offensive close merge delete

Comments

Can you provide the input for partial_fraction for 4 and 6?

rws ( 2014-12-16 17:07:16 +0100 )edit

For reference: other ask-sage questions about infinite sums
http://ask.sagemath.org/questions/sco... .

slelievre ( 2014-12-16 18:57:18 +0100 )edit

@rws It's there, input is M[k] for k =0..6. For example input 4 = M[4] = x^4 + 6*x^2 + 2.

Peter Luschny ( 2014-12-17 11:27:47 +0100 )edit

@Peter Luschny: I think @rws meant "can you provide p in each case, before applying .partial_fraction()". My long comment disguised as an answer does that.

slelievre ( 2014-12-17 19:20:50 +0100 )edit

@slelievre This occurred to me right after I had written my comment. See my edit to the question above, which I wrote immediately after the comment and you apparently overlooked. It is easy to get: take the solutions of Maple and press 'simplify'. At least this gives you the correct input.

Peter Luschny ( 2014-12-18 17:24:47 +0100 )edit

add a comment

2 Answers

Sort by » oldest newest most voted

answered 2014-12-19 08:27:07 +0100

rws
1149 ●2 ●20 ●44 http://google.com/+Ral...

updated 2014-12-19 08:30:19 +0100

It's not partial_fraction that does the unexpected here but symbolic_sum, and actually it seems correct. Granted, it's not obvious, but you must assume abs(x)<1 beforehand, for convergence reasons:

  sage: assume(abs(x)<1)

With that, I get from your code:

sage: for k in (1..6):
        p = symbolic_sum(x^j*M[k](x=j), j, 0, oo)
        print [k], p.partial_fraction()
....:     
[1] 1/(x - 1) + 1/(x - 1)^2
[2] -2/(x - 1) - 3/(x - 1)^2 - 2/(x - 1)^3
[3] 4/(x - 1) + 10/(x - 1)^2 + 12/(x - 1)^3 + 6/(x - 1)^4
[4] -9/(x - 1) - 33/(x - 1)^2 - 62/(x - 1)^3 - 60/(x - 1)^4 - 24/(x - 1)^5
[5] 21/(x - 1) + 111/(x - 1)^2 + 300/(x - 1)^3 + 450/(x - 1)^4 + 360/(x - 1)^5 + 120/(x - 1)^6
[6] -51/(x - 1) - 378/(x - 1)^2 - 1412/(x - 1)^3 - 3000/(x - 1)^4 - 3720/(x - 1)^5 - 2520/(x - 1)^6 - 720/(x - 1)^7

edit flag offensive delete link

Comments

Well done! This is also quite fast. In addition, it also works for k = 0, if the first entry of M is SR(1).

slelievre ( 2014-12-19 10:18:59 +0100 )edit

Very nice! But then: a symbolic sum does, by definition, make no assumptions on convergence! So this is a serious bug nevertheless. By the way, motivation comes from the Motzkin polynomials, see A247497.

Peter Luschny ( 2014-12-20 08:37:36 +0100 )edit

I don't think so. It's just a power series, and the only criterion for being a '''formal''' power series is the assumption given. However, to have something like formal_sum(...) would be worthwhile in Sage, and this could be then made very fast for polynomials.

rws ( 2014-12-20 16:00:03 +0100 )edit

add a comment

answered 2014-12-16 18:54:36 +0100

slelievre
17704 ●22 ●160 ●350 http://carva.org/samue...

Not able to answer, I am adding some comments which would not fit in a comment.

The TypeError exception for k = 0 comes from trying to evaluate M[k](x=j), which becomes 1(x=j), but the integer 1 is not callable.

One idea to fix this is to make the first entry in your list M an element of the symbolic ring.

sage: from sage.calculus.calculus import symbolic_sum
sage: x, j = SR.var('x, j')
sage: M = [SR(1), x, x^2 + 1, x^3 + 3*x, x^4 + 6*x^2 + 2, x^5 + 10*x^3 + 10*x, x^6 + 15*x^4 + 30*x^2 + 5]

But this leads to another failure, with a ValueError exception. I don't know how to work around that.

For k = 1 to 6, here are more step by step results. I can't explain the failures. Note the suprising output for k = 4 and the atrocious one for k = 6.

sage: def show_steps(k):
....:     p = symbolic_sum(x^j*M[k](x=j), j, 0, oo)
....:     print '\n-----   k =', k, '  -----'
....:     print '\np =', p
....:     print '\nf =', p.partial_fraction()

sage: for k in (1..6): show_steps(k)

-----   k = 1   -----

p = x/(x^2 - 2*x + 1)

f = 1/(x - 1) + 1/(x - 1)^2

-----   k = 2   -----

p = -(2*x^2 - x + 1)/(x^3 - 3*x^2 + 3*x - 1)

f = -2/(x - 1) - 3/(x - 1)^2 - 2/(x - 1)^3

-----   k = 3   -----

p = 2*(2*x^3 - x^2 + 2*x)/(x^4 - 4*x^3 + 6*x^2 - 4*x + 1)

f = 4/(x - 1) + 10/(x - 1)^2 + 12/(x - 1)^3 + 6/(x - 1)^4

-----   k = 4   -----

p = 2*hypergeometric((1, -1/2*sqrt(4*sqrt(7) - 12) + 1, 1/2*sqrt(4*sqrt(7) - 12) + 1, -1/2*sqrt(-4*sqrt(7) - 12) + 1, 1/2*sqrt(-4*sqrt(7) - 12) + 1), (-sqrt(-sqrt(7) - 3), sqrt(-sqrt(7) - 3), -sqrt(sqrt(7) - 3), sqrt(sqrt(7) - 3)), x)

f = 2*hypergeometric((1, -1/2*sqrt(-4*sqrt(7) - 12) + 1, 1/2*sqrt(-4*sqrt(7) - 12) + 1, -1/2*sqrt(4*sqrt(7) - 12) + 1, 1/2*sqrt(4*sqrt(7) - 12) + 1), (-sqrt(-sqrt(7) - 3), sqrt(-sqrt(7) - 3), -sqrt(sqrt(7) - 3), sqrt(sqrt(7) - 3)), x)

-----   k = 5   -----

p = 3*(7*x^5 + 2*x^4 + 22*x^3 + 2*x^2 + 7*x)/(x^6 - 6*x^5 + 15*x^4 - 20*x^3 + 15*x^2 - 6*x + 1)

f = 21/(x - 1) + 111/(x - 1)^2 + 300/(x - 1)^3 + 450/(x - 1)^4 + 360/(x - 1)^5 + 120/(x - 1)^6

-----   k = 6   -----

p = -((sqrt(1/2)*(((15/2*I*sqrt(11) - 105/2)^(8/3)*(-700*I*sqrt(11) + 1900) - (175*(11*sqrt(11))*sqrt(3) - 875*I*sqrt(11) + 5225*I*sqrt(3) + 2375)*(15/2*I*sqrt(11) - 105/2)^(7/3) + (15750*I*sqrt(11) - 42750)*(15/2*I*sqrt ...

(more)

edit flag offensive delete link

Comments

So it's rather a symbolic_sum problem, i.e., Maxima. SymPy is no help:

sage: sympy.Sum(x^j*M[4](x=j), (j, 0, oo))
Sum(x**j*(j**4 + 6*j**2 + 2), (j, 0, oo))
sage: sympy.Sum(x^j*M[4](x=j), (j, 0, oo)).doit()
6*Piecewise((x*(-x - 1)/(x - 1)**3, Abs(x) < 1), (Sum(j**2*x**j, (j, 0, oo)), True)) + Piecewise((x*(-x**3 - 11*x**2 - 11*x - 1)/((-x + 1)**2*(x**3 - 3*x**2 + 3*x - 1)), Abs(x) < 1), (Sum(j**4*x**j, (j, 0, oo)), True)) + 2*Piecewise((1/(-x + 1), Abs(x) < 1), (Sum(x**j, (j, 0, oo)), True))

rws ( 2014-12-19 07:51:56 +0100 )edit

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Problem with partial_fraction

Comments

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Problem with partial_fraction edit

Comments

2 Answers

Comments

Comments

Your Answer

Question Tools

Stats

Related questions

Problem with partial_fraction