# Numerically obtaining the inverse Laplace transform of data

I have been using several Mathematica packages to do numerical inverse Laplace transforms on known (expressible in closed form) expressions, $\tilde{f}(s)$. I am now being confronted with the more difficult problem where I can only obtain $\tilde{f}(s)$ as a table like

{{s_1,f_1},{s_2, f_2}, {s_3, f_3 ... etc.}


I have not been able to find any packages for doing the inverse Laplace transform on tabulated data. I saw one ancient discussion in a forum that says the inverse Fourier transform can be used but I don't see how that's possible.

I would greatly appreciate if anyone had done something like this before, or knew of an appropriate package, or had any suggestions for a workaround.

EDIT: this is how I'm getting tabulated results for the transformed function: I have

$\tilde{f}(s) =\frac{1}{s} \frac{1-\tilde{g}(s)}{1-a \tilde{g}(s)}$ where $\tilde{g}(s) = \int_{-\infty}^0 h(s,x) dx$. I know what $h(s,x)$ is, but the integral over $x$ cannot be done analytically so I am doing it numerically which is why I end up with tabular data for $\tilde{f}(s)$.

-
Is it common to get tabulated data of the laplace transform? I find it weird –  Rojo Apr 26 '12 at 16:26
@R.M and is the problem of finding a numerical approximation of a sampled Z-transform's inverse Z-transform easier? –  Rojo Apr 26 '12 at 16:36
@Rojo I have edited the question to show why I am getting tabulated data. Maybe there is a way to get around it, but I don't see it –  BeauGeste Apr 26 '12 at 16:36
@BeauGeste. Well, if you can choose which points to get, and that integral converges in the imaginary s... Perhaps an inverse dft? –  Rojo Apr 26 '12 at 16:39
It might be better to state your problem instead as "find the inverse Laplace transform of a black-box function". Are you able to evaluate your black-box function at complex arguments, or only real ones? –  Ｊ. Ｍ. Apr 27 '12 at 14:12

Here is my attempt at an answer - I had to make up an example, and obviously much of what follows is dependent on details of this example.

Edit

However, what I believe this example shows quite clearly is that a finite set of tabulated data at discrete points does not suffice to guarantee a good inverse Laplace transform, because the analytic structure of the function that interpolates between these data points is not uniquely determined by a finite number of points. You'll need additional information beyond the data table to get a correct inverse transform.

End edit

Let's start with a known function and its known Laplace transform, so we have something to compare the numerical results to:

originalFunction[t_] := t^4 Sin[t];
l[s_] := Evaluate[LaplaceTransform[originalFunction[t], t, s]];


Now I sample the Laplace transform l at discrete points to simulate the data that would be the given quantities of the problem:

data = Table[{s, l[s]}, {s, -5, 5, .1}];


The numerical inversion of this Laplace transform now can be performed by assuming a fit to the data that has a sufficiently simple functional form that allows us to do the inversion. I'll make the ansatz that we can fit the data with a rational function, leaving the degree of the numerator and denominator as parameters that may have to be adjusted by trial and error:

fit[s_] := Evaluate[
Block[{numeratorN = 5, denominatorN = 8},
rationalFunction =
Total[Array[a, numeratorN + 1] s^Range[0, numeratorN]]/
Total[Array[b, denominatorN + 1] s^Range[0, denominatorN]];
rationalFunction /.
FindFit[data, rationalFunction,
Join[Array[a, numeratorN + 1], Array[b, denominatorN + 1]], s]]];


This fit function looks promising if we check the plot versus the data points:

Show[ListPlot[data, PlotRange -> All],
Plot[fit[s], {s, -5, 5}, PlotRange -> All]]


So we might feel somewhat confident that the inversion will work as follows:

fitInverse[t_] := Evaluate[InverseLaplaceTransform[fit[s], s, t]];


Unfortunately, the result isn't too good:

Plot[{originalFunction[t], Re@fitInverse[t]}, {t, -1, 1},
PlotStyle -> {Directive[Thick, Red], Directive[Blue, Dashed]}]


The red curve is the original function that we're trying to recover with the inverse Laplace transform. This is a problem that will be hard to avoid in practice. In my example, I can indeed make the fit work out much better by just increasing the degree of the denominator to get a faster fall-off at infinity:

fit[s_] := Evaluate[
Block[{numeratorN = 5, denominatorN = 10},
rationalFunction =
Total[Array[a, numeratorN + 1] s^Range[0, numeratorN]]/
Total[Array[b, denominatorN + 1] s^Range[0, denominatorN]];
rationalFunction /.
FindFit[data, rationalFunction,
Join[Array[a, numeratorN + 1], Array[b, denominatorN + 1]], s]]];
Show[ListPlot[data, PlotRange -> All],
Plot[fit[s], {s, -5, 5}, PlotRange -> All]]


In the plot interval chosen here, you won't see any difference in the fit, compared to the first attempt. But the small change is enough to make the inversion work:

fitInverse[t_] := Evaluate[InverseLaplaceTransform[fit[s], s, t]];
Plot[{originalFunction[t], Re@fitInverse[t]}, {t, -1, 1},
PlotStyle -> {Directive[Thick, Red], Directive[Yellow, Dashed]}]


The numerical inverse and the target function lie on top of each other now. So that's how you can do it in principle: Get a good fit to the data that also has the "correct" asymptotic behavior (i.e., falls off in the way you expect based on the physics or other knowledge about your data). Then do the InverseLaplaceTransform on that fit and hope for the best.

-

I have two suggestions. You may:

1. Establish a function that contains the integration $g(s)$. The bad news is that you may need a lot of time to evaluate the integration. In my opinion, in the inverse Laplace calculation, we don't have to give the analytical expression of $F(S)$. If I were you, I may write the function like:

F[s_]:=(g = NIntegrate[h[s,x],{x,-Infinity,0}];(1-g)/(s*(1-a*g)))


I tried this method. If I set h(x,s)=s/(x^3 - 1),doing NIntegrate[h[x,3+4I],{x,-Infinity,0}] for 10000 times cost me 47.6 seconds:

h[x_, s_] := s/(x^3 - 1);
Do[b = NIntegrate[h[x, 3 + 4 I], {x, -Infinity, 0}];, {i, 1, 10000}] // Timing
(* {47.611999999999995, Null} *)


If you use this method, i think the total time for deploying the NIntegrate will drive you crazy! So I use Compaq Visual Fortran and numerical libraries IMSL V4.0 and data Visualization library Matfor V4.1 to give you the exampe in order to save time. Im my example, I set h(x,s)=s/(x^3 - 1), and you may revise my program below.

MODULE GUANGHUI_HE
IMPLICIT NONE
COMPLEX(8) S1

CONTAINS
Function F(S)
USE IMSL !this is numerica libraries IMSL V4.0
IMPLICIT NONE
COMPLEX(8) S, F, g
REAL(8) :: BOUND, ERRABS, ERRREL, RESULT1(2), ERREST
S1 = S
BOUND  = 0d0    !integration interval end point 0
ERRABS = 1d-10  !absolute error set for subroutine 1E-10
ERRREL = 1d-10  !relative error set for subroutine 1E-10
CALL DQDAGI (FR, BOUND, -1, ERRABS, ERRREL, RESULT1(1), ERREST)
!this subroutine DQDAGI is one of imsl for Numerical Integration,
!see http://www.roguewave.com/products/imsl-numerical-libraries/fortran-library.aspx
CALL DQDAGI (FI, BOUND, -1, ERRABS, ERRREL, RESULT1(2), ERREST)
g = RESULT1(1) + (0D0, 1D0)*RESULT1(2)
F = (1d0 - g)/(S*(1D0 - 2D0 * g))
END FUNCTION F

FUNCTION FR(X)
IMPLICIT NONE
REAL(8) :: X, FR
FR = DBLE( S1/(X**3 - 1d0) ) ! "S1/(X**3 - 1d0)" this stands for your h(x,s),you may revise it as you like. Function FR is Real part of h(x,s)
END FUNCTION FR

FUNCTION FI(X)
IMPLICIT NONE
REAL(8) :: X, FI
FI = DIMAG( S1/(X**3 - 1d0) ) ! "S1/(X**3 - 1d0)" this stands for your h(x,s),you may revise it as you like. Function FI is Imaginary part of h(x,s)
END FUNCTION FI

END MODULE GUANGHUI_HE

! now let's start to deploy the tedious F(S), in the main program
PROGRAM MAIN
USE GUANGHUI_HE
USE IMSL, ONLY : DINLAP !module DINLAP is one of IMSL, used for LaplaceInverse
USE fml ! this module is one of Matfor, for data process
use fgl ! this module is one of Matfor, for data Visualization, I will use it to plot below.

IMPLICIT NONE
INTEGER, PARAMETER :: N=2000
!N is the time point number i want to compute.
REAL(8) :: T(N), FINV(N), Pi=3.141592654
!Array T and Finv have length N, used for restoring time and time domain value of F(s) respectively.
REAL(4) :: T0,T1 !T0 and T1 are set for recording time spot.
REAL(8) :: ALPHA=0D0, RELERR=1D-6
INTEGER :: KMAX=10000, I
Type(mfArray) :: V
DO I=1,N
T(I) = 1D-2 + 1D-2 * DBLE(I)
END DO

CALL CPU_TIME(T0)
CALL DINLAP (F, N, T, ALPHA, RELERR, KMAX, FINV)
CALL CPU_TIME(T1)
WRITE(*,*) "Total CPU TIME FOR Laplace Inverse is",T1-T0,"Seconds!"
V = 1D0 - 0.5d0 * mfExp(-((3D0*DSqrt(3d0)* mf(T) )/(4*Pi))) !(1 + (2 Pi s)/(3 Sqrt[3]))/(s (1 + (4 Pi s)/(3 Sqrt[3]))) is the analytical soulution
CALL msPlot(mf(T), mf(Finv))
call msHold('on')
call msPlot(mf(T),V,"r0--")
CALL msViewPause()

END


2. Get the 2D grid data about $\Re(S)$ and $\Im(S)$, to create the 2D interpolation function and treat it as the analytical one, being called by the Laplace inverse subroutine. We can create a function like (pseudocode):

function F(S)
complex S,F
real A, B
A=Re(S); B=Im(S)
F=Interpolation(A, B) (* This is why we need the 2D interpolation function *)
end function

-
could you show a simple example of each? –  acl May 12 '12 at 21:05