Some people may have been wondering about my increasingly frequent absences on the Internet. I am writing this in an attempt to assuage them.

## Quick-and-dirty numerical inversion of the Laplace transform II

February 1, 2017In this post, I will present an extension of the method previously posted in Quick-and-dirty numerical inversion of the Laplace transform.

The method, due to Piessens, is effectively a Gauss-Kronrod extension of the complex Gaussian quadrature method by Salzer. Recall that the Gauss-Kronrod idea is to use two different quadrature abscissa/weight pairs, where one’s abscissas is a subset of the other, to be able to estimate the accuracy of the higher-order method.

Here is the implementation of the extended version of `NInverseLaplaceTransform[]`

:

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

piessensQ[n_Integer, x_] := Module[{cof}, cof = LinearSolve[SparseArray[{{j_, k_} /; j <= k :> 1/((k - j)! Binomial[2 n + k - j, n + k - j])}, {n + 1, n + 1}], -Table[1/(j! Binomial[2 n + j, n + j]), {j, n + 1, 1, -1}]]; Expand[FromDigits[Prepend[Reverse[cof], 1], x]]]

SalzerPiessensRuleData[n_Integer, prec_: MachinePrecision] := Module[{c, pk, pl, sg, sk, sl}, sl = \[FormalX] /. NSolve[HypergeometricPFQ[{-n, n}, {}, \[FormalX]], \[FormalX], prec]; pl = \[FormalX] /. NSolve[piessensQ[n, \[FormalX]], \[FormalX], prec]; sk = HypergeometricPFQ[{1 - n, n - 1}, {}, sl]; sg = (-1)^(n + 1) sl ((2 n - 1)/sk)^2/n; c = n! Binomial[2 n, n]; sk = sg - (2 n - 1) sl/(c n sk piessensQ[n, sl]); pk = 1/(c pl HypergeometricPFQ[{-n, n}, {}, pl] (D[piessensQ[n, \[FormalX]], \[FormalX]] /. \[FormalX] -> pl)); {Join[sl, pl], Join[sk, pk]}]

Options[NInverseLaplaceTransform] = {Method -> Automatic, "Points" -> 6, WorkingPrecision -> MachinePrecision}; NInverseLaplaceTransform[f_, s_, t_, opts : OptionsPattern[]] := Module[{met, xa, wt}, met = Switch[OptionValue[Method], "SalzerPiessens" | Automatic, SalzerPiessensRuleData, "Salzer", SalzerRuleData, _, Return[$Failed]]; {xa, wt} = met[OptionValue["Points"], OptionValue[WorkingPrecision]]; Chop[(wt.Function[s, f][1/xa/t])/t]]

For error estimation purposes, one could then generate the approximant corresponding to both Salzer and Salzer-Piessens for the same setting of `"Points"`

and look at their difference.

Using the Bessel function example from the previous post,

f1[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t, "Points" -> 9, WorkingPrecision -> 30]; f2[t_] = NInverseLaplaceTransform[1/Sqrt[1 + s^2], s, t, Method -> "Salzer", "Points" -> 9, WorkingPrecision -> 30];

GraphicsGrid[{{Plot[{BesselJ[0, t], Re[f1[t]]}, {t, 0, 19}, Frame -> True, ImageSize -> Medium, PlotLegends -> {{TraditionalForm[BesselJ[0, t]], "Salzer-Piessens"}}, PlotStyle -> {AbsoluteThickness[4], Automatic}], SpanFromLeft}, {Plot[Re[f1[t] - f2[t]], {t, 0, 19}, Frame -> True, PlotLabel -> Style["Difference between Salzer-Piessens and Salzer methods", Small], PlotRange -> All], Plot[Re[f1[t] - BesselJ[0, t]], {t, 0, 19}, Frame -> True, PlotLabel -> Style["Difference between Salzer-Piessens and true transform", Small], PlotRange -> All]}}]

## A compact Hilbert curve routine

January 30, 2017Here is a relatively compact *Mathematica* routine for generating the -th iterate of a Hilbert curve in . The algorithm is due to Skilling:

HilbertCurve = Compile[{{x, _Integer}, {b, _Integer}, {n, _Integer}}, Module[{t = BitXor[x, Quotient[x, 2]], p = 2, k, q, r, xx}, xx = Total[Table[BitAnd[Quotient[t, 2^(r (n - 1) + k)], 2^r], {r, b - 1, 0, -1}, {k, n - 1, 0, -1}]]; Do[q = p - 1; Do[t = xx[[k]]; If[BitAnd[t, p] != 0, xx[[1]] = xx[[1]] ~BitXor~ q, t = BitAnd[t ~BitXor~ xx[[1]], q]; xx[[{1, k}]] = BitXor[xx[[{1, k}]], t]], {k, n, 2, -1}]; If[BitAnd[xx[[1]], p] != 0, xx[[1]] = xx[[1]] ~BitXor~ q]; p *= 2, {r, b - 1}]; 2 xx/(p - 1) - 1], RuntimeAttributes -> {Listable}]

I’ve seen a variety of routines for generating the Hilbert curve before, but none were as compact or as general as this one. If you want a version that generates exact rational coordinates, it is easy to take out the requisite parts from the `Compile[]`

function.

Here are the 2D and 3D versions of the curve:

With[{p = 6}, Graphics[{ColorData[97, 1], Line[HilbertCurve[Range[0, 2^(2 p) - 1], p, 2]]}, Frame -> True]]

With[{p = 3}, Graphics3D[Tube[HilbertCurve[Range[0, 2^(3 p) - 1], p, 3]], Axes -> True]]

Here’s a perspective projection of a four-dimensional Hilbert curve:

With[{p = 3, k = 2, f = 4, d = 3/2, r = 1/8}, Graphics3D[{Directive[EdgeForm[], CapForm[None], GrayLevel[1/5], Glow[ColorData["Legacy", "SteelBlue"]], Specularity[ColorData["Legacy", "Gold"], 16]], Tube[(f Delete[#, k])/(d - Extract[#, k]) & /@ HilbertCurve[Range[0, 2^(4 p) - 1], p, 4], r]}, Boxed -> False, Lighting -> "Neutral"]]

(This visualization is admittedly a bit messy-looking; if you have suggestions for making a nice display of the four-dimensional version, tell me about it in the comments!)

## Some terms left unsimplified

January 6, 2017I had been slowly going through some OEIS entries in an attempt to add or optimize the *Mathematica* code found there. One thing that struck me with respect to the sequences that are related to the elliptic integrals/elliptic functions was that they were often left in the form that was spit out by a CAS. I had railed about the deficiencies of current software for generating simple expressions for elliptic integrals many times before, so this was but a confluence with the fact that not many people these days have experience in the analytical manipulation of these functions, as classical and important as they are.

I certainly have my work cut out for me…

## Using the random.org generator in Mathematica

May 15, 2012An answer by the user Emre on the Mathematica StackExchange site introduced me to the random.org random number generation service by Mads Haahr. In the spirit of a previous post of mine, I wanted to be able to use the random.org service through the method plug-in framework for the random number generators of *Mathematica*. I have thus written a barebones set of routines for the purpose:

`RandomOrg /: Random`InitializeGenerator[RandomOrg, ___] := RandomOrg["http://www.random.org/"]`

`RandomOrg[___]["GeneratesRealsQ"] := True;`

RandomOrg[___]["GeneratesIntegersQ"] := True;

`RandomOrg[___]["SeedGenerator"[seed_]] := RandomOrg["http://www.random.org/"]`

`RandomOrg[url_]["GenerateReals"[n_, {a_, b_}, prec_]] := {a + (b - a) Import[url <> "decimal-fractions/?col=1&format=plain&num=" <> ToString[n] <> "&dec=" <> ToString[Round[prec]] <> "&rnd=new", "List"], RandomOrg[url]}`

`RandomOrg[url_]["GenerateIntegers"[n_, {a_, b_}]] := {Import[url <> "integers/?col=1&base=10&format=plain" <> "&min=" <> ToString[a] <> "&max=" <> ToString[b] <> "&num=" <> ToString[n] <> "&rnd=new", "List"], RandomOrg[url]}`

Here are a few examples of its use:

`BlockRandom[SeedRandom[0, Method -> RandomOrg];`

RandomReal[{0, 1}, 10, WorkingPrecision -> 20]]

BlockRandom[SeedRandom[0, Method -> RandomOrg];

RandomInteger[{1, 10}, 10]]

As a tiny reminder, the site has a daily quota in place that limits the amount of random numbers you can generate using their service in a single day. See their website for more details.

## A short note on Costa’s minimal surface

May 14, 2012The other day, I finally managed to simplify Alfred Gray’s parametric equations for Costa’s minimal surface. I might edit this post later, with details on how to manipulate the Weierstrass elliptic functions that show up in the equations, but enjoy these for now:

Here, and are respectively the Weierstrass elliptic function, its derivative, and the Weierstrass zeta function, with invariants and , and and are the usual beta and gamma functions. The invariants are the ones corresponding to the semi-periods and . The parameter ranges are and .

I was very delighted at my result, that I decided to make a stylized plot of the surface in *Mathematica*, using Perlin’s simplex noise for the coloring: