ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Mon, 14 May 2018 03:24:36 -0500"partial_fraction_decomposition" with possibly "complex roots", againhttp://ask.sagemath.org/question/42220/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 sageMon, 30 Apr 2018 03:58:31 -0500http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/Answer by florin for <p>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</p>
<pre><code>L = 2*(s + 3)/(3*s^2 + 13*s + 10)
Ks=FractionField(PolynomialRing(CC,names='s'))
Lr=Ks(L)
</code></pre>
<p>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</p>
<pre><code>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)
</code></pre>
<p>I should add that a different attempt which used to work last year </p>
<pre><code>C=ComplexField(53);
dec=Frac(C['s'])(Lrs).partial_fraction_decomposition();
</code></pre>
<p>gets now same error message. So, this is probably due to an "improvement" of sage</p>
http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?answer=42327#post-id-42327Samuel 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 detailsSun, 13 May 2018 14:04:35 -0500http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?answer=42327#post-id-42327Comment by Emmanuel Charpentier for <p>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,</p>
<pre><code>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)
</code></pre>
<p>Here I'm stuck with the "decimals?" numbers and </p>
<pre><code>inverse_laplace(L,s,u)
</code></pre>
<p>ends again in</p>
<pre><code>TypeError: unable to convert [0.6934316866155686?/(s + 0.7570751092241817?), 0.0477251478733486?/(s + 2.861392366086696?), 0.00884316551108293?/(s + 5.539427261531228?)] to a symbolic expression
</code></pre>
<p>on the other hand</p>
<pre><code>SR(L).partial_fraction()
</code></pre>
<p>fails to find the real roots, and inverse_laplace fails also. Thanks, F</p>
<p>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 :)</p>
<p>I have added a new question "inverse_laplace of a fraction whose denominator roots are real (or complex) to relance this
with more details</p>
http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?comment=42329#post-id-42329This 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](https://en.wikipedia.org/wiki/Casus_irreducibilis)). This deserves further discussion : would you mind asking a new question ? (this will help future users of `ask.sagemath.org` with a similar problem...).Mon, 14 May 2018 03:24:36 -0500http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?comment=42329#post-id-42329Answer by slelievre for <p>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</p>
<pre><code>L = 2*(s + 3)/(3*s^2 + 13*s + 10)
Ks=FractionField(PolynomialRing(CC,names='s'))
Lr=Ks(L)
</code></pre>
<p>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</p>
<pre><code>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)
</code></pre>
<p>I should add that a different attempt which used to work last year </p>
<pre><code>C=ComplexField(53);
dec=Frac(C['s'])(Lrs).partial_fraction_decomposition();
</code></pre>
<p>gets now same error message. So, this is probably due to an "improvement" of sage</p>
http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?answer=42221#post-id-42221Thanks 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)Mon, 30 Apr 2018 04:34:42 -0500http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?answer=42221#post-id-42221Comment by slelievre for <p>Thanks for reporting this issue.</p>
<p>Doing things in a slightly different order works well:</p>
<pre><code>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)])
</code></pre>
<p>Note that in this case where all coefficients are integer,
one could work with exact algebraic numbers:</p>
<pre><code>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)])
</code></pre>
<p>Edit (to answer follow-up question):</p>
<p>In this case one can even work over the rationals.</p>
<pre><code>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)])
</code></pre>
<p>Then you can take Laplace transform of each partial fraction:</p>
<pre><code>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)
</code></pre>
http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?comment=42270#post-id-42270@florin -- in this case you can just work over QQ instead of QQbar. I edited my answer.Sun, 06 May 2018 13:26:23 -0500http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?comment=42270#post-id-42270Comment by florin for <p>Thanks for reporting this issue.</p>
<p>Doing things in a slightly different order works well:</p>
<pre><code>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)])
</code></pre>
<p>Note that in this case where all coefficients are integer,
one could work with exact algebraic numbers:</p>
<pre><code>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)])
</code></pre>
<p>Edit (to answer follow-up question):</p>
<p>In this case one can even work over the rationals.</p>
<pre><code>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)])
</code></pre>
<p>Then you can take Laplace transform of each partial fraction:</p>
<pre><code>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)
</code></pre>
http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?comment=42265#post-id-42265Thanks, 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)
Sun, 06 May 2018 04:46:25 -0500http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?comment=42265#post-id-42265Comment by Emmanuel Charpentier for <p>Thanks for reporting this issue.</p>
<p>Doing things in a slightly different order works well:</p>
<pre><code>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)])
</code></pre>
<p>Note that in this case where all coefficients are integer,
one could work with exact algebraic numbers:</p>
<pre><code>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)])
</code></pre>
<p>Edit (to answer follow-up question):</p>
<p>In this case one can even work over the rationals.</p>
<pre><code>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)])
</code></pre>
<p>Then you can take Laplace transform of each partial fraction:</p>
<pre><code>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)
</code></pre>
http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?comment=42268#post-id-42268Why settle for an approximation (`nearby_rational`) when an exact solution (`radical_expression`) can be attempted ?Sun, 06 May 2018 12:02:28 -0500http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?comment=42268#post-id-42268Answer by Emmanuel Charpentier for <p>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</p>
<pre><code>L = 2*(s + 3)/(3*s^2 + 13*s + 10)
Ks=FractionField(PolynomialRing(CC,names='s'))
Lr=Ks(L)
</code></pre>
<p>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</p>
<pre><code>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)
</code></pre>
<p>I should add that a different attempt which used to work last year </p>
<pre><code>C=ComplexField(53);
dec=Frac(C['s'])(Lrs).partial_fraction_decomposition();
</code></pre>
<p>gets now same error message. So, this is probably due to an "improvement" of sage</p>
http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?answer=42267#post-id-42267Well, 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,
Sun, 06 May 2018 11:44:14 -0500http://ask.sagemath.org/question/42220/partial_fraction_decomposition-with-possibly-complex-roots-again/?answer=42267#post-id-42267