2
$\begingroup$

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])$

enter image description here

$\endgroup$
6
  • 1
    $\begingroup$ 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. $\endgroup$
    – EMP
    Mar 23 at 17:03
  • $\begingroup$ 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. $\endgroup$
    – njuffa
    Mar 23 at 19:11
  • $\begingroup$ @njuffa jup around 35 $\endgroup$ Mar 24 at 14:10
  • $\begingroup$ @EMP so I did it $\endgroup$ Mar 24 at 14:11
  • 3
    $\begingroup$ 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. $\endgroup$
    – user14717
    Mar 24 at 17:58

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.