Ask Your Question
2

"partial_fraction_decomposition" with possibly "complex roots", again

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

florin gravatar image

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 flag offensive close merge delete

3 answers

Sort by ยป oldest newest most voted
1

answered 2018-04-30 04:34:42 -0500

slelievre gravatar image

updated 2018-05-06 13:36:00 -0500

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)
edit flag offensive delete link 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[1],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)

florin gravatar imageflorin ( 2018-05-06 04:46:25 -0500 )edit

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

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2018-05-06 12:02:28 -0500 )edit

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

slelievre gravatar imageslelievre ( 2018-05-06 13:26:23 -0500 )edit
0

answered 2018-05-13 14:04:35 -0500

florin gravatar image

updated 2018-05-14 04:13:58 -0500

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

edit flag offensive delete link 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...).

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2018-05-14 03:24:36 -0500 )edit
0

answered 2018-05-06 11:44:14 -0500

Emmanuel Charpentier gravatar image

updated 2018-05-06 12:00:54 -0500

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,

edit flag offensive delete link more

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

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

Seen: 78 times

Last updated: May 14