# 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. = 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. = 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

edit retag close merge delete

Sort by » oldest newest most voted

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

more

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


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,s,u)


It would be nice if something like

inverse_laplace(LL.change_ring(RR),s,u)


or inverse_laplace(RR(LL),s,u) would work

more

First, 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)

more