inordinatum

Physics and Mathematics of Disordered Systems

Tricks for inverting a Laplace Transform, part IV: Substitutions

with 9 comments

After a few less technical posts recently, I now continue the series of articles on tricks for the inversion of Laplace transforms. You can find the other parts here: part I (guesses based on series expansions), part II (products and convolutions), part III, and part V (pole decomposition).

Today’s trick will be variable substitutions. Let’s say we need to find P(x) so that

\displaystyle \int_0^\infty \mathrm{d}x\, e^{\lambda x}P(x) = \hat{P}(r(\lambda))=e^{a r(\lambda)}\left[1-r(\lambda)\right]^b,   (1)

where r(\lambda) = \frac{1}{2}\left(c-\sqrt{c^2-\lambda}\right), and a,b,c are constants. This example is not as absurd as it may seem; it actually came up recently in my research on a variant of the ABBM model.

There is a general formula permitting one to invert the Laplace transform above, in terms of an integral in the complex plane:

\displaystyle P(x) = \int_\gamma \frac{\mathrm{d} \lambda}{2\pi i} e^{-\lambda x} \hat{P}(r(\lambda)).

This so-called Fourier-Mellin integral runs along a contour \gamma from -i \infty to +i\infty. \gamma can be chosen arbitrarily, as long as the integral converges and all singularities (poles, branch cuts, etc.) of \hat{P} lie on the right of it. Note that my convention for the Laplace variable \lambda is the convention used for generating functions in probability theory, which has just the opposite sign of \lambda of the “analysis” convention for Laplace transforms.

The fact that our Laplace transform \hat{P} is written entirely in terms of the function r(\lambda) suggests applying a variable substitution. Instead of performing the Fourier-Mellin integral over \lambda, we will integrate over r(\lambda). We then have:

\displaystyle P(x) = \int \frac{\mathrm{d}r}{2\pi i} \frac{\mathrm{d}\lambda(r)}{\mathrm{d}r} e^{-\lambda(r) x + a r}\left(1-r\right)^b.  (2)

Solving the definition of r(\lambda) for \lambda, we get

\displaystyle \lambda(r) = 4 (c r -r^2)

and

\displaystyle \frac{\mathrm{d}\lambda}{\mathrm{d}r} = 4 c - 8r.

Inserting this into (2), we have

\displaystyle P(x) = 4\int_{\gamma'} \frac{\mathrm{d}r}{2\pi i} (c-2r)e^{4 r^2 x - (4cx-a)r}\left(1-r\right)^b.

The only singularity of the integrand is at r=1, for b \in \mathbb{R} \backslash \mathbb{N}. We can thus choose the contour \gamma' to go parallel to the imaginary axis, from 1-i \infty to 1+i\infty. Mathematica then knows to how to evaluate the resulting integral, giving a complicated expression in terms of Kummer’s confluent hypergeometric function \,_1 F_1(a,b,z) = M(a,b,z). However, a much simpler expression is obtained if one introduces an auxiliary integral instead:

\displaystyle g(q,c,d):= \int_{-\infty}^{\infty} \mathrm{d}h\, e^{-q h^2}\left(c-i h\right)^d.

Mathematica knows a simple expression for it:

\displaystyle g(q,c,d) = \sqrt{\pi} q^{-(d+1)/2}U(-d/2;1/2;c^2 q).

U is Tricomi’s confluent hypergeometric function, which is equivalent to Kummer’s confluent hypergeometric function but gives more compact expressions in our case.
Using this auxiliary integral, (2) can be expressed as

\displaystyle P(x) = \left[\frac{4}{\pi}g\left(4x,1+\frac{a-4cx}{8x},b+1\right)+\frac{c-4}{\pi}g\left(4x,1+\frac{a-4cx}{8x},b\right)\right]\exp\left[-\frac{(a-4cx)^2}{16x}\right].

Simplifying the resulting expressions, one obtains our final result for P(x):

\displaystyle   \begin{array}{rcl} P(x) & = & \frac{4}{\sqrt{\pi}}(4x)^{-1-b/2}\left[U\left(-\frac{1+b}{2};\frac{1}{2};\frac{(a+4x(2-c))^2}{16x}\right) + \right. \\  & & \left.+(c-2)\sqrt{x}U\left(-\frac{1+b}{2};\frac{1}{2};\frac{(a+4x(2-c))^2}{16x}\right)\right]\exp\left[-\frac{(a-4cx)^2}{16x}\right].  \end{array} (3)

Equivalently, the confluent hypergeometric functions can be replaced by Hermite polynomials:

\displaystyle    P(x) =  \frac{8}{\sqrt{\pi}}(8x)^{-1-b/2}\left[H_{1+b}\left(\frac{a+4x(2-c)}{4\sqrt{x}}\right) + 2(c-2)\sqrt{x}H_b\left(\frac{a+4x(2-c)}{4\sqrt{x}}\right)\right]\exp\left[-\frac{(a-4cx)^2}{16x}\right].

For complicated Laplace transforms such as these, I find it advisable to check the final result numerically. In the figure below you see a log-linear plot of the original expression (1) for \hat{P}(\lambda), and a numerical evaluation of \int_0^\infty \mathrm{d}x\, e^{\lambda x}P(x), with P(x) given by (3). You can see that they match perfectly!

Numerical verification of the Laplace transform. Yellow line: Original expression \hat{P}(\lambda) as given in (1). Red crosses: Numerical evaluation of \int_0^\infty \mathrm{d}x\,e^{\lambda x}P(x) for the final expression P(x) given in (3).

Advertisements

Written by inordinatum

November 4, 2012 at 2:14 pm

9 Responses

Subscribe to comments with RSS.

  1. […] posts on Laplace Transform inversion. You can find the subsequent articles here: part II, part III, part IV. […]

  2. […] posts on Laplace Transform inversion. You can find the subsequent articles here: part I, part III, part IV. […]

  3. […] This is a continuation of my articles on methods for inverting Laplace transforms. You can find the previous parts here: part I (guesses based on series expansions), part II (products and convolutions), part III, part IV (substitutions). […]

  4. […] have a look at part I (guesses based on series expansions), part II (products and convolutions), part IV (substitutions), and part V (pole […]

  5. Typo: e^(lambda x) –> e^(-lambda x)

    a

    July 6, 2015 at 3:02 pm

  6. Typo?:r(lambda)=1/2 -> r(lambda)=1/(2lambda)

    a

    July 6, 2015 at 3:17 pm

    • You are right… Need to check the +- from Baskara…

      a

      July 6, 2015 at 3:22 pm

      • I’ve re-checked it and I still think that if you start with r(lambda)=1/2(c-\sqrt(c^2-\lambda)) and solve this equation for lambda you will get lambda=4cr-4r^2 as the only solution (no +- ambiguity…)

        inordinatum

        July 7, 2015 at 9:25 pm


Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: