Ask Your Question
0

Is it possible to change Polynomial Ring in the middle of a computation?

asked 2018-05-14 19:03:38 +0200

florin gravatar image

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

3 Answers

Sort by » oldest newest most voted
0

answered 2018-05-15 00:18:18 +0200

tmonteil gravatar image

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

answered 2018-05-15 10:55:02 +0200

florin gravatar image

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[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 work

edit flag offensive delete link more
0

answered 2018-05-16 17:49:12 +0200

nbruin gravatar image

updated 2018-05-16 19:53:19 +0200

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

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-05-14 19:03:38 +0200

Seen: 551 times

Last updated: May 16 '18