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

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)

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

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,