Plotting the Gauss map

A recent post looked at an example from one of Michael Trott’s tomes. This post looks at another example from the same tome.

Trott made a contour plot of the Gauss map

f(x) = \frac{1}{x} - \left\lfloor \frac{1}{x}\right\rfloor

over the complex plane. I copied his code (almost) and reproduced his plot.

        Abs[1/(x + I y) - Floor[1/(x + I y)]], 
        {x, -1.1, 1.1}, {y, -1.1, 1.1}, 
        PlotPoints -> 40, ColorFunction -> Hue, 
        ContourStyle -> {Thickness[0.00001]}

This produced the following image.

The original code set PlotPoints to 400, and it was taking forever. I got impatient and set the parameter to 40. Looking at the image in the book, it seems the plot with the parameter set to 400 is very similar, except there’s a black tangle in the middle rather than a white space.

In 2004 when Trott wrote his book, Mathematica did not have a command for plotting functions of a complex variable directly, and so Trott rolled his own as a function of two real variables.

Here’s a plot of the same function using the relatively new ComplexPlot function.

Here’s the code that produced the plot.

        (1/z - Floor[1/z]), 
        {z, -1.1 - 1.1 I, 1.1 + 1.1 I}, 
        ColorFunction -> Hue, PlotPoints -> 400

Related posts

3D bifurcation diagram

The following 2D bifurcation diagram is famous. You’ve probably seen it elsewhere.

logistic bifurcation diagram

If you have seen it, you probably know that it has something to do with chaos, iterated functions, fractals, and all that. If you’d like to read in more detail about what exactly the plot means, see this post.

I was reading Michael Trott’s Mathematica Guidebook for Numerics and ran across a 3D version of the above diagram. I’d never seen that before. He plots the iterations of the system

x ← a – y(β x + (1 – β) y)
yx + y²/100

The following plot is for β = 1/4 using the code included in the book. (I changed the parameter α in the book to β because the visual similarity between a and α was a little confusing.)

Trott cites four papers regarding this iteration. I looked at a couple of the papers, and they contain similar systems but aren’t quite the same. Maybe his example is a sort of synthesis of the examples he found in the literature.

Nonlinear phase portrait

I was reading through Michael Trott’s Mathematica Guidebook for Programming and ran across the following plot.

I find the image aesthetically interesting. I also find it interesting that the image is the phase portrait of a differential equation whose solution doesn’t look that interesting. That is, the plot of (x(t), x ‘(t)) is much more visually interesting than the plot of x(t).

The differential equation whose phase portrait is plotted above is

x''(t) + \frac{1}{20}x'(t)^3 + \frac{1}{5} x(t) = \frac{1}{3}\cos(et)

with initial position 1 and initial velocity 0. It’s a familiar damped, driven harmonic oscillator, except the equation is nonlinear because the derivative term is cubed.

Here’s a plot of the solution as a function of time.

The solution basically looks like the solution of the linear case, except the solutions are more jagged near the local maxima and minima. In fact, the plot looks so jagged that I wondered whether this was an artifact of needing more plotting points. But adding more points didn’t make much difference. Interestingly, although this plot looks jagged, the phase portrait is smooth.

I produced the phase portrait by copying Trott’s code and making a couple small modifications.

    sol = NDSolve[{(*differential equation*)
        x''[t] + 1/20 x'[t]^3 + 1/5 x[t] == 
        1/3 Cos[E t],(*initial conditions*)x[0] == 1, x'[0] == 0}, 
        x, {t, 0, 360}, MaxSteps -> Infinity]

    ParametricPlot[Evaluate[{x[t], x'[t]} /. sol], {t, 0, 360}, 
        Frame -> True, Axes -> False, PlotPoints -> 3600]

Apparently the syntax of NDSolve has changed slightly since Trott’s book was published in 2004. The argument x in the code above was written x[t] in Trott’s original code and this did not work in the current version of Mathematica. I also simplified the call to ParametricPlot slightly, though the original code would work.

Polynomial analog of the Collatz conjecture

The Collatz conjecture, a.k.a. the 3n + 1 problem, a.k.a. the hailstone conjecture, asks whether the following sequence always terminates.

Start with a positive integer n.

  1. If n is even, set nn /2. Otherwise n ← 3n + 1.
  2. If n = 1, stop. Otherwise go back to step 1.

The Collatz conjecture remains an unsolved problem, though there has been progress toward a proof. Some people, however, are skeptical whether the conjecture is true.

This post will look at a polynomial analog of the Collatz conjecture. Instead of starting with a positive integer, we start with a polynomial with coefficients in the integers mod m.

If the polynomial is divisible by x, then divide by x. Otherwise, multiply by (x + 1) and add 1. Here x is analogous to 2 and (x + 1) is analogous to 3 in the (integer) Collatz conjecture.

Here is Mathematica code to carry out the polynomial Collatz operation.

    c[p_, m_] := 
            If[ (p /. x -> 0) == 0, 
                 (x + 1) p + 1

If m = 2, the process always converges. In fact, it converges in at most n² + 2n steps where n is the degree of the polynomial [1].

Here’s an example starting with the polynomial x7 + x3 + 1.

    Nest[c[#, 2] &, x^7 + x^3 + 1, 15]

This returns 1, and so in this instance 15 iterations are enough.

If m = 3, however, the conjecture is false. In [1] the authors report that the sequence starting with x² + 1 is periodic with period 6.

The following code produces the sequence of values.

    NestList[c[#, 3] &, x^2 + 1, 6]

This returns the sequence

  • 1 + x2
  • 2 + x + x2 + x3
  • 2 x2 + 2 x3 + x4
  • 2 x + 2 x2 + x3
  • 2 + 2 x + x2
  • x + x3
  • 1 + x2

Related posts

[1] Kenneth Hicks, Gary L. Mullen, Joseph L. Yucas and Ryan Zavislak. A Polynomial Analogue of the 3n + 1 Problem. The American Mathematical Monthly, Vol. 115, No. 7. pp. 615-622.

Approximation error over unit disk

The previous post ended by plotting the error in approximating exp(x) with (2+x)/(2 – x):

error in linear and bilinear approximations to exp

The error is much larger on the right end than the left. That made me wonder what the error might look like in the complex plane. Does the slice along real axis exhibit the minimum and maximum error, or are the other places in the unit disk where the error is even smaller or even larger?

Here’s a plot of the error over the unit square.

Here’s the Mathematica code that produced the plot.

    ComplexPlot3D[ Exp[z] - (2 + z)/(2 - z), 
        {z, -1 - I, 1 + I}, PlotRange -> All]

It looks like the error might indeed be minimized and maximized along the real axis.

Incidentally, we can tell the error has a zero of order three at the origin because the sequence of colors, representing the phase, goes through three cycles.

We know the maximum and minimum of the absolute value of the error over the unit disk both occur on the boundary, i.e. on the unit circle, by the maximum modulus theorem. So we plot the error along the rim of the unit disk by letting z = exp(2πit) for 0 ≤ t ≤ 1.

Here’s the code that produced the plot.

    Plot[Abs[f[Exp[2 Pi I t]]], {t, 0, 1}]

It appears the error is maximized when t is 0 or 1, i.e. when z = 1, and minimized when t = 1/2, i.e. when z = -1.

Simple derivation of exponential approximation

I was watching one of Brian Douglas’ videos on control theory (Discrete Control #5) and ran into a simple derivation of an approximation I presented earlier.

Back in April I wrote several post on simple approximations for log, exp, etc. In this post I gave an approximation for the exponential function:

\exp(x) \approx \frac{2 + x}{2 - x}

The control theory video arrives at the same approximation as follows.

\begin{align*} \exp(x) &= \exp(x/2)\, \exp(x/2) \\ &= \frac{\exp(x/2)}{\exp(-x/2)} \\ &\approx \frac{1 + x/2}{1 - x/2} \\ &= \frac{2 + x}{2-x} \end{align*}

As I believe I’ve suggested before here, in a derivation like the one above, where you have mostly equalities and one or two approximations, pay special attention to the approximation steps. The approximation step above uses a first order Taylor approximation in the numerator and denominator.

The plot below shows that the approximation above (the bilinear approximation) is more accurate than doing a single Taylor approximation, approximating exp(x) by 1 + x (linear approximation).

exp, linear approx, bilinear approx

Here’s a plot focusing on the error in the bilinear and linear approximations.

error in linear and bilinear approximations to exp

The bilinear approximation is hard to tell from 0 in the plot above for x up to 0.5.

The derivation above is simple, but why is the result so good? An explanation in terms of Padé approximation is given here.

Continued fraction from entropy source testing

NIST publication 800-90B, Recommendations for the Entropy Sources Used for Random Bit Generation, contains an interesting continued fraction.

Define a function

F\left(\frac{1}{z}\right) = \Gamma(3,z) z^{-3} e^z


\Gamma(a, b) = \int_b^\infty t^{a-1} e^{-t} \, dt

is the incomplete gamma function.

NIST 800-90B gives [1] a continued fraction implement F, but the continued fraction includes a parameter k for no apparent reason. NIST cites a paper that in turn cites Abramowitz and Stegun formula 6.5.31.

It appears that the k in NIST 800-90B corresponds to a-1 in A&S 6.5.31, where a is the first argument to the incomplete gamma function and the exponent on 1/z, and so a = 3.

The continued fraction in A&S is

\cfrac{1}{z+\cfrac{1-a}{1+\cfrac{1}{z+\cfrac{2-a}{1+ \cfrac{2}{z +\ddots}}}}}

However, it’s not entirely clear how to fill in the dots. Apparently the denominators alternate between z and 1, but what about the numerators. The pattern we are give is

1, 1-a, 1, 2-a, 2

That’s as far as A&S goes. My initial guess was that we should ignore the 1 at the beginning of the series. In that case the rest of the series would be

3-a, 3, 4-a, 4, …

and that turns out to be right. If you don’t ignore the first 1, it can be hard to figure out how it fits into the pattern (because it doesn’t!).

How many terms of the continued fraction do we need? It depends on the size of the argument z. For large z, say z > 4, the terms shown above give a good approximation. But for smaller z, the approximation is terrible. For example, F(1) = 5, but the continued fraction evaluated at 1 gives -3/2.

The NIST paper requires evaluating F for values around 1, and it would take some work to figure out how far to carry the continued fraction before it is reliable in that range.

Related posts

[1] Why does NIST include a way to compute F? The incomplete gamma function is difficult to implement well and some mathematical libraries don’t include it. Having a stand-alone implementation of this function might make a testing program more self-contained.

Predator-Prey period

The Lotka-Volterra equations are a system of nonlinear differential equations for modeling a predator-prey ecosystem. After a suitable change of units the equations can be written in the form

\begin{align*} \frac{dx}{dt} &= \phantom{-}ax(y-1) \\ \frac{dy}{dt} &= -by(x-1) \end{align}

where ab = 1. Here x(t) is the population of prey at time t and y(t) is the population of predators. For example, maybe x represents rabbits and y represents foxes, or x represents Eloi and y represents Morlocks.

It is well known that the Lotka-Volterra equations have periodic solutions. It is not as well known that you can compute the period of a solution without having to first solve the system of equations.

This post will show how to compute the period of the system. First we’ll find the period by solving the equations for a few different initial conditions, then we’ll show how to directly compute the period from the system parameters.

Phase plot

Here is a plot of (x(t), y(t)) showing that the solutions are periodic.

Predator-prey phase plots for varying initial conditions

And here’s the Python code that made the plot above.

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.integrate import solve_ivp
    # Lotka-Volterra equations
    def lv(t, z, a, b):
        x, y = z
        return [a*x*(y-1), -b*y*(x-1)]
    begin, end = 0, 20
    t = np.linspace(begin, end, 300)
    styles = ["-", ":", "--"]
    prey_init = [2, 4, 8]
    a = 1.5
    b = 1/a
    for style, p0 in zip(styles, prey_init):
        sol = solve_ivp(lv, [begin, end], [p0, 1], t_eval=t, args=(a, b))
        plt.plot(sol.y[0], sol.y[1], style)
    plt.legend([f"$p_0$ = {p0}" for p0 in prey_init])    

Note that the derivative of x is zero when y = 1. Since our initial condition sets y(0) = 1, we’re specifying the maximum value of x. (If our initial values of x were less than 1, we’d be specifying the minimum of x.)

Time plot

The components x and y have the same period, and that period depends on a and b, and on the initial conditions. The plot below shows how the period increases with x(0), which as we noted above is the maximum value of x.

Predator-prey solutions over time for varying initial conditions

And here’s the code that made the plot.

    fig, ax = plt.subplots(3, 1)
    for i in range(3):
        sol = solve_ivp(lv, [begin, end], [prey_init[i], 1], t_eval=t, args=(a, b))
        ax[i].plot(t, sol.y[0], "-")
        ax[i].plot(t, sol.y[1], "--")    
        ax[i].set_title(f"x(0) = {prey_init[i]}")

Finding the period

In [1] the author develops a way to compute the period of the system as a function of its parameters without solving the differential equations.

First, calculate an invariant h of the system:

h = b(x - \log(x) - 1) + a(y - \log(y) - 1)

Since this is an invariant we can evaluate it anywhere, so we evaluate it at the initial conditions.

Then the period only depends on a and h. (Recall we said we can scale the equations so that ab = 1, so a and b are not independent parameters.)

If h is not too large, we can compute the approximate period using an asymptotic series.

P(h) = 2\pi\left( 1 + \frac{\sigma}{6} h + \frac{\sigma^2}{144} h^2 + \left(\frac{\sigma}{360} - \frac{139\sigma^3}{38880}\right) h^3 + \cdots \right)

where σ = (a + b)/2 = (a² + 1)/2a.

We use this to find the periods for the example above.

    def find_h(a, x0, y0):
        b = 1/a
        return b*(x0 - np.log(x0) - 1) + a*(y0 - np.log(y0) - 1)

    def P(h, a):
        sigma = 0.5*(a + 1/a)
        s = 1 + sigma*h/6 + sigma**2*h**2/144
        return 2*np.pi*s

    print([P(find_h(1.5, p0, 1), 1.5) for p0 in [2, 4, 8]])

This predicts periods of roughly 6.5, 7.5, and 10.5, which is what we see in the plot above.

When h is larger, the period can be calculated by numerically evaluating an integral given in [1].

Related posts

[1] Jörg Waldvogel. The Period in the Volterra-Lotka Predator-Prey Model. SIAM Journal on Numerical Analysis, Dec., 1983, Vol. 20, No. 6, pp. 1264-1272

Error estimates for splines with more boundary conditions

Yesterday I wrote about rates of convergence for natural cubic splines. This brief post reports similar results for more boundary conditions.

As explained in the earlier post, a cubic spline that interpolates a function f at n+1 points satisfies 4n -2 equations in 4n variables. Two more equations are necessary to uniquely determine the interpolating cubic spline.

A natural cubic spline S over [a, b] satisfies

S”(a) = 0
S”(b) = 0.

There are other ways of coming up with two more equations. You might match first derivatives at the end points

S‘(a) = f‘(a)
S‘(b) = f‘(b)

or second derivatives

S”(a) = f”(a)
S”(b) = f”(b).

Or if f is periodic you could require first and second derivatives of the spline to match at a and b.

S‘(a) = S‘(b)
S”(a) = S”(b).

In each case, you get the same order of convergence as we reported earlier for the natural cubic spline, assuming f is 4 times continuously differentiable over [a, b]. The error goes down like O(1/n4).

Source: D. Kershaw. A Note on the Convergence of Interpolatory Cubic Splines. SIAM Journal on Numerical Analysis, Mar., 1971, Vol. 8, No. 1 pp. 67-74.


When we use natural cubic splines to fit exp(x) over [0, 2π] in the earlier post, the error was fairly large. We get better results when we specify the derivatives at the two ends of the interval. Here’s what we get for 11 points:

The error is about 7 times smaller than with natural cubic splines.

And here’s what we get for 21 points:

In this case the error is about 8x smaller than with natural cubic splines.

When we go from 11 to 21 points, specifying the derivatives on each end, the error drops by about a factor of 15, close value of 16 we’d expect in the limit.

Incidentally, to specify the boundary conditions of the spline in Python we need to add an option parameter bc_type in the call to CubicSpline. Otherwise the code is the same as before. In our example

    bc_type=((1, 1), (1, np.exp(2*np.pi)))

We set bc_type equal to a pair of pairs. The first pair gives a boundary condition on the left end and the second gives a boundary condition on the right. The first argument of the inner pair in both cases is 1, because we’re specifying 1st derivatives. Use a 2 to specify 2nd derivatives. The second argument to the pair the value of the derivative: exp(0) = 1 on the left, and exp(2π) on the right.

Random Fourier series

A theorem by Paley and Wiener says that a Fourier series with random coefficients produces Brownian motion on [0, 2π]. Specifically,

W(t) = \frac{Z_0}{\sqrt{2\pi}}\,t + \frac{2}{\sqrt{\pi}} \sum_{n=1}^\infty \frac{Z_n}{n} \sin\left(\frac{nt}{2}\right)

produces Brownian motion on [0, 2π]. Here the Zs are standard normal (Gaussian) random variables.

Here is a plot of 10 instances of this process.

Plot of 10 random Fourier series

Here’s the Python code that produced the plots.

    import matplotlib.pyplot as plt
    import numpy as np
    def W(t, N):
        z = np.random.normal(0, 1, N)
        s = z[0]*t/(2*np.pi)**0.5
        s += sum(np.sin(0.5*n*t)*z[n]/n for n in range(1, N))*2/np.pi**0.5
        return s
    N = 1000
    t = np.linspace(0, 2*np.pi, 100)
    for _ in range(10):
        plt.plot(t, W(t, N))

Note that you must call the function W(t, N) with a vector t of all the time points you want to see. Each call creates a random Fourier series and samples it at the points given in t. If you were to call W with one time point in each call, each call would be sampling a different Fourier series.

The plots look like Brownian motion. Let’s demonstrate that the axioms for Brownian motion are satisfied. In this post I give three axioms as

  1. W(0) = 0.
  2. For s > t ≥ 0. W(s) – W(t) has distribution N(0, st).
  3. For vust ≥ 0. W(s) – W(t) and W(v) – W(u) are independent.

The first axiom is obviously satisfied.

To demonstrate the second axiom, let t = 0.3 and s = 0.5, just to pick two arbitrary points. We expect that if we sample W(s) – W(t) a large number of times, the mean of the samples will be near 0 and the sample variance will be around 0.2. Also, we expect the samples should have a normal distribution, and so a quantile-quantile plot would be nearly a straight 45° line.

To demonstrate the third axiom, let u = 1 and v = √7, two arbitrary points in [0, 2π] larger than the first two points we picked. The correlation between our samples from W(s) – W(t) and our samples from W(v) – W(u) should be small.

First we generate our samples.

    N = 1000
    diff1 = np.zeros(N)
    diff2 = np.zeros(N)
    x = [0.3, 0.5, 1, 7**0.5]
    for n in range(N):
        w = W(x)
        diff1[n] = w[1] - w[0]
        diff2[n] = w[3] - w[2]

Now we test axiom 2.

    from statsmodels.api import qqplot
    from scipy.stats import norm

    print(f"diff1 mean = {diff1.mean()}, var = {diff1.var()}.")
    qqplot(diff1, norm(0, 0.2**0.5), line='45')

When I ran this the mean was 0.0094 and the variance was 0.2017.

Here’s the  Q-Q plot:

quantile-quantile plot to test normal distribution

And we finish by calculating the correlation between diff1 and diff2. This was 0.0226.

Related posts