# Revision history [back]

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)]

"""
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
if (len(v)>1): ## Multi-variable case : TBD
raise ValueError("Multi-variable case not (yet) handled")
v=v[0]
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))

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)


HTH,

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)]

"""
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
if (len(v)>1): ## Multi-variable case : TBD
raise ValueError("Multi-variable case not (yet) handled")
v=v[0]
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))

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,