Ask Your Question
1

inverse_laplace of a fraction whose denominator has real roots (or complex)

asked 2018-05-14 11:04:41 +0100

florin gravatar image

Hi, this is the simplest and very important method of Laplace inversion . It can be applied to any Laplace transform (LT), by starting with a Pade approximation of the LT, then partial fractions and inversion. For some reason, a program I had written last year stopped working

def Ruin(Fx, rho, x, u):#assumes rational survival function Fx, i.e. a hyperexponential density
    var('s')
    L_F=laplace(Fx,x,s) 
    m1=L_F(s=0) #some algebraic manipulations known as Pollaczek-Khinchine formula
    fe=L_F/m1
    Fe=factor((1-fe)/s)
    L_rui=rho*Fe/(1-rho*fe)
    C=ComplexField(53);
    dec=Frac(C['s'])(L_rui).partial_fraction_decomposition();
    n=len(dec)
    par=[inverse_laplace(dec[1][i],s,u) for i in [0..n]];
    rui=sum(par)
    return rui,  L_rui, fe

It's pretty easy to repair the program when roots are rational -- see question "partial_fraction_decomposition" with possibly "complex roots", again. For nonrational roots, I proposed there the following test case

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[0])
inverse_laplace(LL[0],s,u)

Thanks in advance :)

edit retag flag offensive close merge delete

2 Answers

Sort by » oldest newest most voted
1

answered 2018-05-14 13:30:47 +0100

Emmanuel Charpentier gravatar image

Okay. Now, we have two parts :

  1. How to get an exact partial fraction decomposition of L
  2. How to get ins inverse Laplace transform.

Let's be lazy and first solve the second one. We'll do it symbolically.

sage: var("a,b,c,d,e,s,t")
(a, b, c, d, e, s, t)
sage: Ls=(s-a)*(s-b)/((s-c)*(s-d)*(s-e))
sage: inverse_laplace(Ls,s,t)
(a*b - (a + b)*c + c^2)*e^(c*t)/(c^2 - c*d - (c - d)*e) - (a*b - (a + b)*d + d^2)*e^(d*t)/(c*d - d^2 - (c - d)*e) + (a*b - (a + b)*e + e^2)*e^(e*t)/(c*d - (c + d)*e + e^2)

which is $\frac{{\left(a b - {\left(a + b\right)} c + c^{2}\right)} e^{\left(c t\right)}}{c^{2} - c d - {\left(c - d\right)} e} - \frac{{\left(a b - {\left(a + b\right)} d + d^{2}\right)} e^{\left(d t\right)}}{c d - d^{2} - {\left(c - d\right)} e} + \frac{{\left(a b - {\left(a + b\right)} e + e^{2}\right)} e^{\left(e t\right)}}{c d - {\left(c + d\right)} e + e^{2}}$.

Our problem is now reduced to find exact (if possible closed_form) expressions of $a,~b,~c,~d$ and $e$, respective roots of the numerator and denominator of L.

The second part of my answer to your previous question shows how to obtain the radical expression of an element of F if such an expression can be computed by Sage.

This can be applied either to the result of L.partial_fraction_decomposition() or to the factorization of the numerator or denominator of L, with identification to the monomials into which its numerator and denominator can be factorized into.

But I'm not sure that such explicit representation of the resulting function is worthy...

Back to your original problem : I'll have a look at this (survival analysis is an important case in biostatistics). What do you want to show ? I'd also like to question the (implicit) assumption of interest in the (total) survival function : in competitive risks models, there are more interesting things to do...

edit flag offensive delete link more
0

answered 2018-05-14 11:10:44 +0100

florin gravatar image

updated 2018-05-14 17:07:35 +0100

in competitive risks models, there are more interesting things to do OK, I'm dealing with a risk expert :) which University are you based in?

I am just trying now to program in Sage a pedagogic exercise: compute the ruin probability for the Cramer-Lundberg model with a) hyperexponential claims; this should be easy, except that I am novice enough to Sage to have forgotten how these 1.33456? numbers are called, and how to convert an expression involving many of these numbers to a type which will be accepted by the command inverse_laplace

b) uniform claims (so the LT is not rational); after rationalizing the LT by Pade, this reduces to previous; it's a one liner in Mathematica, but takes hours in sage due to the crude help system.

edit flag offensive delete link more

Comments

The Pade of the uniform program is

var('s'); rho=2/3;L_f=(1-exp(-s))/(s) #LT of uniform density
L_F=(1-L_f)/s #LT of uniform survival function
#Compute Pollackek-Khinchine formula L_rui for the Laplace transform of ruin probability
m1=limit(L_F,s=0) #the s is removable singularity
fe=L_F/m1
Fe=(1-fe)/s 
L_rui=rho*Fe/(1-rho*fe)
show(L_rui)
n=2;
tay=taylor(L_rui,s,0,2*n+1)
ta=tay.power_series(QQ)
Lp=ta.pade(1,2)
Lp

with answer

(5/8*s + 15/2)/(s^2 + 45/4*s + 45/4)

it only remains to parfrac this. Thanks, Florin

florin gravatar imageflorin ( 2018-05-14 17:35:58 +0100 )edit

the conversion is doable by RR(). Is there an easy way to apply RR to all parts of an expression?

florin gravatar imageflorin ( 2018-05-14 17:59:50 +0100 )edit

OK I got it, I just need to replace QQbar by RR. Silly me !!!!

florin gravatar imageflorin ( 2018-05-14 18:07:14 +0100 )edit

which University are you based in ?

I'm working for the Assistance Publique - Hôpitaux de Paris (i. e. University Hospitals in Paris).

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2018-05-15 11:05:01 +0100 )edit

Your Answer

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

Add Answer

Question Tools

Stats

Asked: 2018-05-14 11:04:41 +0100

Seen: 589 times

Last updated: May 14 '18