Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
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
1  
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? –  J. M. Apr 27 '12 at 14:12

2 Answers 2

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
    

    the black line stands for numerical solution, and red one is exact one! They almost overlap each other.

  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
    
share|improve this answer
    
could you show a simple example of each? –  acl May 12 '12 at 21:05

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

fit 1

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]}]

inverse 1

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

fit 1

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]}]

inverse 2

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.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.