I have a problem with Laplace inversion and my function is not numerically stable for the Laplac 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 instab
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)