# 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

Posted in Maths, Technical Tricks

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