Laplace transform using numerical integration in Python has very poor precision

Question

I have written a function to compute the Laplace transform of a function using scipy.integrate.quad. It is not a very sophisticated function and currently performs poorly on the probability density function of an Erlang distribution.

I have included all my work below. I first compute the Laplace transform and then the inverse in order to compare it to the original p.d.f. of the Erlang. I use mpmath for this. The mpmath.invertlaplace is not the problem as it manages to convert the closed-form Laplace transform back to the original p.d.f. quite perfectly.

Please help me understand what the problem is with my numerical laplace transform. I get the following error but have not been able to resolve it.

IntegrationWarning: The occurrence of roundoff error is detected, which prevents the requested tolerance from being achieved. The error may be underestimated. a=0,b=np.inf,limit=limit)[0]

PLOT

CODE

import numpy as np
import matplotlib.pyplot as plt
import math as m
import mpmath as mp
import scipy.stats as st
from scipy.integrate import quad


def get_laplace(func,limit=10000):
    
    '''
    Returns laplace transfrom function
    '''
    
    def laplace(s):
        '''Numerical laplace transform'''
        # Seperate into real and imaginary parts
        x = np.real(s)
        y = np.imag(s)
        
        def real_func(t):
            return m.exp(-x*t)*m.cos(y*t)*func(t)
        
        def imag_func(t):
            return m.exp(-x*t)*m.sin(y*t)*func(t)
        
        real_integral = quad(real_func,
                             a=0,b=np.inf,limit=limit)[0]
        
        imag_intergal = quad(real_func,
                             a=0,b=np.inf,limit=limit)[0]
        
        return complex(real_integral,-imag_intergal)
    
    return laplace
        

def L_erlang(s,lam,k):
    '''
    Closed form laplace transform of Erlang or Gamma distribution.
    '''
    return (lam/(lam+s))**k


if __name__ == '__main__':
    
    # Setup Erlang probability density function
    k = 5
    lam = 1
    pdf = st.erlang(a=k,scale=1/lam).pdf
    
    # Laplace transforms
    Lnum = get_laplace(pdf)          # numerical approximation
    L = lambda s: L_erlang(s,lam,k)  # closed form
    
    
    # Use mpmath library to perform inverse laplace
    
    # Invserse transfrom on numerical laplace function
    invLnum = lambda t: mp.invertlaplace(Lnum,t,
                                         method='dehoog',
                                         dps=5, 
                                         degree=3)
    
    # Invserse transfrom on closed-form laplace function
    invL = lambda t: mp.invertlaplace(L,t,
                                      method='dehoog',
                                      dps=5, 
                                      degree=3)
    
    # Grid to visualise
    T = np.linspace(0.1,10,10)
    # Get values of inverse transforms
    lnum = np.array([invLnum(t) for t in T])
    l = np.array([invL(t) for t in T])
    
    # Plot
    plt.plot(T,lnum,label='from numerical laplace')
    plt.plot(T,l,label='from closed-form laplace')
    plt.plot(T,pdf(T),label='original pdf',linestyle='--')
    plt.legend(loc='best')
    plt.show()

UPDATE

After two cups of VERY strong coffee, I managed to see the obvious mistake and make the code work. It quite embarrassing actually. Have a look at this line of code:

imag_intergal = quad(real_func,
                     a=0,b=np.inf,limit=limit)[0]

Hmmm, real_func hey? So it should read:

imag_intergal = quad(imag_func_func,
                     a=0,b=np.inf,limit=limit)[0]

The one gets this lovely plot:

Conclusion So why go through all the trouble to numerically perform the Laplace transform of something that we have a closed form solution for. That is because the interest lies somewhere else. We do not have a closed expression for the conditional future lifetime distribution which is similar to the hazard function. Lets call it h. Then for the Erlang distribution erl = st.erlang(a=k,scale=1/lam) that has been active for tau units of time we have h = lambda t: erl.pdf(t+tau)/erl.sf(tau). This distribution can be used as a holding time in a Semi-Markov Model (SMP). To analyse the transient behaviour of the SMP one use Laplace transforms. Usually only pdfs are used but now I can use hazard functions. Its pretty cool because it means one can model transient behaviour without assuming everything is new.

Answer 1

Good day. I am repeating the Updates section from the original questions as this is the solution to the questions. This way, the questions can be marked as resolved.

UPDATE

After two cups of VERY strong coffee, I managed to see the obvious mistake and make the code work. It quite embarrassing actually. Have a look at this line of code:

imag_intergal = quad(real_func,
                     a=0,b=np.inf,limit=limit)[0]

Hmmm, real_func hey? So it should read:

imag_intergal = quad(imag_func_func,
                     a=0,b=np.inf,limit=limit)[0]

The one gets this lovely plot:

Conclusion So why go through all the trouble to numerically perform the Laplace transform of something that we have a closed form solution for. That is because the interest lies somewhere else. We do not have a closed expression for the conditional future lifetime distribution which is similar to the hazard function. Lets call it h. Then for the Erlang distribution erl = st.erlang(a=k,scale=1/lam) that has been active for tau units of time we have h = lambda t: erl.pdf(t+tau)/erl.sf(tau). This distribution can be used as a holding time in a Semi-Markov Model (SMP). To analyse the transient behaviour of the SMP one use Laplace transforms. Usually only pdfs are used but now I can use hazard functions. Its pretty cool because it means one can model transient behaviour without assuming everything is new.