A little personal update

July 14, 2017

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

Read the rest of this entry »


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

February 1, 2017

In 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"],
               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]],
                    PlotStyle -> {AbsoluteThickness[4], Automatic}],
              {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, 2017

Here is a relatively compact Mathematica routine for generating the b-th iterate of a Hilbert curve in [0,1]^n. 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, 2017

I 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…

Yes, I’m still here…

December 21, 2014

…and to those guys who showed concern and worry during my long-ish absence, I’m fine. Two important things, though: the hard drive containing a lot of my unpublished research crashed, including a good amount that I have not been able to back up; and, I have been quite busy with a (relatively) new job over the past few months, details of which I cannot discuss here.

In short, I’m fine, and I hope you guys will still be here when I’m actually back. Happy Holidays!

~ Jan

Using the random.org generator in Mathematica

May 15, 2012

An 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, 2012

The 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:

\displaystyle x=\frac{\pi u}{2}-\Re\left(\frac{\zeta(u+iv;g_2,0)}{2}-\frac{\pi\,\wp(u+iv;g_2,0)}{\wp^\prime(u+iv;g_2,0)}\right) \\ y=\frac{\pi v}{2}+\Im\left(\frac{\zeta(u+iv;g_2,0)}{2}+\frac{\pi\,\wp(u+iv;g_2,0)}{\wp^\prime(u+iv;g_2,0)}\right) \\ z=\sqrt{\frac{\pi}{8}}\log\left|\frac1{\frac12+\frac{\wp(u+iv;g_2,0)}{\sqrt{g_2}}}-1\right|

Here, \wp, \wp^\prime and \zeta are respectively the Weierstrass elliptic function, its derivative, and the Weierstrass zeta function, with invariants g_2=\left(\frac12\mathrm{B}\left(\frac14,\frac14\right)\right)^4=\frac{\Gamma(1/4)^8}{16\pi^2} and 0, and \mathrm{B}(x,y) and \Gamma(x) are the usual beta and gamma functions. The invariants are the ones corresponding to the semi-periods \omega_1=\frac12 and \omega_3=\frac{i}{2}. The parameter ranges are 0 < u < 1 and 0 < v < 1.

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:

Costa's minimal surface