3D bifurcation diagram

Posted on 31 July 2021 by John

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)
y ←x + 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.

Related posts: More on dynamical systems

Where has all the productivity gone?

Posted on 31 July 2021 by John

Balaji Srinivasan asks in a Twitter thread why we’re not far more productive given the technology available. Here I collect the five possible explanations he mentions.

The Great Distraction.
All the productivity we gained has been frittered away on equal-and-opposite distractions like social media, games, etc.
The Great Dissipation.
The productivity has been dissipated on things like forms, compliance, process, etc.
The Great Divergence.
The productivity is here, it’s just only harnessed by the indistractable few.
The Great Dilemma.
The productivity has been burned in bizarre ways that require line-by-line “profiling” of everything.
The Great Dumbness.
The productivity is here, we’ve just made dumb decisions in the West while others have harnessed it.

If I had to choose one of the five, I’d lean toward The Great Dissipation, inventing new tasks to absorb new capacity. This is what happened with the introduction of household appliances. Instead of spending less time doing laundry, for example, we do laundry more often.

Maybe we’re seeing that technological bottlenecks were not as important as we thought.

For example, it’s easier to write a novel using Microsoft Word than using a manual typewriter, but not that much easier. MS Word makes the physical work easier, but most of the effort is mental. (And while moving from Smith Corona 1950 to Word 95 is a big improvement, moving from Word 95 to Word 365 isn’t.)

Technology calls our bluff. Improvements in technology show us that technology wasn’t the obstacle that we thought it was.

Alt tags on tweet images

Posted on 30 July 2021 by John

I learned this morning via a comment that Twitter supports alt text descriptions for images. I didn’t think that it did, and said that it didn’t, but someone kindly corrected me.

When I post equations as images on this site, I always include the LaTeX source code in an alt tag. That way someone using a screen reader can determine the content of the equation. It also helps me if I need to go back and change an equation. I’d like to do the same on Twitter.

Unfortunately, it seems support for this feature is inconsistent. Maybe it’s new. The software I use to manage my Twitter accounts apparently doesn’t offer a way to add alt text.

When I use Twitter via its website, I am able to write alt text but not able to read it. Maybe you have to have accessibility features turned on, which would be unfortunate. People who do not use screen readers occasionally benefit from being able to read photo descriptions. I could imagine, for example, that someone might be curious to see the LaTeX code I used to create an equation image.

I tried a couple experiments, one on my personal account and one on my AlgebraFact account. On the latter, I posted an image of the quadratic formula with the text description

    x = \frac{-b \pm \sqrt{b^2 - 4ac}}{2a}

When I look at the tweet while logged into my AlgebraFact account, I see a little black box in the lower left corner of the image with “ALT” in white letters.

I do not see the same box on an image I posted from my personal account. When I log into my personal account, I no longer see the ALT box on the AlgebraFact tweet but I see one on an image I posted.

Question and request

It appears that with default settings, users cannot see image descriptions. What do you have to do to see the descriptions?

I intend to always put LaTeX source in the alt tag of equation images on this site. If you run across an equation without alt text, please let me know.

Nonlinear phase portrait

Posted on 29 July 2021 by John

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

Posted on 27 July 2021 by John

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.

If n is even, set n ← n /2. Otherwise n ← 3n + 1.
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_] := 
        PolynomialMod[
            If[ (p /. x -> 0) == 0, 
                 p/x, 
                 (x + 1) p + 1
            ], 
            m
        ]

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 x⁷ + x³ + 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 + x²
2 + x + x² + x³
2 x² + 2 x³ + x⁴
2 x + 2 x² + x³
2 + 2 x + x²
x + x³
1 + x²

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

Major memory system telephone keypad

Posted on 26 July 2021 by John

This weekend I had to enter an alphabetic passcode on a numeric keypad. The keypad used the same letter-to-digit convention as a phone, but the letters were not printed on the keypad. That made me think about how much better the Major system is.

I wondered what phone keypads would look like if they used the Major memory system, and so I made the image below.

0: S Z, 1: T D ð θ, 2: N ŋ, 3: M, 4: R, 5: L, 6: ʤ ʧ ʃ ʒ, 7: K G, 8: F V, 9: P B

The Major memory system is a way of encoding numbers as words to make them easier to memorize. The system associates a consonant sound with each digit; you’re free to insert any vowels you like. For example, if you wanted to memorize 745, you might encode it as gorilla.

Note that gorilla decodes to 745 and not 7455 because the word has only one L sound, even though it is spelled with two Ls.

One nice feature of the Major system is that if you multiply a number by 10, you can pluralize its mnemonic. For example, a possible mnemonic for 7450 is gorillas.

The Major system emphasizes sounds because humans remember sounds more easily than symbols. If you have a photographic memory for symbols, just memorize the digits and don’t bother with any mnemonic system.

Some of the sounds associated with digits are not represented by a single letter in English and so the keypad above contains a few IPA (International Phonetic Alphabet) symbols. The number 1 is encoded by any of the sounds “t”, “d”, or “th.” The IPA symbol θ represents the th sound in think and the ð represents the th sound in this. The symbol ŋ as a possible encoding for 2 represents the ng sound in sing. And the number 6 can be encoded as one of several similar sounds: ch, sh, soft g, or soft z.

The conventional phone keypad looks simpler: 2 = A, B, or C, 3 = D, E, or F, etc. It’s the kind of thing James Scott would call “legible,” something that looks simple on paper and warms a bureaucrat’s heart, but doesn’t necessarily work well in practice. The sounds associated with the letters for a given digit have nothing in common, so a number can be represented by dissimilar sounds, and similar sounds do not represent the same number.

Encoding telephone numbers as words is rarely possible using the conventional keypad letters. Phone numbers that do correspond to memorable words are highly valued. Every letter has to correspond to a digit, and it matters how the word is spelled.

The Major system is much more flexible since you’re free to supply vowels as you wish, and you can choose from a wide variety of words that spell a single consonant sound in different ways.

Approximation error over unit disk

Posted on 25 July 2021 by John

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

Posted on 24 July 2021 by John

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.

Upper case, lower case, title case

Posted on 22 July 2021 by John

Converting text to all upper case or all lower case is a fairly common task.

One way to convert text to upper case would be to use the tr utility to replace the letters a through z with the letters A through Z. For example,

    $ echo Now is the time | tr '[a-z]' '[A-Z]'
    NOW IS THE TIME

You could convert to lower case by reversing the arguments to tr.

The approach above works if your text consists of only unadorned Roman letters. But it wouldn’t work, for example, if you gave it a jalapeño or π:

    $ echo jalapeño π | tr '[a-z]' '[A-Z]'
    JALAPEñO π

Using the character classes [:lower:] and [:upper:] won’t help either.

Tussling with Unicode

One alternative would be to use the uc command from the Unicode::Tussle package [1] I mentioned a few days ago. There’s also a lc counterpart, and a tc for title case. These utilities handle far more than Roman letters.

    $ echo jalapeño π | uc
    JALAPEÑO Π

Unicode capitalization rules are a black hole, but we’ll just look at one example and turn around quickly before we cross the event horizon.

Suppose you want to send all the letters in the Greek word σόφος to upper case.

    $ echo σόφος | uc
    ΣΌΦΟΣ

Greek has two lower case forms of sigma: ς at the end of a word and σ everywhere else. But there’s only one upper case sigma, so both get mapped to Σ. This means that if we convert the text to upper case and then to lower case, we won’t end up exactly where we started.

    $ echo σόφος | uc | lc
    σόφοσ

Note that the lc program chose σ as the lower case of Σ and didn’t take into account that it was at the end of a word.

[1] “Tussle” is an acronym for Tom [Christiansen]’s Unicode Scripts So Life is Easier.

Continued fraction from entropy source testing

Posted on 22 July 2021 by John

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$

where

$\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.

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