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

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


edit retag close merge delete

Sort by » oldest newest most voted

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

more

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.

more

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


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


it only remains to parfrac this. Thanks, Florin

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

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

which University are you based in ?

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