Click here to Skip to main content
Click here to Skip to main content

Numerical Laplace Transforms and Inverse Transforms in C#

, 12 Apr 2008 CPOL
Rate this:
Please Sign up or sign in to vote.
Presents a C# class for calculating Laplace tranforms and inverse transforms.

Contents

Introduction

In modeling many physical processes described by differential or partial differential equations, solutions can be obtained analytically using Laplace transforms. For most of the more interesting problems, however, inverting the solution is often a problem. Additionally, there are other classes of problems where manipulation of Laplace transforms is simple, but the function of interest is not easily transformed. This article presents a C# implementation of the Stehfest algorithm for Laplace transform inversion as well as an algorithm to compute the Laplace transform integral, numerically, for an arbitrary function.

The source code and a test program for checking the Laplace transform calculations are available from the links above. A screenshot of the test program is shown below.

LaplaceTest.jpg

Why Laplace Transforms?

Many physical processes are described by the diffusion equation including temperature distributions, diffusion of chemicals, and fluid flow in porous media. My particular interest is in modeling oil reservoirs, so the diffusion or diffusivity equation is of fundamental importance. In cases where a reservoir can be considered uniform and homogeneous, the time variable in the partial differential equation can be transformed using Laplace transforms, eliminating the time derivatives, and converting a partial differential equation into an ordinary differential equation, which is fairly easily solved. Problems involving Laplace transforms also arise in computing statistical functions and other areas of math.

For reference, the diffusivity equation used for radial flow of a single fluid near a well in a homogeneous reservoir can be written as follows, where pD represents the dimensionless pressure, rD represents the dimensionless radial distance from a well, and tD represents dimensionless time. To actually solve the equation is beyond the scope of this article; however, a suitable boundary and initial conditions are needed, and the result is typically the Laplace transform of the pressure distribution in an oil reservoir. For anyone interested in oil field applications, please see the classic paper by Van Everdingen and Hurst.

DiffusivityEquation.gif

In other cases, we may already know the solution of the equation for some ideal situation, for example, when the flow from a well is constant. In the real world, however, things are rarely constant, and to determine the actual solution, it is necessary to calculate the convolution of the ideal solution with a boundary condition that varies as a function of time. In Laplace transforms, that simply involves multiplying the Laplace transforms of two functions, but again, the result is the Laplace transform of the desired solution, and we need to convert that back to real values. For reference, the convolution consists of an integral, expressed as follows:

Convolution.gif

In addition, there are cases where an ideal solution cannot be determined explicitly using Laplace transforms, but requires the use of finite difference approximations or other numerical techniques, for example, to account for some strange geometry, irregular boundaries, or other strangeness. In these cases, we might have the basic solution in terms of some tabulated results or a fitted curve, and then need to calculate the convolution to adjust for variations not included in the numerical model. Note that in complex cases, calculating the basic solution may take hours of computer time, so running a new finite difference calculation isn't an option, especially if we need to display interactive graphics. An example of this is in the paper by Fair (1981).

A Laplace transform is an integral transform. The transform and the corresponding inverse transform are defined as follows:

LaplaceTransform.gif

A complete description of the transforms and inverse transforms is beyond the scope of this article. Anyone needing more information can refer to the "bible" of numerical mathematics, Abramowitz and Stegun (1970). A fairly complete tabulation of analytical Laplace transforms is also included.

Since the most common programming problem involves taking inverse Laplace transforms, I'll discuss that first, then tackle the problem of numerically calculating the Laplace transform integral.

The Inverse Laplace Transform

As you can see from the equation defining the inverse Laplace transform, direct calculation using brute force is formidable, because it involves calculating a complex path integral. Fortunately, there are some efficient numerical methods available for computing the inverse transform. Of these methods, one of the best was presented by Harald Stehfest (1970) and has become known as the Stehfest Algorithm. The algorithm itself is documented in the reference, but contains a typographical error that was later corrected. Stehfest points out that the algorithm requires an even number of coefficients, and is limited to functions that do not oscillate rapidly and that have no discontinuities. Fortunately, solutions to many practical equations rarely suffer from those problems. A web page [^] maintained by Dr. Peter Valko also shows some other algorithms that may be useful, especially when the Stehfest Algorithm cannot be applied. The equations used for the Stehfest Algorithm for Laplace transform inversion are as follows:

StehfestFormula.gif

My implementation of the Stehfest Algorithm is fairly straightforward. Since the coefficients, Vi, do not depend on the function to be inverted, they can all be calculated once before doing any inversions. The C# code is shown below.

class Laplace 
{ 
    const int DefaultStehfest = 14; 
    public delegate double FunctionDelegate(double t); 
    static double[] V; // Stehfest coefficients 
    static double ln2; // log of 2
 
    public static void InitStehfest(int N) 
    { 
        ln2 = Math.Log(2.0); 
        int N2 = N / 2; 
        int NV = 2 * N2; 
        V = new double[NV]; 
        int sign = 1; 
        if ((N2 % 2) != 0) 
            sign = -1; 
        for (int i = 0; i < NV; i++) 
        { 
            int kmin = (i + 2) / 2; 
            int kmax = i + 1; 
            if (kmax > N2) 
                kmax = N2; 
            V[i] = 0; 
            sign = -sign; 
            for (int k = kmin; k <= kmax; k++) 
            { 
                V[i] = V[i] + (Math.Pow(k, N2) / Factorial(k)) * (Factorial(2 * k) 
                     / Factorial(2 * k - i - 1)) / Factorial(N2 - k) 
                     / Factorial(k - 1) / Factorial(i + 1 - k); 
            } 
            V[i] = sign * V[i]; 
        } 
    }
 
    public static double InverseTransform(FunctionDelegate f, double t) 
    { 
        double ln2t = ln2 / t; 
        double x = 0; 
        double y = 0; 
        for (int i = 0; i < V.Length; i++) 
        { 
            x += ln2t; 
            y += V[i] * f(x); 
        } 
        return ln2t * y; 
    } 

    public static double Factorial(int N) 
    { 
        double x = 1; 
        if (N > 1) 
        { 
            for (int i = 2; i <= N; i++) 
                x = i * x; 
        } 
        return x; 
    } 
}

I've found that a default of N = 14 coefficients gives reasonable results in C#, but depending on the hardware you use and the specific functions being inverted, more or less coefficients may yield better accuracy. Stehfest notes that the optimum number of coefficients mainly depends on the number of digits used in the floating point calculations, and above a certain limit, the numerical truncation errors start to swamp any increase in accuracy. For the test project, N = 14 yields between 7 and 8 decimal digits of precision.

To implement the inverse Laplace transform, I defined a static class Laplace. The static constructor calculates the default number of Stehfest coefficients the first time the class is referenced. Private static members of the class are used to save the Stehfest coefficients, V, and also the value of Ln(2) to avoid having to recalculate it every time. In case it's necessary to change the number of coefficients at run time, the InitStehfest(int N) method can be called to recalculate. Note that there is also a Factorial(int N) method included, since that is needed to calculate the Stehfest coefficients.

The function containing the Laplace transform is defined as a delegate, so that it can be passed as a parameter to the InverseTransform(FunctionDelegate f, double t) method that calculates the inverse transform at the specified value of t. As an example of using the transform inversion, the following code snippet defines a Laplace transform function equal to Ln(s)/s, then calculates the inverse transform at a value of t = 2.

    double f(double s) 
    { 
        return -Math.Log(s)/s; 
    }
     ... 
    double LinvCalc = Laplace.InverseTransform(f, 2.0);

The Laplace Transform

To calculate the Laplace transform of a function, it's necessary to calculate the Laplace integral defined above. The main problems involved lie in two areas:

  1. the integral has an upper limit of infinity and
  2. the exponential function inside the integral requires special consideration if accuracy is to be maintained.

Since directly calculating infinite integrals is impossible on a finite computer, the first "trick" I use is to change the Laplace transform integral into an equivalent finite integral by substituting u = e-t. Using this variable transformation, when t = 0, u = 1 and when t approaches infinity, u approaches 0. This is shown in the following equations:

VariableTransform.gif

Essentially, what I've done is map the interval between 0 and infinity onto the interval between 0 and 1, thereby avoiding having to deal with infinity. If you take a close look at the equation, though, you'll note that there's still a potential problem, since it requires evaluating Ln(u = 0) and also u(s - 1) when u = 0. If s < 1, there's also the possibility of a singularity. Fortunately, evaluating Ln(u = 0) isn't a real problem, since if the Laplace transform exists, the limiting value at t = infinity or u = 0 must be exactly 0. The problem with a singularity for s < 1 may, however, still be a problem. Evaluating Laplace transforms at small values of s, corresponding to infinite times, is a notoriously difficult problem that I've chosen to ignore, since I rarely need to work with such small values of s.

My implementation of the transform calculation is fairly straightforward. As shown in the following code, I added a double Transform(FunctionDelegate F, double s) method that calculates the integral using the supplied function, F, at the specified value of the Laplace parameter, s. The formula Math.Pow(u, s - 1) * F(-Math.Log(u)) corresponds to the transformed integral kernel, and u is the independent variable. The while loop implements Simpson's Rule[^] to calculate the integral, but for efficiency, I've written the algorithm using du as half of the interval usually shown in standard Simpson's rule formulas. This allowed me to factor a 2 out of the loop, and multiply once at the end, rather than every loop iteration. When you're dealing with large loops and trying to maintain precision, every little bit helps.

    public static double Transform(FunctionDelegate F, double s)
    {
        const int DefaultIntegralN = 5000;
        double du = 0.5 / (double)DefaultIntegralN;
        double y =  - F(0) / 2.0;
        double u = 0;
        double limit = 1.0 - 1e-10;
        while (u < limit)
        {
            u += du;
            y += 2.0 * Math.Pow(u, s - 1) * F(-Math.Log(u));
            u += du;
            y += Math.Pow(u, s - 1) * F(-Math.Log(u));
        }
        return 2.0 * y * du / 3.0;
    }

Note also that I've set the loop limit to 1.0 - 1e-10, rather than 1.0. When dealing with a sum of small numbers, numerical precision can haunt you, so I've learned to never count on the sums of floating point numbers being exactly equal to anything.

Using the Laplace transform method is simple. Simply define a function that returns a value, given a value of the variable, t, then call the Transform method to retrieve the Laplace transform for any given value of the Laplace variable, s. The following code snippet shows how to compute the Laplace transform of the function e-2t for a Laplace variable value of s = 2.

    double F(double t)
    {
        return 1.0 / Math.Exp(2.0*t);
    }
     ... 
    double LCalc = Laplace.Transform(F, 2);

After running some tests, I concluded that 5000 intervals in the Simpson's Rule calculation are enough to ensure good precision. If you play with this algorithm, you may want to check that. In the case of "strange" functions, you may need to use more divisions to avoid losing precision. Conversely, you may be able to use fewer intervals and gain some speed. For the test case shown in the example program, 5000 integral divisions yield around 9 or 10 decimal digits of precision.

Comments

I hope these algorithms will be of interest to someone and keep them from having to reinvent the wheel. Using analytical solutions to physical problems is still often a viable method, even though the current trend is to move more towards finite difference methods. Analytical solutions, however, are often orders of magnitude faster and more accurate. An example of a commercial application using these techniques can be seen on my website[^].

If anyone has any ideas for improving the performance and flexibility of the algorithms presented in this article, I'd be very pleased to hear about them. Specifically, there could be improvements when calculating the Laplace transforms of functions with singularities and discontinuities, although those aren't usually associated with physical processes.

To me, it's interesting to see how things have changed over the years. I originally implemented these algorithms back before there were PCs, and a Laplace transform inversion took about 10 minutes on a TI-99 programmable calculator. Later, the same algorithms were coded in Turbo Pascal 1.0, and would run in less than a minute on an 8086. Now, I run thousands per second on my laptop. To be fair, though, I have to point out that the main improvements have come by optimizing the functions, since many of the problems of interest to me involve Modified Bessel functions. As a point of reference, engineers who worked on the calculations presented in the Van Everdingen and Hirst paper have told me that to compute the functions tabulated in the paper, they had a bank of over a dozen secretaries armed with tables of mathematical functions and mechanical calculators, and they could turn out about one column in the tables per day. In comparison, I now turn out the equivalent of one column in a few milliseconds.

If those were the good old days, no thank you!

References

  • Van Everdingen, A. F. and W. Hurst, "The Application of the Laplace Transformation to Flow Problems in Reservoirs," Trans., AIME (194) 186, 305-324.
  • Fair, Walter B. Jr., "Pressure Buildup Analysis with Wellbore Phase Redistribution," SPEJ (April 1981), 259-270.
  • Abramowitz, Milton, and Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, Washington, 1970.
  • Stehfest, H., "Algorithm 368 - Numerical Inversion of Laplace Transforms," Communications of the ACM, (1970) 13, 47.
  • Valko, Peter, Texas A&M University.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Walt Fair, Jr.
Engineer Comport Computing
United States United States
Walt has been playing with software since around 1967 and has generated more runtime errors than the average village idiot. He is a CEO, Petroleum Engineer, software consultant, janitor, and now a graduate student again. Rather than sleep, he also plays with algorithms and systems for technical computing, develops software for engineering evaluations and is an avid amateur radio operator.
 
Walt was admitted back to UT Austin and is actually attempting to complete a PhD in engineering, thereby proving that he is crazier than the average old fart.
 
And now UT has gone and admitted Walt to PhD candidacy, proving that old guys can still ... what was he doing again?

Comments and Discussions

 
QuestionWas of much help Pinmemberabukamaranga1-Jun-14 12:04 
AnswerRe: Was of much help PinsubeditorWalt Fair, Jr.1-Jun-14 14:56 
AnswerRe: Was of much help PinmemberAlex Konnen1-Jun-14 23:07 
GeneralRe: Was of much help PinsubeditorWalt Fair, Jr.2-Jun-14 4:15 
QuestionCompiling and executing PinmemberDerekEccles14-Jan-13 9:40 
AnswerRe: Compiling and executing PinsubeditorWalt Fair, Jr.14-Jan-13 10:18 
GeneralRe: Compiling and executing PinmemberDerekEccles15-Jan-13 8:46 
QuestionLaplace transform PinmemberR K Srivastava13-Jun-12 23:23 
AnswerRe: Laplace transform PinsubeditorWalt Fair, Jr.14-Jun-12 3:58 
QuestionThanks a lot Pinmembermeiriweixin22-Jul-11 6:25 
AnswerRe: Thanks a lot PinsubeditorWalt Fair, Jr.23-Jul-11 3:11 
GeneralMy vote of 3 Pinmembermounte@cognimatics26-Jan-11 21:07 
GeneralMy vote of 5 PinmemberLT Mason5-Dec-10 22:57 
GeneralMy vote of 1 Pinmembermikerobe6-Aug-09 13:24 
GeneralRe: My vote of 1 PinsubeditorIndivara14-Dec-10 3:17 
GeneralRe: My vote of 1 Pinmembermounte@cognimatics26-Jan-11 7:03 
GeneralRe: My vote of 1 PinsubeditorIndivara26-Jan-11 14:12 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

| Advertise | Privacy | Mobile
Web02 | 2.8.141022.2 | Last Updated 12 Apr 2008
Article Copyright 2008 by Walt Fair, Jr.
Everything else Copyright © CodeProject, 1999-2014
Terms of Service
Layout: fixed | fluid