"partial_fraction_decomposition" with possibly "complex roots", again

edit

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[1], s, u)
sage: u = SR.var('u')
sage: inverse_laplace(b[0], s, u)
4/7*e^(-u)
sage: inverse_laplace(b[1], s, u)
2/21*e^(-10/3*u)

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()[0].radical_expression()
            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))

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,