# Numerical instability in the inverse Laplace transform

I have a problem with Laplace inversion and my function is not numerically stable for the Laplace inverse, but I do not understand the cause of this problem. Here is my code and graph of this problem. Does anyone have any idea? I have also tried all 3 Dehoog Stehfest and Talbot, they all have a different type of instabillity

import mpmath as mp
import numpy as np
import matplotlib.pyplot as plt
def Laplace_Max(s):
#mp.dps = 100
A= s**(5/4)/ (
0.00480931 +
0.0244077*s**(1/4) +
0.0129056*mp.exp(-35*s)*s**(1/4) +
0.00329167*mp.exp(0.707997*s)*s**(1/4) * mp.gammainc(0.0, 35.708*s,mp.inf, regularized=True) -
0.00530593*mp.gammainc(1.25, 35*s,mp.inf,   regularized=True)
)
return A
t_R_inv = np.linspace(1,1200,12000)
Max_R=np.zeros(len(t_R_inv))
for i in range(len(t_R_inv)):
Max_R[i]=mp.invertlaplace(Laplace_Max, t_R_inv[i], method = 'stehfest', dps = 100, degree = 1000)
plt.plot(Max_R)


Update The inverse laplace of this function!

$$\frac{1}{s}^{\frac{5}{4}}\cdot (0.00480931 + 0.02440766 \cdot s^{\frac{1}{4}} + 0.0129056 \cdot e^{-35 \cdot s} \cdot s^{\frac{1}{4}}+ 0.003291670 \cdot e^{0.7079967\cdot s} \cdot s^{\frac{1}{4}} \cdot \Gamma[0, 35.7079967 \cdot s] - 0.00530593 \cdot\Gamma[1.25, 35\cdot s])$$

• Please write mathematically what the problem is that you are trying to solve and what these solution methods you mention are that way people have an easier time helping you.
– EMP
Mar 23, 2022 at 17:03
• From the graph it is hard to see where exactly the region of instability is. It seems to be around $t=30$? Please provide a numerical range for this. Mar 23, 2022 at 19:11
• @njuffa jup around 35 Mar 24, 2022 at 14:10
• The inverse Laplace transform is ill-conditioned: dl.acm.org/doi/abs/10.1016/j.amc.2013.03.112; maybe you can fight it a bit by (say) high precision. Mar 24, 2022 at 17:58
• If the problem is ill-conditioned, you can't really expect that different methods are gonna help. Mar 24, 2022 at 18:00