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.Wed, 16 May 2018 10:49:12 -0500Is it possible to change Polynomial Ring in the middle of a computation?http://ask.sagemath.org/question/42341/is-it-possible-to-change-polynomial-ring-in-the-middle-of-a-computation/ Hi, I'm trying to invert the "Pollaczek-Khinchine" Laplace transform when it is rational
This works for me at degree 2:
var('x,s')
Fx = (1/6*exp(-2*x)+5/6*exp(-6*x));rho=2/3
print('Hyperexponential claims:',Fx)
R.<s> = PolynomialRing(QQbar)#when all coefficients are not integer, use CC
FF = R.fraction_field()
L_F=laplace(Fx,x,s)#Laplace transform of F
#Compute Pollackek-Khinchine (PK) formula L_rui for the Laplace transform (LT) of ruin probability
m1=L_F(s=0)
fe=L_F/m1
Fe=(1-fe)/s
L_rui=rho*Fe/(1-rho*fe)
show(L_rui.simplify_full())
inverse_laplace(SR(L_rui),s,u)
but not at degree 3, since I do not know how to use partial_fraction_decomposition, and then to switch to RR numbers and then invert .
If I start in R.<s> = PolynomialRing(RR), for an already known LT, everything is fine. But, a certain simplification by s in PK formula will become impossible due to rounding errors, so I am forced to start with R.<s> = PolynomialRing(QQbar)
After obtaining the partial_fraction_decomposition, I must apply RR to all numbers , but I do not manage to do it. Without that conversion, inverse_laplace won't workMon, 14 May 2018 12:03:38 -0500http://ask.sagemath.org/question/42341/is-it-possible-to-change-polynomial-ring-in-the-middle-of-a-computation/Answer by nbruin for <p>Hi, I'm trying to invert the "Pollaczek-Khinchine" Laplace transform when it is rational</p>
<p>This works for me at degree 2:</p>
<pre><code>var('x,s')
Fx = (1/6*exp(-2*x)+5/6*exp(-6*x));rho=2/3
print('Hyperexponential claims:',Fx)
R.<s> = PolynomialRing(QQbar)#when all coefficients are not integer, use CC
FF = R.fraction_field()
L_F=laplace(Fx,x,s)#Laplace transform of F
#Compute Pollackek-Khinchine (PK) formula L_rui for the Laplace transform (LT) of ruin probability
m1=L_F(s=0)
fe=L_F/m1
Fe=(1-fe)/s
L_rui=rho*Fe/(1-rho*fe)
show(L_rui.simplify_full())
inverse_laplace(SR(L_rui),s,u)
</code></pre>
<p>but not at degree 3, since I do not know how to use partial_fraction_decomposition, and then to switch to RR numbers and then invert .
If I start in R.<s> = PolynomialRing(RR), for an already known LT, everything is fine. But, a certain simplification by s in PK formula will become impossible due to rounding errors, so I am forced to start with R.<s> = PolynomialRing(QQbar)
After obtaining the partial_fraction_decomposition, I must apply RR to all numbers , but I do not manage to do it. Without that conversion, inverse_laplace won't work</s></s></p><s><s>
</s></s> http://ask.sagemath.org/question/42341/is-it-possible-to-change-polynomial-ring-in-the-middle-of-a-computation/?answer=42354#post-id-42354(suggestion: edit your question to include clarifying information rather than put it in an answer)
It seems homomorphisms for fraction fields are a bit underdeveloped. However, for polynomials things work fine and since `numerator` and `denominator` allow you to take apart a fraction field element, you can just assemble the result you want yourself.
sage: RRs = PolynomialRing(RR,'s')
sage: conv = lambda f: RRs(f.numerator())/RRs(f.denominator())
sage: [conv(c) for c in LL]
[0.533227478329869/(s + 1.08766810903621),
0.0780427896513687/(s + 2.96794729592164),
0.0137297320187619/(s + 5.57596354241057)]
sage: [inverse_laplace(conv(c),s,u) for c in LL]
[54528517/102261266*e^(-13571505/12477616*u),
3364341/43108928*e^(-163049898/54936925*u),
6233434/454009881*e^(-9174183/1645309*u)]
That last answer looks very fishy, though: you're inputting something with floats and you're getting back rational numbers. That can't be trustworthy. Something in the code has replaced floats by rationals that are close to it (that's what maxima does more often). But it means you don't know how to interpret those rational numbers.
Much more useful is
sage: var('a,b,s,u')
(a, b, s, u)
sage: inverse_laplace(a/(s+b),s,u)
a*e^(-b*u)
which shows you what the transform actually does on expressions of the form that you have.
Depending on what you want to do next, it could be that
def ilt(f):
if f.numerator().degree() > 0 or f.denominator().degree() >1 or f.denominator().leading_coefficient() != 1:
raise ValueError
return f.numerator()*exp(-f.denominator().constant_coefficient()*s)
sage: [ilt(f) for f in LL]
[0.5332274783298694?*e^(-1.087668109036214?*s),
0.07804278965136865?*e^(-2.967947295921641?*s),
0.01372973201876195?*e^(-5.575963542410567?*s)]
gives a more useful result. (Be careful: you do end up with an element in SR and elements of Qbar usually don't play nice with that. So you might STILL want to convert the relevant numbers to actual floats)Wed, 16 May 2018 10:49:12 -0500http://ask.sagemath.org/question/42341/is-it-possible-to-change-polynomial-ring-in-the-middle-of-a-computation/?answer=42354#post-id-42354Answer by florin for <p>Hi, I'm trying to invert the "Pollaczek-Khinchine" Laplace transform when it is rational</p>
<p>This works for me at degree 2:</p>
<pre><code>var('x,s')
Fx = (1/6*exp(-2*x)+5/6*exp(-6*x));rho=2/3
print('Hyperexponential claims:',Fx)
R.<s> = PolynomialRing(QQbar)#when all coefficients are not integer, use CC
FF = R.fraction_field()
L_F=laplace(Fx,x,s)#Laplace transform of F
#Compute Pollackek-Khinchine (PK) formula L_rui for the Laplace transform (LT) of ruin probability
m1=L_F(s=0)
fe=L_F/m1
Fe=(1-fe)/s
L_rui=rho*Fe/(1-rho*fe)
show(L_rui.simplify_full())
inverse_laplace(SR(L_rui),s,u)
</code></pre>
<p>but not at degree 3, since I do not know how to use partial_fraction_decomposition, and then to switch to RR numbers and then invert .
If I start in R.<s> = PolynomialRing(RR), for an already known LT, everything is fine. But, a certain simplification by s in PK formula will become impossible due to rounding errors, so I am forced to start with R.<s> = PolynomialRing(QQbar)
After obtaining the partial_fraction_decomposition, I must apply RR to all numbers , but I do not manage to do it. Without that conversion, inverse_laplace won't work</s></s></p><s><s>
</s></s> http://ask.sagemath.org/question/42341/is-it-possible-to-change-polynomial-ring-in-the-middle-of-a-computation/?answer=42350#post-id-42350The working code for n=3 is
A=[1/6,2/6,3/6];ex=[2,6,4];rho=5/8;var('x s u')
F=sum(A[i]*exp(-x*ex[i]) for i in [0..2])#sum of exponentials
R.<s> = PolynomialRing(QQbar)
FF = R.fraction_field()
L_F=sum(A[i] /(s+ex[i]) for i in [0..2]);#Laplace transform of F
#Compute Pollackek-Khinchine formula L_rui for the Laplace transform of ruin probability
m1=L_F(s=0)
fe=L_F/m1
Fe=(1-fe)/s
show(Fe) #note the s is simplified when working in QQbar
L_rui=rho*Fe/(1-rho*fe)
R.change_ring(RR)
FF = R.fraction_field()
wh,LL=L_rui.partial_fraction_decomposition()
show(LL[0])
with output
[0.5332274783298694?s+1.087668109036214?,0.07804278965136865?s+2.967947295921641?,0.01372973201876195?s+5.575963542410567?]
I need now to "remove" somehow the ? before taking
inverse_laplace(LL[0],s,u)
It would be nice if something like
inverse_laplace(LL[0].change_ring(RR),s,u)
or `inverse_laplace(RR(LL[0]),s,u)`
would workTue, 15 May 2018 03:55:02 -0500http://ask.sagemath.org/question/42341/is-it-possible-to-change-polynomial-ring-in-the-middle-of-a-computation/?answer=42350#post-id-42350Answer by tmonteil for <p>Hi, I'm trying to invert the "Pollaczek-Khinchine" Laplace transform when it is rational</p>
<p>This works for me at degree 2:</p>
<pre><code>var('x,s')
Fx = (1/6*exp(-2*x)+5/6*exp(-6*x));rho=2/3
print('Hyperexponential claims:',Fx)
R.<s> = PolynomialRing(QQbar)#when all coefficients are not integer, use CC
FF = R.fraction_field()
L_F=laplace(Fx,x,s)#Laplace transform of F
#Compute Pollackek-Khinchine (PK) formula L_rui for the Laplace transform (LT) of ruin probability
m1=L_F(s=0)
fe=L_F/m1
Fe=(1-fe)/s
L_rui=rho*Fe/(1-rho*fe)
show(L_rui.simplify_full())
inverse_laplace(SR(L_rui),s,u)
</code></pre>
<p>but not at degree 3, since I do not know how to use partial_fraction_decomposition, and then to switch to RR numbers and then invert .
If I start in R.<s> = PolynomialRing(RR), for an already known LT, everything is fine. But, a certain simplification by s in PK formula will become impossible due to rounding errors, so I am forced to start with R.<s> = PolynomialRing(QQbar)
After obtaining the partial_fraction_decomposition, I must apply RR to all numbers , but I do not manage to do it. Without that conversion, inverse_laplace won't work</s></s></p><s><s>
</s></s> http://ask.sagemath.org/question/42341/is-it-possible-to-change-polynomial-ring-in-the-middle-of-a-computation/?answer=42345#post-id-42345First, you should provide the whole code in degree 3 so that we see where the problem is !
That said, if `P` is an element of the polynomial ring `R` defined over `QQbar`, you can transform it into a (floating-point) complex polynomial as follows:
sage: P.change_ring(CC)
Mon, 14 May 2018 17:18:18 -0500http://ask.sagemath.org/question/42341/is-it-possible-to-change-polynomial-ring-in-the-middle-of-a-computation/?answer=42345#post-id-42345