Quirks in Mathematica’s administrative division data for Mexico

If you ask Mathematica for a list of Mexican states via

    CountryData["Mexico", "RegionNames"]

you will get a list of strings:

    "Aguascalientes", "Baja California", ..., "Zacatecas"}

However, when you try to turn this into a list of objects representing these states via

    states = Entity["AdministrativeDivision", {#, "Mexico"}] & /@ 
                 CountryData["Mexico", "RegionNames"]

something strange happens. Some items in the list are turned into useful objects, and some are uninterpreted symbols.

For example, Aguascalientes is recognized as an administrative division, but Baja California is not. It recognizes Oaxaca but not Nuevo Leon. The pattern is that states with a space are not recognized. There is an inconsistency in Mathematica: output names do not always match input names. To create the object representing Baja California, you need to pass in the string BajaCalifornia with no space.

    Entity["AdministrativeDivision", {"BajaCalifornia", "Mexico"}]

OK, so let’s remove spaces before we try to create a list of geographic objects.

    names = StringReplace[#, " " -> ""] & /@ 
               CountryData["Mexico", "RegionNames"]

This mostly works, but it trips up on Mexico City. The output name for the region is Ciudad de México, but Mathematica does not recognize CiudaddeMéxico as an administrative division. Mathematica does recognize MexicoCity as the name of a city but not as the name of an administrative division.

Changing CiudaddeMéxico to MexicoCity in the list of names did not fix the problem. But when I directly edited the list of state objects by replacing the uninterpreted value with the output running

    Entity["AdministrativeDivision", {"MexicoCity", "Mexico"}]

by itself everything worked. Then I was able to find a Traveling Salesman tour as in earlier posts (Africa, Americas, Eurasia and Oceania, Canada).

Traveling Salesman tour of Mexico

The tour is

  1. Baja California
  2. Baja California Sur
  3. Sinaloa
  4. Durango
  5. Zacatecas
  6. Aguascalientes
  7. Nayarit
  8. Jalisco
  9. Colima
  10. Michoacán
  11. México
  12. Mexico City
  13. Morelos
  14. Guerrero
  15. Oaxaca
  16. Chiapas
  17. Tabasco
  18. Campeche
  19. Quintana Roo
  20. Yucatán
  21. Veracruz
  22. Puebla
  23. Tlaxcala
  24. Hidalgo
  25. Querétaro
  26. Guanajuato
  27. San Luis Potosí
  28. Tamaulipas
  29. Nuevo León
  30. Coahuila
  31. Chihuahua
  32. Sonora

The tour is 8,343 kilometers.

Costas arrays in Mathematica

A couple days ago I wrote about Costas arrays. In a nutshell, a Costas array of size n is a solution to the n rooks problem, with the added constraint that if you added wires between the rooks, no two wires would have the same length and slope. See the earlier post for more details.

The earlier post implemented the Lempel algorithm in Python. Here we’ll implement it in Mathematica. Lempel’s algorithm says to start with parameters p and a where p is a prime and a is a primitive root mod p [1]. Then you can create a Costas array of size n = p-2 by filling in the (i, j) square if and only if

ai + aj = 1 mod p.

We can implement this in Mathematica by

    fill[i_, j_, p_, a_] := If[Mod[a^i + a^j, p] == 1, 1, 0]
    m[p_, a_] := Table[fill[i, j, p, a], {i, 1, p - 2}, {j, 1, p - 2}]

We could visualize the Costas array generated by p = 11 and a = 2 using ArrayPlot. The code

    ArrayPlot[m[11, 2], Mesh -> All]

Generates the following image.

Note that the image is symmetric with respect to the main diagonal. That’s because our algorithm is symmetric in i and j. But Costas arrays are not always symmetrical, so this underscores the point that there is no known scalable algorithm for finding all Costas arrays.

In this example we started with knowing that 2 was a primitive root mod 11. We could have had Mathematica pick a primitive root mod p for us by calling PrimitiveRoot[p]

Just out of curiosity, let’s redo the example above, but instead of testing whether

ai + aj = 1 mod p.

we’ll just use ai + aj mod p.

The code

    ArrayPlot[Table[Mod[2^i + 2^j, 11], {i, 1, 9}, {j, 1, 9}]]

produces the following image.

[1] Lempel’s algorithm generalizes to finite fields of any order, not just integers modulo a prime.

Finding similar world flags with Mathematica

A week ago I posted some pairs of similar flags on Twitter, and later I found that Mathematica’s CountryData database contains flag descriptions. So I thought I’d use the flag descriptions to see which flags Mathematica things are similar.

For example, the FlagDescription attribute for Chad in Mathematica is

Three equal vertical bands of blue (hoist side), yellow, and red; similar to the flag of Romania; also similar to the flags of Andorra and Moldova, both of which have a national coat of arms centered in the yellow band; design was based on the flag of France.

I had Mathematica output a list of countries and flag descriptions, then searched the output for the word “similar.” I then made the following groupings based on the output [1].

Chad / Romania


Bolivia / Ghana

Bolivia   Ghana

Colombia / Ecuador

Equador   Columbia

India / Niger


Ireland / Côte d’Ivoire

Ireland   Ivory Coast

El Salvador / Nicaragua / Honduras

El Salvador    

Egypt / Iraq / Syria / Yemen



Luxembourg / The Netherlands


Andorra / Moldova

Indonesia / Monaco

Indonesia   Monaco


Each flag has an emoji, so here are the groupings above using emoji icons

  • 🇹🇩 🇷🇴
  • 🇧🇴 🇬🇭
  • 🇨🇴 🇪🇨
  • 🇮🇳 🇳🇪
  • 🇮🇪 🇨🇮
  • 🇸🇻 🇳🇮 🇭🇳
  • 🇪🇬 🇮🇶 🇸🇾 🇾🇪
  • 🇱🇺 🇳🇱
  • 🇦🇩 🇲🇩
  • 🇮🇩 🇲🇨

Related posts

[1] The groupings are based on Mathematica’s output, but I did some editing. Strictly following Mathematica’s descriptions would have been complicated. For example, Mathematica’s description might say A is similar to B, but not say B is similar to A. Or it might cluster four flags together that could better be split into two pairs.

Adjacency networks

Suppose you want to color a map with no two bordering regions having the same color. If this is a map on a plane, you can do this using only four colors, but maybe you’d like to use more.

You can reduce the problem to coloring the nodes in a graph. Each node corresponds to a region, and there is an edge between two nodes if and only if their corresponding regions share a border.

Here is a sort of topologists’s or graph theorist’s view of the continental United States.

This was created using the following sample code from the Mathematica documentation.

    RelationGraph[MemberQ[#2["BorderingStates"], #1] &, 
            EntityClass["AdministrativeDivision", "ContinentalUSStates"]]]

You can recognize Maine in the graph because it’s the only state that only borders one other state. Alaska is also easy to locate. Exercise for the reader: mentally add Hawaii to the graph.

The analogous graph for Texas counties took much longer to draw: there are 49 continental US states but 254 Texas counties.

This was created with the following code.

    RelationGraph[MemberQ[#2["BorderingCounties"], #1] &, 
        EntityList[EntityClass["AdministrativeDivision", "USCountiesTexas"]]]

You can find El Paso county in the top left; it only borders one county just as Maine only borders one state.

Related posts

A traveling salesman tour of Africa

Suppose you’d like to tour Africa, visiting each country once, then returning to your starting point, minimizing the distance traveled.

Here’s my first attempt at a solution using Mathematica, based on an example in the documentation for FindShortestTour.

    africa = CountryData["Africa"]
    GeoGraphics[{Thick, Red, GeoPath[africa[[%[[2]]]]]}]

This produced the following map:

Hmm. Maybe I should have been more specific about what I mean by “Africa.” My intention was to find a tour of continental Africa, i.e. not including islands. This means I needed to remove several items from Mathematica’s list of African countries. Also, I had in mind sovereign states, not territories of overseas states and not disputed territories.

After doing this, the map is more like what I’d expect.

The tour is then

  1. Algeria
  2. Tunisia
  3. Libya
  4. Egypt
  5. Chad
  6. Central African Republic
  7. Democratic Republic of the Congo
  8. Burundi
  9. Rwanda
  10. Uganda
  11. South Sudan
  12. Sudan
  13. Eritrea
  14. Djibouti
  15. Somalia
  16. Ethiopia
  17. Kenya
  18. Tanzania
  19. Malawi
  20. Zambia
  21. Mozambique
  22. Zimbabwe
  23. Eswatini
  24. Lesotho
  25. South Africa
  26. Botswana
  27. Namibia
  28. Angola
  29. Republic of the Congo
  30. Gabon
  31. Equatorial Guinea
  32. Cameroon
  33. Nigeria
  34. Niger
  35. Mali
  36. Burkina Faso
  37. Benin
  38. Togo
  39. Ghana
  40. Ivory Coast
  41. Liberia
  42. Sierra Leone
  43. Guinea
  44. Guinea-Bissau
  45. Gambia
  46. Senegal
  47. Mauritania
  48. Morocco

The initial tour, including islands, foreign territories, and Western Sahara, was 23,744 miles or 38,213 kilometers. The second tour was 21,074 miles or 33915 kilometers.

Here’s a tour of just the islands, excluding foreign territories.

The order of the tour is

  1. Cape Verde
  2. Seychelles
  3. Mauritius
  4. Madagascar
  5. Comoros
  6. São Tomé and Príncipe

This tour is 13,034 miles or 20,976 kilometers.

Update: See the next two posts for tours of the Americas and Eurasia and Oceania.

Related posts

Random Blaschke products and Mathematica binding

A Blaschke product is a function that is the product of Blaschke factors, functions of the form

b(z; a) = |a|  (a – z) / a (1 – a*z)

where the complex number a lies inside the unit circle and a* is the complex conjugate of a.

I wanted to plot Blaschke products with random values of a using Mathematica, and I ran into a couple items of interest.

First, although Mathematica has a function RandomComplex for returning complex numbers chosen by a pseudorandom number generator, this function selects complex numbers uniformly over a rectangle and I wanted to select uniformly over a disk. This is easy enough to get around. I wrote my own function by selecting a magnitude and phase at random.

    rc[] := Sqrt[RandomReal[]] Exp[-2 Pi I RandomReal[]]

Now if I reuse the definition of a Blaschke factor from an earlier post

    b[z_, a_] := (Abs[a]/a) (a - z)/(1 - Conjugate[a] z)

I can define a product of two random Blaschke factors as follows.

    f[z_] := b[z, rc[]]  b[z, rc[]]

However, this may not do what you expect. If you plot the function twice, you’ll get different results! It’s a matter of binding order. At what point in the process are the two random values of a chosen and fixed? The answer is some time between the definition of f and executing the plotting function. If the plotting function generated new values of a every time it needed to evaluate f the result would be a hot mess.

On the other hand, if we leave out the colon then the function f behaves as expected.

    f[z_] = b[z, rc[]]  b[z, rc[]]

Now random values are generated by each call to rc, and these values are frozen in the definition of f.

Here’s a Blaschke product with 10 random parameters.

    f[z_] = Product[b[z, rc[]], {i, 1, 10}]

The code

    ComplexPlot[f[z], {z, -1 - I, 1 + I}]

produces the following plot.

If I plot this function again, say changing the plot range, I’ll plot the same function. But if I execute the line of code defining f again, and call ComplexPlot again, I’ll generate a new function.

Incidentally, if I plot this same function using ComplexPlot3d I get the following.

Why does it look like a bowl? As explained in the post on Blaschke factors, each factor is zero at its parameter a and has a pole at the reflection of a in the unit disk.

The function plotted above has zeros inside the unit disk near the boundary. That means it also has poles outside the unit disk near the boundary on the other side.

Inversion in a circle

Inversion in the unit circle is a way of turning the circle inside-out. Everything that was inside the circle goes outside the circle, and everything that was outside the circle comes in.

Not only is the disk turned inside-out, the same thing happens along each ray going out from the origin. Points on that ray that are inside the circle go outside and vice versa. In polar coordinates, the point (r, θ) goes to (1/r, θ).

Complex numbers

In terms of complex numbers, inversion in the unit circle amounts to taking the reciprocal and the conjugate (in either order, because these operations commute). This is the same as dividing a complex number by the square of its magnitude. Proof:

z \bar{z} = |z|^2 \implies \frac{1}{\bar{z}} = \frac{z}{|z|^2}

There are two ways to deal with the case z = 0. One is to exclude it, and the other is to say that it maps to the point at infinity. This can be made rigorous by working on the Riemann sphere rather than the complex plane. More on that here.

Inverting a hyperbola

The other day Alex Kantorovich pointed out on Twitter that “the perfect figure 8 (or infinity) is simply the inversion through a circle of a hyperbola.” We’ll demonstrate this with Mathematica.

What happens to a point on the hyperbola

x^2 - y^2 = 1

when you invert it through a circle? If we think of x and y as the real and imaginary parts of a complex number, the discussion above shows that the point (x, y) goes to the same point divided by its length squared.

(x, y) \mapsto \left( \frac{x}{x^2 + y^2}, \frac{y}{x^2 + y^2} \right)

Let (u, v) be the image of the point (x, y).

\begin{align*} u &= \frac{x}{x^2 + y^2} \\ v &= \frac{y}{x^2 + y^2} \end{align*}

Then we have

u^2 + y^2 = \frac{x^2 + v^2}{(x^2 + y^2)^2} = \frac{1}{x^2 + y^2}


u^2 - v^2 = \frac{x^2 - y^2}{(x^2 + y^2)^2} = \frac{1}{(x^2 + y^2)^2}

because x² – y² = 1. Notice that the latter is the square of the former, i.e.

u^2 - v^2 = (u^2 + v^2)^2

Now we have everything to make our plot. The code

    x^2 + y^2 == 1, 
    x^2 - y^2 == 1, 
    x^2 - y^2 == (x^2 + y^2)^2}, 
    {x, -3, 3}, {y, -3, 3}]

produces the following plot.

The blue circle is the first contour, the orange hyperbola is the second contour, and the green figure eight is the third contour.

Related posts

Catenary kiln

Will Buckner sent me an email with the following question recently. (I’m sharing this with permission.)

I am building a kiln using a catenary arch. The rear wall and front wall/door will be vertical and fill in the space under the arch, which has the dimensions of 41″W x 39.5″H. I need the area within this arch in order to calculate how many bricks I need to construct the two walls. Can you help me calculate that area?

Here’s my solution. I fit a parabola and a catenary in part to check my work and in part to see how different they are.


We’ll start with a parabola as our warmup because it’s easier.

Let w = 41 and h = 39.5. We want a quadratic function y(x) such that y(0) = h and yw/2) = 0.

Assume y has the form

y(x) = ba

y(0) = h leads to the conclusion that b = h, and yw/2) = 0. leads to

a(w/2) = h

And this leads to the final form.

y = h(1 – (2x/w)²)


Our catenary will have the form

y(x) = ba cosh(x/a)

and again y(0) = h and yw/2) = 0. These two requirements lead to the equations

ba = h


b – a cosh(w/2a) = 0.


b = h + a


h + a (1 – cosh(w/2a) = 0.

The latter equation cannot be solved in closed form, but it’s easy enough to solve numerically. When I asked Mathematica to compute the root with

    FindRoot[39.5 + a - a Cosh[20.5/a] == 0, a]

it failed. Then I plotted the function, saw that the root was around 8.5, then gave Mathematica a hint.

    FindRoot[39.5 + a - a Cosh[20.5/a] == 0, {a, 8.5}]

Then Mathematica came back with

a = 8.475359432102444.

So the equation of our catenary is

y = h + a(1 – cosh(0.5 x/a))

where h and a are given above.


Here’s a plot of both curves. The blue curve inside is the parabola, and the gold curve on the outside is the catenary.

Here’s the code that made the plot.

    Plot[{h (1 - (2 x/w)^2), h + a - a Cosh[x/a] }, {x, -w/2, w/2}]


And now for the area, the thing I was asked to find. For the parabola

    Integrate[h (1 - (2 x/w)^2), {x, -w/2, w/2}]

returns 1079.67. For the catenary,

    Integrate[h + a (1 -  Cosh[x/a]), {x, -w/2, w/2}]

returns 1166.56.

The number of bricks this will take is roughly the area of the arch divided by the area of a brick, but being more precise is a complicated question that depends on how much mortar you use, and how you use plan to use rectangular bricks to make a curved arch.

More catenary posts

A closer look at zero spacing

A few days ago I wrote about Kneser’s theorem. This theorem tells whether the differential equation

u ″(x) + h(x) u(x) = 0

will oscillate indefinitely, i.e. whether it will have an infinite number of zeros.

Today’s post will look at another theorem that gives more specific information about the spacing of the zeros.

Kneser’s theorem said that the growth rate of h determines the oscillatory behavior of the solution u. Specifically, if h grows faster than ¼x-2 then oscillations will continue, but otherwise they will eventually stop.

Zero spacing

A special case of a theorem given in [1] says that an upper bound on h gives a lower bound on zero spacing, and a lower bound on h gives an upper bound on zero spacing.

Specifically, if h is bound above by M, then the spacing between zeros is no more than π/√M. And if h is bound below by M′, then the spacing between zeros is no less than π/√M′.

Let’s look back a plot from the post on Kneser’s theorem:

The spacing between the oscillations appears to be increasing linearly. We could have predicted that from the theorem in this post. The coefficient of the linear term is roughly on the order of 1/x² and so the spacing between the zeros is going to be on the order of x.

Specifically, in the earlier post we looked at

h(x) = (1 + x/10)p

where p was 1.9 and 2.1, two exponents chosen to be close to either side of the boundary in Kneser’s theorem. That post used q and t rather than h and x, but I changed notation here to match [1].

Suppose we’ve just seen a zero at x*. Then from that point on h is bounded above by h(x*) and so the distance to the next zero is bounded below by π/√h(x*). That tells you where to start looking for the next zero.

Our function h is bounded below only by 0, and so we can’t apply the theorem above globally. But we could look some finite distance ahead, and thus get a positive lower bound, and that would lead to an upper bound on the location of the next zero. So we could use theory to find an interval to search for our next zero.

More general theorem

The form of the theorem I quoted from [1] was a simplification. The full theorem considers the differential equation

u ″(x) + g(xu′(x) +  h(x) u(x) = 0

The more general theorem looks at upper and lower bounds on

h(x) – ½ g′(x) – ¼ g(x)².

but in our case g = 0 and so the hypotheses reduced to bounds on just h.


Exercise 3 from [1] says to look at the zeros of solution to

u ″(x) + ½x²  u ′ + (10 + 2x) u(x) = 0.

Here we have

h(x) – ½ g ′(x)- ¼ g(x)² = 10 + 2xx – 1  = 9 + x

and so the lower bound of 9 tells us zeros are spaced at most π/3 apart.

Let’s look at a plot of the solution using Mathematica.

    s = NDSolve[{u''[x] + 0.5 x^2 u'[x] + (10 + 2 x) u[x] == 0, 
                 u[0] == 1, u'[0] == 1 }, u, {x, 0, 5}]
    Plot[Evaluate[u[x] /. s], {x, 0, 5}]

This produces the following.

There are clearly zeros between 0 and 1, between 1 and 2, and between 2 and 3. It’s hard to see whether there are any more. Let’s see exactly where the roots are that we’re sure of.

    FindRoot[u[x] /. s, {x, 0, 1}]
    {x -> 0.584153}

    FindRoot[u[x] /. s, {x, 1, 2}]
    {x -> 1.51114}

    FindRoot[u[x] /. s, {x, 2, 3}]
    {x -> 2.41892}

The distance between the first two zeros is 0.926983 and the distance between the second and third zeros is 0.90778. This is consistent with our theorem because π/3 > 1 and the spacing between our zeros is less than 1.

Are there more zeros we can’t see in the plot above? When I asked Mathematica to evaluate

    FindRoot[u[x] /. s, {x, 2, 4}]

it failed with several warning messages. But we know from our theorem that if there is another zero, it has to be less than a distance of π/3 from the last one we found. We can use this information to narrow down where we want Mathematica look.

    FindRoot[u[x] /. s, {x, 2.41892, 2.41892 + Pi/3}]

Now Mathematica says there is indeed another solution.

    {x -> 3.42523}

Is there yet another solution? If so, it can’t be more than a distance π/3 from the one we just found.

    FindRoot[u[x] /. s, {x, 3.42523, 3.42523 + Pi/3}]

This returns 3.42523 again, so Mathematica didn’t find another zero. Assuming Mathematica is correct, this means there cannot be any more zeros.

On the interval [0, 5] the quantity in the theorem is bound above by 14, so an additional zero would need to be at least π/√14 away from the latest one. So if there’s another zero, it’s in the interval [4.26486, 4.47243]. If we ask Mathematica to find a zero in that interval, it fails with warning messages. We can plot the solution over that interval and see that it’s never zero.

[1] Protter and Weinberger. Maximum Principles in Differential Equations. Springer-Verlag. 1984. Page 46.

Quintic trinomial root

This post looks at an exercise from Special Functions, exercise 6 in Appendix E.

Suppose that xm+1 + axb = 0. Show that

x = \frac{b}{a} - \frac{b^{m+1}}{a^{m+2}} + \frac{2m + 2}{2!} \frac{b^{2m+1}}{a^{2m+3}} - \frac{(3m+2)(3m+3)}{3!} \frac{b^{3m+1}}{a^{3m+4}} + \cdots

Use this formula to find a solution to x5 + 4x + 2 = 0 to four decimal places of accuracy. When m = 0 this series reduces to the geometric series. Write this sum as a hypergeometric series.

There are several things I find interesting about this exercise.

  • It’s a possibly useful result.
  • It’s an application of Lagrange’s inversion formula; that’s where the series comes from.
  • The equation has m+1 roots. What’s special about the root given by the series?
  • The sum is not obviously hypergeometric.
  • What happens when the sum diverges? Is it a useful asymptotic series? Does the analytic continuation work?

I want to address the last two points in this post. Maybe I’ll go back and address the others in the future. To simplify things a little, I will assume m = 4.

The parameters of the hypergeometric series are not apparent from the expression above, nor is it even apparent how many parameters you need. It turns out you need m upper parameters and m-1 lower parameters. Since we’re interested in m = 4, that means we have a hypergeometric function of the form 4F3.

x = \frac{b}{a}\,\, {}_4 F_3\left(\frac{1}{5},\frac{2}{5},\frac{3}{5},\frac{4}{5};\frac{1}{2}, \frac{3}{4}, \frac{5}{4} ; -\frac{3125b^4}{256a^5}\right)

We can evaluate this expression in Mathematica as

    f[a_, b_] := (b/a) HypergeometricPFQ[
        {1/5, 2/5, 3/5, 4/5}, 
        {1/2, 3/4, 5/4}, 
        -3125 b^4 / (256 a^5)

When we evaluate f[4, -2] we get -0.492739, which is the only real root of

x5 + 4x + 2 = 0.

Recall that the sign on b is negative, so we call our function with b = -2.

Now lets, try another example, solving for a root of

x5 + 3x – 4 = 0.

If we plug a = 3 and b = 4 into the series at the top of the post, the series diverges immediately, not giving a useful asymptotic series before it diverges.

The series defining our hypergeometric function diverges for |z| > 1, and we’re evaluating it at z = -3125/243. However, the function can be extended beyond the unit disk by analytic continuation, and when we ask Mathematica to numerically evaluate the root with by running

    N{ f[3, 4] ]

we get 1, which is clearly a root of x5 + 3x – 4.

Related posts