Ask Your Question

# "partial_fraction_decomposition" with possibly "complex roots", again

Hi I come back to this question, even though it's been answered before, since I am still not able to make it work. I should mention maybe that I am just trying to teach undergraduate students to invert rational Laplace transforms (for myself I am able to afford a Mathematica licence, but it would be nice to be able to show students that such simple things may be done nowadays for free). Following a previous answer, I tried

L = 2*(s + 3)/(3*s^2 + 13*s + 10)
Ks=FractionField(PolynomialRing(CC,names='s'))
Lr=Ks(L)


Of course, with quadratic rational roots, I could do this by hand, but the purpose here is to do it when you do not know the roots. The first two commands work, but the third has error message

TypeError: ('cannot convert {!r}/{!r} to an element of {}', 2*(s + 3)/(3*s^2 + 13*s + 10), 1.00000000000000, Fraction Field of Univariate Polynomial Ring in s over Complex Field with 53 bits of precision)


I should add that a different attempt which used to work last year

C=ComplexField(53);
dec=Frac(C['s'])(Lrs).partial_fraction_decomposition();


gets now same error message. So, this is probably due to an "improvement" of sage

edit retag close merge delete

## 3 answers

Sort by » oldest newest most voted

Thanks for reporting this issue.

Doing things in a slightly different order works well:

sage: R.<s> = PolynomialRing(CC)
sage: F = R.fraction_field()
sage: L = 2*(s + 3)/(3*s^2 + 13*s + 10)
sage: L
(2.00000000000000*s + 6.00000000000000)/(3.00000000000000*s^2 + 13.0000000000000*s + 10.0000000000000)
sage: L.partial_fraction_decomposition()
(0,
[0.571428571428571/(s + 1.00000000000000),
0.0952380952380953/(s + 3.33333333333333)])


Note that in this case where all coefficients are integer, one could work with exact algebraic numbers:

sage: R.<s> = PolynomialRing(QQbar)
sage: F = R.fraction_field()
sage: L = 2*(s + 3)/(3*s^2 + 13*s + 10)
sage: L
(2*s + 6)/(3*s^2 + 13*s + 10)
sage: L.partial_fraction_decomposition()
(0, [0.571428571428572?/(s + 1), 2/21/(s + 10/3)])


Edit (to answer follow-up question):

In this case one can even work over the rationals.

sage: R.<s> = PolynomialRing(QQ)
sage: F = R.fraction_field()
sage: L = 2*(s + 3)/(3*s^2 + 13*s + 10)
sage: L
(2*s + 6)/(3*s^2 + 13*s + 10)
sage: f = L.partial_fraction_decomposition()
sage: f
(0, [4/7/(s + 1), 2/21/(s + 10/3)])


Then you can take Laplace transform of each partial fraction:

sage: a, b = f
sage: a
0
sage: b
[4/7/(s + 1), 2/21/(s + 10/3)]
sage: inverse_laplace(b, s, u)
sage: u = SR.var('u')
sage: inverse_laplace(b, s, u)
4/7*e^(-u)
sage: inverse_laplace(b, s, u)
2/21*e^(-10/3*u)

more

## Comments

Thanks, Samuel :) Starting with the declaration works better :) It's still a bit unfortunate though that the first coefficient 4/7 is not obtained as fraction. I need then an inverse_laplace and this seems to work only for second component

var('s,u')
R.<s> = PolynomialRing(QQbar)
F = R.fraction_field()
L = 2*(s + 3)/(3*s^2 + 13*s + 10)
whole,LL=L.partial_fraction_decomposition()
show(LL)
inverse_laplace(LL,s,u)


It would be nice if there was a simple way to convert all the ?numbers in my list of fractions to closest fractions (apply nearby_rational? to all parts of an expression)

Why settle for an approximation (nearby_rational) when an exact solution (radical_expression) can be attempted ?

@florin -- in this case you can just work over QQ instead of QQbar. I edited my answer.

Samuel and Emanuel, thanks for teaching me the magic of QQbar and SR! one more question, please. How to do this inverse_laplace of a fraction whose denominator roots may be real (and maybe even complex)? For example,

var('s,u')
R.<s> = PolynomialRing(QQbar)
F = R.fraction_field()
L=3/4*(19*s^2 + 156*s + 284)/(19*s^3 + 174*s^2 + 422*s + 228)
whole,LL=L.partial_fraction_decomposition()
show(LL)


Here I'm stuck with the "decimals?" numbers and

inverse_laplace(L,s,u)


ends again in

TypeError: unable to convert [0.6934316866155686?/(s + 0.7570751092241817?), 0.0477251478733486?/(s + 2.861392366086696?), 0.00884316551108293?/(s + 5.539427261531228?)] to a symbolic expression


on the other hand

SR(L).partial_fraction()


fails to find the real roots, and inverse_laplace fails also. Thanks, F

The problem is really how to do partial_fractions and then extract the answer, but I added more of the story to make it clear this is a numeric question. For a first answer, I do not need 100 digits precision :)

I have added a new question "inverse_laplace of a fraction whose denominator roots are real (or complex) to relance this with more details

more

## Comments

This is probably due to the fact that the roots of the denominator a real, but can't be expressed as a real radical expression (see Wikipedia's casus irreductibilis). This deserves further discussion : would you mind asking a new question ? (this will help future users of ask.sagemath.org with a similar problem...).

Well, there are cases where brute force helps :

var('s,u')
R.<s> = PolynomialRing(QQbar)
F = R.fraction_field()
L = 2*(s + 3)/(3*s^2 + 13*s + 10)
whole,LL=L.partial_fraction_decomposition()
print LL
[0.571428571428572?/(s + 1), 2/21/(s + 10/3)]

def get_rad_exp(l):
"""
Returns a sum of polynomial fractions expressed as a sum of
fractions of radical expressions, if possible.
"""
def gre_f(f):
"""
Returns a polynomial fraction expressed as a ration of
radical expressions if possible.
"""
def gre_p(p):
"""
Returns a polynomial whose coefficients are expressed as
radical expressions, if possible
"""
v=p.variables()
## print "p = ", p, ", v = ", v
if len(v)==0: ## Constant
return p.coefficients().radical_expression()
if (len(v)>1): ## Multi-variable case : TBD
raise ValueError("Multi-variable case not (yet) handled")
v=v
L=p.coefficients(sparse=False)
D=range(p.degree()+1)
return sum(map(lambda c,d:c.radical_expression()*v^d, L, D))
return gre_p(f.numerator())/gre_p(f.denominator())
return sum(map(gre_f, l))

print get_rad_exp(LL)
2/21/(s + 10/3) + 4/7/(s + 1)


Which is, unless I'm mistaken, the desired results...

And, BTW,

print mathematica.Apart(L).sage()
2/7/(3*s + 10) + 4/7/(s + 1)


Finally :

print inverse_laplace(get_rad_exp(LL),s,u)
4/7*e^(-u) + 2/21*e^(-10/3*u)


should be what you aimed to do...

But things do not need to be so complicated :

L.partial_fraction_decomposition()
(0, [0.571428571428572?/(s + 1), 2/21/(s + 10/3)])
SR(L).partial_fraction()
2/7/(3*s + 10) + 4/7/(s + 1)
sage: inverse_laplace(SR(L).partial_fraction(),s,u)
4/7*e^(-u) + 2/21*e^(-10/3*u)


Even simpler :

Lt=SR(repr(L).replace("s","t"));Lt
2*(t + 3)/(3*t^2 + 13*t + 10)
inverse_laplace(Lt,t,u)
4/7*e^(-u) + 2/21*e^(-10/3*u)


HTH,

more

## Your Answer

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

Add Answer

## Stats

Asked: 2018-04-30 03:58:31 -0500

Seen: 112 times

Last updated: May 14 '18