inordinatum

Physics and Mathematics of Disordered Systems

Posts Tagged ‘laplace transform

Tricks for inverting a Laplace Transform, part V: Pole Decomposition

with 4 comments

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

1. Result

In this post I will explain how to find the inverse P(x) of the following Laplace transform:

\displaystyle \int_0^\infty \mathrm{d}x\, e^{-sx}P(x) = \frac{\tanh \sqrt{2s}}{\sqrt{2s}}.   (1)

The solution is given in terms of the Jacobi theta function \theta_2 (which can be re-expressed in terms of Jacobi elliptic functions, though I won’t discuss that here in detail):

\displaystyle P(x) = \frac{1}{2}\theta_2 \left(0; e^{-x\pi^2/2}\right).   (2)

By taking derivatives or integrals with respect to s, which are equivalent to multiplication or division of P by x, one can obtain many other Laplace transform identities, including for example

\displaystyle \int_0^\infty \mathrm{d}x\, e^{-sx}\frac{1}{2x}\theta_2 \left(0; e^{-x\pi^2/2}\right) = \ln \cosh \sqrt{2s}.

If anyone manages to find a systematic list of those, I’d be very grateful. But for now let’s just see how one obtains (2).

2. Derivation

First, we use the classic decomposition of trigonometric functions in infinite products due to Euler:

\displaystyle   \begin{array}{rl}  \sinh z &= z \prod_{n=1}^{\infty}\left(1+\frac{z^2}{\pi^2 n^2}\right) \\  \cosh z &= \prod_{n=1}^{\infty}\left(1+\frac{z^2}{\pi^2 \left(n-\frac{1}{2}\right)^2}\right)  \end{array}  .

From the second identity, we can obtain a partial fraction decomposition of \tanh z (following this post on StackExchange):

\displaystyle   \begin{array}{rl}  \tanh z = & \partial_z \ln \cosh z = \partial_z \sum_{n=1}^{\infty} \ln \left(1+\frac{z^2}{\pi^2 \left(n-\frac{1}{2}\right)^2}\right)   \\  = & \sum_{n=1}^{\infty} \frac{\frac{2z}{\pi^2 \left(n-\frac{1}{2}\right)^2}}{1+\frac{z^2}{\pi^2 \left(n-\frac{1}{2}\right)^2}} \\  = & 2z \sum_{n=1}^{\infty} \frac{1}{\frac{\pi^2(2n-1)^2}{4}+z^2}.  \end{array}  .

Applying this to the right-hand side of (1), we obtain a sum over simple poles:

\displaystyle \int_0^\infty \mathrm{d}x\, e^{-sx}P(x) = \frac{\tanh \sqrt{2s}}{\sqrt{2s}} = 2\sum_{n=1}^{\infty} \frac{1}{\frac{\pi^2(2n-1)^2}{4}+2s}

The Laplace inverse of a simple pole is just an exponential, \int_0^\infty \mathrm{d}x\, e^{-sx}\,e^{-p x}=\frac{1}{p+s}. By linearity of the Laplace transform, we can invert each summand individually, and obtain an infinite sum representation for P(x):

\displaystyle P(x) = \sum_{n=1}^{\infty} \exp\left(-\frac{\pi^2(2n-1)^2}{8}x\right)

This sum can now be evaluated with Mathematica‘s Sum command, or by hand using the representation of theta functions in terms of the nome, for argument z=0 \Leftrightarrow w=1 and q=e^{-\pi^2 x/2}. This finally gives the solution as claimed above:

\displaystyle P(x) = \frac{1}{2}\theta_2 \left(0; e^{-x\pi^2/2}\right).   (2)

Have fun!

Advertisements

Written by inordinatum

April 15, 2013 at 10:27 pm

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

Written by inordinatum

November 4, 2012 at 2:14 pm

Tricks for inverting a Laplace Transform, part III

with 6 comments

This is a continuation of the series of articles on Laplace transforms… You can also have a look at part I (guesses based on series expansions), part II (products and convolutions), part IV (substitutions), and part V (pole decomposition).

This time we will cover the following inverse Laplace problems:

\displaystyle \int_0^\infty e^{\lambda x}P(x)\,dx = e^{\frac{1}{a+b \lambda}}

The basic result that can be used for such problems is the Laplace transform of the Bessel function:

\int_0^\infty e^{\lambda x}\left[\frac{1}{\sqrt{x}}I_1(2\sqrt{x})+\delta(x)\right]dx = e^{\frac{1}{\lambda}}

Now, using

\displaystyle LT\{f(x)e^{cx}\}(s)=\int_0^\infty \mathrm{d}x\,e^{s x}f(x)e^{c x} = LT\{f(x)\}(s+c),

and by rescaling \lambda, one finds that \int_0^\infty e^{\lambda x}P(x)\,dx = e^{\frac{1}{a+b \lambda}} is solved by:

\displaystyle P(x)= \left[\frac{1}{\sqrt{bx}}I_1\left(2\sqrt{\frac{x}{b}}\right)+\delta(x)\right]e^{-\frac{a}{b}x}

Have fun!

Written by inordinatum

October 15, 2011 at 11:17 am

Tricks for inverting a Laplace Transform, part II: Products and Convolutions

with 4 comments

EDIT: In the meanwhile, I have continued the series of posts on Laplace Transform inversion. You can find the subsequent articles here: part I (guesses based on series expansions), part III, part IV (substitutions), part V (pole decomposition). Enjoy!

Following the previous post on inverting Laplace transforms, here is another trick up the same alley. This one actually considers a generalization of the previous case

Find P(x) such that \int_0^\infty P(x)e^{\lambda x}dx=\left(1-q\lambda\right)^a \left(1-\lambda\right)^b.

As usual, the built-in InverseLaplaceTransform function from Mathematica 8 fails to give a result. To obtain a closed formula manually, note that each of the factors can be easily inverted:
\int_0^\infty e^{\lambda x}R_{a,q}(x)\,dx=(1-q\lambda)^a
has the solution
R_{a,q}(x)=\frac{e^{-\frac{x}{q}}\left(\frac{x}{q}\right)^{-1-a}}{q \Gamma(-a)}.

Hence, using the fact that Laplace transforms of convolutions give products, the solution for P(x) can be written as a convolution:
P(x) = \int_0^x dx' R_{a,q}(x')R_{b,1}(x-x').

Computing the integral gives the following expression for P(x) in terms of the hypergeometric function \,_1 F_1:

P(x) = \frac{e^{-x}q^a x^{-1-a-b}}{\Gamma(-a-b)} \,_1 F_1\left(-a,-a-b,\frac{q-1}{q}x\right)

Enjoy!!

Written by inordinatum

September 3, 2011 at 10:34 am

Tricks for inverting a Laplace Transform, part I: Guesses based on Series

with 4 comments

EDIT: In the meanwhile, I have continued the series of posts on Laplace Transform inversion. You can find the subsequent articles here: part II (products and convolutions), part III, part IV (substitutions), part V (pole decomposition). Enjoy!

Working with probability distributions one frequently meets the problem of inverting a Laplace transform, which is notoriously difficult. An example that I met recently is the following:

Find P(x) such that \int_0^\infty P(x)e^{\lambda x}dx=\left(\frac{1-q\lambda}{1-\lambda}\right)^v.

The built-in InverseLaplaceTransform function from Mathematica 8 fails to give a result. However, a closed formula for P(x) can be obtained using a trick: Consider integer v, for which the right-hand side becomes a rational function. Take its Laplace inverse, and then extend the result analytically for non-integer v>0.

Since this may be useful for other cases, too, here a detailed account of how one proceeds:

For small integer values of v, one obtains (e.g. using Mathematica):

v=0: P(x)=\delta(x)
v=1: P(x)=q \delta(x)+(1-q)e^{-x}
v=2: P(x)=q^2 \delta(x)+2q(1-q)e^{-x}+q(1-q)^2 x e^{-x}
v=3: P(x)=q^3 \delta(x)+3q^2(1-q)e^{-x}+3q(1-q)^2 x e^{-x}+\frac{1}{2}(1-q)^3 x^2 e^{-x}
v=4: P(x)=q^4 \delta(x)+4q^3(1-q)e^{-x}+6q^2(1-q)^2 x e^{-x}+2q(1-q)^3 x^2 e^{-x}+\frac{1}{6}(1-q)^4 x^3 e^{-x}
v=5: P(x)=q^5 \delta(x)+5q^4(1-q)e^{-x}+10q^3(1-q)^2 x e^{-x}+5q^2(1-q)^3 x^2 e^{-x}+\frac{5}{6}q(1-q)^4 x^3 e^{-x}+\frac{1}{24}(1-q)^5 x^4 e^{-x}

By trying out, one can now guess a general formula for integer v:

P(x)=q^v \delta(x)+e^{-x}\sum_{k=1}^v \frac{1}{(k-1)!}\binom{v}{k}q^{v-k} (1-q)^k x^{k-1}

For example using Mathematica’s Sum command, one can re-write this in terms of a hypergeometric function:

P(x)=q^v \delta(x) + e^{-x}(1-q)q^{-1+v}v \,_1 F_1(1-v,2,-\frac{1-q}{q}x)

This expression now has a well-defined meaning for any v>0, and is the result we looked for. Actually, Mathematica even knows how to take its Laplace transform, which can be used to check its validity.

More tricks for similar problems to come later!

Written by inordinatum

July 28, 2011 at 6:44 pm