Numerical instability in the inverse Laplace transform

168 views Asked by At

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)

enter image description here

0

There are 0 answers