Quick-and-dirty numerical inversion of the Laplace transform

May 4, 2012

I present here a set of Mathematica routines implementing a quick-and-dirty method for numerically evaluating the inverse Laplace transform. As it is well-known, the problem is a numerically ill-conditioned one, and thus checking the sensibility of your results is very much important here. The value of the following method, however, compared to more elaborate transform inversion methods, is that if one wants a quick estimate for arguments not too far from the imaginary axis, the approximation derived from this method does not take too much effort to evaluate.

The method, due to Herbert Salzer, is a complex-valued n-point Gaussian quadrature, where the nodes are the (complex!) roots of (a special case of) the Bessel polynomial. Due to this, if you want to use this method, the function you wish to transform should be able to evaluate complex arguments.

The method below does not give an error estimate; what you can try for verifying the accuracy of your results is to evaluate for an initial number of points (the default is ten points in the routine below), and then increase the number of points taken, and then compare the results you have. (You might also wish to use a sufficiently high value of the option WorkingPrecision.)

Here are the the Mathematica routines:

SalzerRuleData[n_Integer, prec_:MachinePrecision] :=
Block[{x}, Transpose[Map[{#, (-1)^(n + 1) # ((2 n - 1)/HypergeometricPFQ[{1 - n, n - 1}, {}, #])^2/n} &,
x /. NSolve[HypergeometricPFQ[{-n, n}, {}, x], x, prec]]]]

Options[NInverseLaplaceTransform] = {Points -> 10,
WorkingPrecision -> MachinePrecision};
NInverseLaplaceTransform[f_, s_, t_, opts:OptionsPattern[]] := Module[{xa, wt},
{xa, wt} = SalzerRuleData[OptionValue[Points], OptionValue[WorkingPrecision]];
Chop[(wt.Function[s, f][1/xa/t])/t]]

As an example, we try the inverse Laplace transform of the function \frac1{(1+s)^2}, which is t\exp(-t):

f[t_] = NInverseLaplaceTransform[(1 + s)^(-2), s, t];


With[{t = N[2]}, {t Exp[-t], f[t], t Exp[-t] - f[t]}]
{0.270671, 0.270671 + 5.68434*10^-14 I, 2.91759*10^-8 - 5.68434*10^-14 I}

(The tiny imaginary part is an artifact of the evaluation; you can use Chop[] if you wish.)

As another example, here is a comparison of the approximate inverse transform of \frac1{\sqrt{1+s^2}} with the true inverse transform, J_0(t) (where J_0(t) is the Bessel function of the first kind):

f[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t, WorkingPrecision -> 25];

GraphicsRow[{Plot[{f[t], BesselJ[0, t]}, {t, 0, 15}, Frame -> True],
Plot[Re[f[t] - BesselJ[0, t]], {t, 0, 15}, Frame -> True, PlotRange -> All]}]

Bessel function via inverse Laplace transformation

The plot on the left shows that the true function and the approximation are almost indistinguishable within the range plotted. The plot on the right shows the difference between the approximate and the true function.


Bunnies, teapots, and noise

March 16, 2012

I should be able to have regular Internet access soon. For now, I’d like to share some of the experiments I did with coloring 3D objects with the help of Ken Perlin’s simplex noise in Mathematica. (Relatedly, I’ve found that the implementation of improved Perlin noise supplied in the Mathematica help file has a few subtle flaws, which I was able to fix. That would be for another blog entry, though.)

Here are some of the results of coloring the Utah teapot with simplex noise:

Utah Teapot with simplex noise coloring

Here are some results of coloring the Stanford bunny with simplex noise:

Stanford Bunny with simplex noise coloring

I never would have thought that Mathematica would ever be able to do graphics as fancy as this. I’ll only note that none of these used the Texture[] construct; both models in fact used simplex noise as solid textures (i.e. an RGBColor[] associated with each point in the model).

I might write a (series of) more descriptive post(s) on these things much later; for now, these pictures, entirely generated in Mathematica (no post-processing at all!) will have to suffice.


On making simulated landscapes

March 1, 2012

This entry is based on an idea by Paul Bourke for generating simulated landscapes through the application of the two-dimensional discrete Fourier transform to an array of uniformly-distributed random numbers. I was also inspired by some of the color gradients featured in GRASS, and decided to see if I could make a somewhat realistic color gradient for depicting a simulated landscape. Here is the color gradient that resulted from my many experiments:

topoFake[u_?NumericQ] :=
Blend[{{0, RGBColor[0, 1/10, 1/4]}, {1/5, RGBColor[0, 3/5, 3/5]}, {1/5, RGBColor[206/255, 197/255, 12/17]}, {1/4, RGBColor[39/85, 9/26, 5/44]}, {11/40, RGBColor[14/25, 27/40, 2/9]}, {23/80, RGBColor[9/14, 19/25, 12/85]}, {1/3, RGBColor[7/15, 29/42, 1/255]}, {3/8, RGBColor[3/10, 3/5, 5/98]}, {17/40, RGBColor[0, 38/85, 0]}, {19/40, RGBColor[0, 6/17, 0]}, {21/40, RGBColor[0, 3/10, 0]}, {23/40, RGBColor[0, 1/4, 0]}, {5/8, RGBColor[0, 1/8, 0]}, {5/8, RGBColor[39/85, 9/26, 5/44]}, {43/60, RGBColor[3/5, 25/54, 1/6]}, {97/120, RGBColor[19/23, 7/10, 1/3]}, {151/180, RGBColor[7/8, 3/4, 13/36]}, {313/360, RGBColor[39/44, 15/17, 2/5]}, {9/10, RGBColor[13/17, 9/11, 16/51]}, {9/10, RGBColor[1/2, 1/2, 1/2]}, {91/100, RGBColor[1/2, 1/2, 1/2]}, {91/100, RGBColor[1, 50/51, 50/51]}, {19/20, RGBColor[1, 50/51, 50/51]}, {19/20, RGBColor[3/4, 3/4, 3/4]}, {24/25, RGBColor[3/4, 3/4, 3/4]}, {24/25, RGBColor[1, 50/51, 50/51]}, {1, RGBColor[1, 50/51, 50/51]}}, u] /; 0

Color gradient for topoFake[]

I wanted to try it out on a simpler function first, and decided to use the peaks() function of MATLAB as a test case:

Plot3D[3 (1 - x)^2 E^(-x^2 - (y + 1)^2) - 10 (x/5 - x^3 - y^5) E^(-x^2 - y^2) - 1/3 E^(-(x + 1)^2 - y^2), {x, -3, 3}, {y, -3, 3},
Axes -> None, BoundaryStyle -> None, Boxed -> False,
ColorFunction -> (topoFake[#3] &), Mesh -> False,
PerformanceGoal -> "Quality", PlotPoints -> 85, PlotRange -> All,
ViewPoint -> {-1.3, -2.4, 2.}]

peaks() function, colored with topoFake[]

Not too shabby, eh? Here’s what results if we try out the gradient on a modification of Bourke’s idea:

With[{n = 256},
BlockRandom[SeedRandom[20, Method -> "Legacy"];
ListPlot3D[Abs[InverseFourier[
Fourier[RandomReal[{0, 1}, {n, n}]] * Table[
Sech[(i/n - 0.5)/0.025] * Sech[(j/n - 0.5)/0.025],
{i, n}, {j, n}]]],
Axes -> None, Boxed -> False, BoundaryStyle -> None,
BoxRatios -> {1, 1, 1/10}, ColorFunction -> (topoFake[#3] &),
Mesh -> False, PerformanceGoal -> "Quality"]]]

topoFake[]-colored landscape

The resulting archipelago looks a tad chaotic, but it’s a start. I’m still looking for a way to modify Bourke’s technique to produce bigger “islands”.


Faking marbles

February 29, 2012

Perlin noise embedded on a sphere

Still on the matter of Perlin noise, I finally managed to figure out how to wrap Perlin noise on a sphere (that is, as a Texture[]). It was a bit inconvenient that it seems (at least, to me) that not one of the usual references on Perlin noise mentioned the explicit formulae needed for a proper spherical embedding of Perlin noise. Making use of the functions already defined in the doc page for Compile[], here are two Mathematica functions you can use for generating textures that can be wrapped on a sphere:

sphericalPerlin1 = With[{octaves = 4},
Compile[{{theta, _Real}, {phi, _Real}, {amplitude, _Real},
{frequency, _Real}, {gain, _Real}, {lacunarity, _Real}},
Module[{noiseVal = 0.0, x, y, z, sf, freq = frequency,
amp = amplitude},
sf = Sin[phi];
x = sf*Cos[theta] + 1; y = sf*Sin[theta] + 1; z = Cos[phi] + 1;
Do[
noiseVal += amp*signedNoise[x*freq, y*freq, z*freq];
freq *= lacunarity; amp *= gain,
{octaves}
];
noiseVal
], CompilationOptions -> {"InlineExternalDefinitions" -> True}
]]

sphericalPerlin2 = With[{octaves = 4},
Compile[{{theta, _Real}, {phi, _Real}, {amplitude, _Real},
{frequency, _Real}, {gain, _Real}, {lacunarity, _Real}},
Module[{noiseVal = 0.0, x, y, z, sf, freq = frequency,
amp = amplitude},
sf = Sin[phi];
x = sf*Cos[theta] + 1; y = sf*Sin[theta] + 1; z = Cos[phi] + 1;
Do[
noiseVal += amp*Abs[signedNoise[x*freq, y*freq, z*freq]];
freq *= lacunarity; amp *= gain,
{octaves}
];
noiseVal
], CompilationOptions -> {"InlineExternalDefinitions" -> True}
]]

As an example, here is how to use sphericalPerlin2[] along with DensityPlot[] to generate a texture:

DensityPlot[sphericalPerlin2[th, ph, 0.5, 2.0, 0.5, 2.0],
{th, 0, 2 Pi}, {ph, 0, Pi},
AspectRatio -> Automatic, Frame -> None, PerformanceGoal -> "Quality", PlotPoints -> 200, PlotRange -> All]

and here is the result, after adding appropriate coloring (with thanks to Rob Collyer for the color gradient):

fiery marble map

In fact, this is the texture I used for upper leftmost “marble” in the first image.


On faking a speckled paint finish

February 15, 2012

Perlin color swatches

I’ve been experimenting a lot with Perlin noise recently in Mathematica, thanks to the example implementation given in the documentation for Compile[].

(I must however note that the method implemented in the Mathematica help file is not actually the “classic” Perlin noise. In particular, the compiled fade[] routine is the Hermite quintic produced by InterpolatingPolynomial[{{{0}, 0, 0, 0}, {{1/2}, 1/2}, {{1}, 1, 0, 0}}, x], not the classical one that makes use of InterpolatingPolynomial[{{{0}, 0, 0}, {{1/2}, 1/2}, {{1}, 1, 0}}, x]. The Mathematica code also makes use of the centers of the edges of a cube as the gradients, stored in the grad variable within the function dot[], as proposed in Perlin’s update.)

The grid of swatches depicted above is just one sample of Perlin noise, colored differently with 49 of the gradient color schemes available in ColorData["Gradients"], and rendered with ReliefPlot[]. The choice of coloring scheme makes (at least to me) a heap of difference in the visual effect of Perlin noise. I like how the outputs look like samples of the fancy speckled paint I sometimes see in hardware stores.

I still haven’t fully digested the theory behind Perlin noise (though Gustavson’s notes seem to be one of the better treatments I’ve seen), but my experiments with the routines thus far have been quite encouraging.


On emulating the Texas Instruments random number generator

February 12, 2012

I recently needed to emulate the random number generator internally used by Texas Instruments calculators in Mathematica. Thanks to this United-TI forum post, I found that the algorithm being used is a combined multiplicative linear congruential generator due to Pierre L’Ecuyer. The particular generator used is listed on page 747 of the article; as noted there, the combined MLCG has a period of 2305842648436451838\approx 2.30584\times10^{18}.

I figured this might be useful to other people, so I’m posting the routines needed to get RandomReal[] and related functions to emulate the rand() function (and ilk) on TI calculators:

TexasInstruments /: Random`InitializeGenerator[TexasInstruments, ___] := TexasInstruments[12345, 67890]

TexasInstruments[___]["GeneratesRealsQ"] := True;
TexasInstruments[___]["GeneratesIntegersQ"] := True;

TexasInstruments[___]["SeedGenerator"[seed_]] := If[seed == 0, TexasInstruments[12345, 67890],
TexasInstruments[Mod[40014 seed, 2147483563],
Mod[seed, 2147483399]]];

TexasInstruments[s1_, s2_]["GenerateReals"[n_, {a_, b_}, prec_]] :=
Module[{p = s1, q = s2, temp},
{a + (b - a) Table[p = Mod[40014 p, 2147483563]; q = Mod[40692 q, 2147483399]; temp = (p - q)/2147483563; If[temp < 0, temp += 1]; temp, {n}], TexasInstruments[p, q]}
]

TexasInstruments[s1_, s2_]["GenerateIntegers"[n_, {a_, b_}]] :=
Module[{p = s1, q = s2, temp},
{a + Floor[(b - a + 1) Table[p = Mod[40014 p, 2147483563]; q = Mod[40692 q, 2147483399]; temp = (p - q)/2147483563; If[temp < 0, temp += 1]; temp, {n}]], TexasInstruments[p, q]}
]

This is based on the PRNG method plug-in framework in Mathematica that is described here. For instance, here is the Mathematica equivalent of executing 0→rand:rand(10) (for Zilog Z80 calculators like the TI-83 Plus) or RandSeed 0:rand() (for Motorola 68k calculators like the TI-89):

(* 0→rand:rand(10) *)
BlockRandom[SeedRandom[0, Method -> TexasInstruments];
RandomReal[{0, 1}, 10, WorkingPrecision -> 20]]

Here’s corresponding code for randInt() and randM() on the Z80 calculators:

(* 0->rand:randInt(1,10,10) *)
BlockRandom[SeedRandom[0, Method -> TexasInstruments];
RandomInteger[{1, 10}, 10]]

(* 0->rand:randM(3,3) *)
BlockRandom[SeedRandom[0, Method -> TexasInstruments];
Reverse[Reverse /@ RandomInteger[{-9, 9}, {3, 3}]]]

Due to differences in precision, there is bound to be some discrepancy in the least significant figures of the results from the calculators and the results from the Mathematica emulation; still, I believe that this is a serviceable fake.

You can download a Mathematica notebook containing these definitions and tests here.


Fly me to the moon

February 6, 2012

I had always wanted a nice animation depicting periodic solutions of the restricted three-body problem, and was always a skosh underwhelmed by the animations I’ve found. I then decided, on a quiet afternoon, to see if the new (to me) capabilities of Mathematica were up to the task, and I was not at all disappointed:

visualizing the Arenstorf orbit(click on the image to see a bigger version of the animation)

This is one of my favorite Arenstorf orbits (see also these slides); I went with the route taken in Hairer/Nørsett/Wanner and used the Bulirsch-Stoer extrapolation method for numerically solving the underlying differential equation (Method -> "Extrapolation" in NDSolve[]). (Note that for this animation, I used the Arenstorf orbit depicted in the Bulirsch-Stoer paper instead of the slightly more complicated path depicted in Hairer/Nørsett/Wanner; I can be persuaded to produce an animation for that version, though… ;) )

It wasn’t much trouble to depict the Earth, the Moon, and a spaceship. For ease of visualization, I (grudgingly) decided to exaggerate the sizes of the Moon and the spaceship; if I went with the actual scale relative to the Earth, you’d probably have a hard time seeing those. Mathematica‘s Texture[] function, with texture maps obtained from here, as well as ExampleData[{"Geometry3D", "SpaceShuttle"}] for the spaceship, were useful in this regard.

The one thing I did have trouble with is depicting the starry backdrop of outer space in Graphics3D[], so I elected to skip it. Still, the animation looks quite picturesque, no?


Added 2/7/2012:

I finally decided to do the “four loop” version depicted in Hairer/Nørsett/Wanner as well:

four-looped Arenstorf orbit(click on the image to see a bigger version of the animation)