Numerical Mellin and Two-Sided Laplace Inversion

There are only a handful of methods for numerical Mellin (or equivalently, two-sided Laplace) inversion that have survived due to their computational accuracy and efficiency. One of these such methods is accredited to Dishon and Weiss [1]. The method is quite accurate despite being a numerical approximation and requiring the evaluation of an infinite series.

We implement example 1 in their paper below:


#### DISHON AND WEISS ###

#packages
library(fAsianOptions)

#global variables
N=100
x=seq(0.1,10,by=0.1)
k=1:N
t=-log(x)
a=2
T1=20

#function definitions
F_s <- function(s){(2^(-s))*cgamma(s)}
f_x <- function(x){exp(-2*x)}

#DW routine
DW <- function(t){
step <- function(r){Re(F_s(a+pi*1i*r/T1)*cos(pi*r*t/T1))-Im(F_s(a+pi*1i*r/T1)*sin(pi*r*t/T1))}
f_a <- exp(a*t)/(2*T1)*(F_s(a)+ 2*sum(sapply(k,step)))
return(Re(f_a))}

#main
V <- sapply(t,DW)
U <- sapply(x,f_x)
plot(V, xlab='x', main="f(x)=exp(-2x)")
lines(V, col='blue')
lines(U, col='red')
plot(log(abs(V-U)), ylab="Log Abs Error", xlab='x', main="Error via Inversion: f(x)=exp(-2x)")
lines(log(abs(V-U)), col='blue')

dwactdwerror

 

The other examples in [1] can be inverted by changing the functions above.

[1] Dishon, M., and G.H. Weiss, 1978: Numerical inversion of Mellin and two-sided Laplace transforms. J. Comput. Phys., 28, 129-132, doi:10.1016/0021-9991(78)90053-0.