The Bring sextic of equilateral pentagons​
by Lyle Ramshaw​
May 14, 2022
This is the calculations of the paper “The Bring sextic of equilateral pentagons” in detail.
Saving this notebook from within Mathematica — even without having done any evaluations or edits — writes out, as the file “EquilateralPentagonMap.m”, a Mathematica package that exports the conformal mapping projection of the paper as a routine named PoincareCoordinates. That routine is both specified and implemented in section [A9] below. (The package "EquilateralPentagonMap.m" consists of the cells in section [A9] that are InitializationCells, with gray backgrounds.)
Both this notebook and that package assume that the data file “ApproximatingPolynomials.m” will be available in the directory from which the notebook or package is read.

[A1]: The singular points of
E
4

Figure 3 depicts the polygon space
E
4
of equilateral quadrilaterals. It consists of three loops of decomposables, each of which is a geodesic under the bag-of-edges metric
e
4
. But what does
e
4
say about the three intersection points?
At a singular point of a polygon space, all of the
n
-gon’s edges lie along a single line. Perturbing their phases thus leads to a failure-to-close vector that is always orthogonal to that line, to first order. With only one dimension of failure-to-close vectors, rather than two, the dimension of the tangent space jumps up, at a singular point, from
n-3
to
n-2
; for
E
4
, from a tangent line to a tangent plane. And the metric
e
4
gives that entire tangent plane a Euclidean structure.
Consider the bottom intersection
θ=(0,π,0,π)
, where the four edges of the quadrilateral go back and forth over the same segment twice each way. We can leave
θ
by opening up that quadrilateral into an L that is traveled forward and backward, with the angle of the L at either the left or the right. There are tangent vectors to
E
4
at
θ
that point either way — for example,
∂
left
=(0,0,1,1)
and
∂
right
=(0,-1,-1,0)
. Their sum
∂
left
+
∂
right
=(0,-1,0,1)
is also a tangent vector to
E
4
, lying in the tangent plane to
E
4
at the singular point
θ
. Note, though, that we can’t move more than an infinitesimal distance in the direction
∂
left
+
∂
right
while remaining in
E
4
. Opening up on both sides simultaneously while keeping the edge lengths fixed produces quadrilaterals that close up only to first order, not to second order. To have them close up, looking like
, we’d have to lengthen their crossing edges or shorten their noncrossing edges. Such bowtie-shaped quadrilaterals are modeled in the Grassmannian

G
(2,4)
, since we can balance the lengthening and shortening so as to preserve the required fixed perimeter.
We can determine the angle between the vectors
∂
left
and
∂
right
by applying, to them and to their sum, our formula in (1) for the bag-of-edges metric
p
r
, where
re(4)
1
4
,
1
4
,
1
4
,
1
4
. Taking the six pairs
(i,j)
in (1) in lexicographic order, we have
p
r
(
∂
left
)(0+1+1+1+1+0)/161/4
and
p
r

∂
right
(1+1+0+0+1+1)/161/4
, while
p
r

∂
left
+
∂
right
(1+0+1+1+4+1)/161/2
. It follows that

∂
left
==
∂
right
1/2
, while

∂
left
+
∂
right
1
2
, which implies that
∂
left
and
∂
right
are orthogonal.
Another way to see that this must be so is to reason as follows: For the manifold-with-singularities
(
E
4,
e
4
)
to have 48 isomorphisms, we must be able to map the intersection
θ
to any of the three intersections, say
′
θ
, and to map the tangent plane at
θ
to the tangent plane at
′
θ
in one of eight ways, leaving a binary choice about how to map the two halves of the top loop. In order to have eight different ways to map the tangent plane, all three intersections must be orthogonal.

[A2]: Figures 2 depicts
+
E
5

We here verify that Figure 2 depicts the polygon space
E
5
with the orientation
+
E
5
, the orientation suggested by the argument in Section 7.2. Since Figures 7 and 8 have the same orientation as Figure 2 where they overlap, they also depict
+
E
5
.​ Let’s consider the isosceles trapezoid along the horizontal axis in Figure 2, which has an orange top. If we rotate that trapezoid clockwise (CW) to make its double-length edge horizontal, it has phase vector
(2π/3,π,-2π/3,0,0)
, taking the edges in rainbow order from red to purple.
◼
  • The tangent vector
    (0,1,0,0,1)
    points along the ditri-pentagon loop that runs diagonally up and to the right. As we move in that direction, we rotate the orange and purple edges CCW in lockstep, passing through a child’s house and then leaving the decatract of convex-and-CCW pentagons at the purple-topped isosceles trapezoid.
    • ◼
  • The tangent vector
    (0,-1,0,-1,0)
    points along the ditri-pentagon loop that runs diagonally up and to the left. As we move in that direction, we rotate the orange and blue edges CW in lockstep, passing through a child’s house and leaving the decatract at the blue-topped trapezoid. So the vectors
    (0,1,0,0,1)
    and
    (0,-1,0,-1,0)
    , in that order, are a positive basis for the tangent space to
    E
    5
    , when that tangent space has the orientation of Figure 2.
    • ◼
  • The vector
    (1,1,1,1,1)
    is a positive basis for the real line
    
    of rigid rotations.
    • ◼
  • The orange edge
    ε
    2
    is CCW from the red edge
    ε
    1
    . Hence, rotating
    ε
    2
    a bit CCW gives a failure-to-close vector that is CCW of the analogous vector that would result from rotating
    ε
    1
    a bit CCW. Thus, the vectors
    (1,0,0,0,0)
    and
    (0,1,0,0,0)
    , in that order, are a positive basis for the plane
    
    of failures to close.
    • Figure 2 depicts
    +
    E
    5
    because this 5-by-5 determinant is positive:
    In[]:=
    Det[{{0,1,0,0,1},​​{0,-1,0,-1,0},​​{1,1,1,1,1},​​{1,0,0,0,0},​​{0,1,0,0,0}}]
    Out[]=
    1

    [A3]: Relations among the vertex angles

    We begin with the two equations in (3):
    In[]:=
    cosSum=1+Cos[
    τ
    3
    ]+Cos[
    τ
    3
    ]Cos[
    τ
    4
    ]-Sin[
    τ
    3
    ]Sin[
    τ
    4
    ]+Cos[
    τ
    2
    ]+Cos[
    τ
    2
    ]Cos[
    τ
    1
    ]-Sin[
    τ
    2
    ]Sin[
    τ
    1
    ]
    Out[]=
    1+Cos[
    τ
    2
    ]+Cos[
    τ
    1
    ]Cos[
    τ
    2
    ]+Cos[
    τ
    3
    ]+Cos[
    τ
    3
    ]Cos[
    τ
    4
    ]-Sin[
    τ
    1
    ]Sin[
    τ
    2
    ]-Sin[
    τ
    3
    ]Sin[
    τ
    4
    ]
    In[]:=
    sinSum=Sin[
    τ
    3
    ]+Sin[
    τ
    3
    ]Cos[
    τ
    4
    ]+Cos[
    τ
    3
    ]Sin[
    τ
    4
    ]-Sin[
    τ
    2
    ]-Sin[
    τ
    2
    ]Cos[
    τ
    1
    ]-Cos[
    τ
    2
    ]Sin[
    τ
    1
    ]
    Out[]=
    -Cos[
    τ
    2
    ]Sin[
    τ
    1
    ]-Sin[
    τ
    2
    ]-Cos[
    τ
    1
    ]Sin[
    τ
    2
    ]+Sin[
    τ
    3
    ]+Cos[
    τ
    4
    ]Sin[
    τ
    3
    ]+Cos[
    τ
    3
    ]Sin[
    τ
    4
    ]
    In[]:=
    algSubs=JoinTableCos[
    τ
    k
    ]1-
    2
    t
    k
    1+
    2
    t
    k
    ,{k,4},​​TableSin[
    τ
    k
    ]2
    t
    k
    1+
    2
    t
    k
    ,{k,4}
    Out[]=
    Cos[
    τ
    1
    ]
    1-
    2
    t
    1
    1+
    2
    t
    1
    ,Cos[
    τ
    2
    ]
    1-
    2
    t
    2
    1+
    2
    t
    2
    ,Cos[
    τ
    3
    ]
    1-
    2
    t
    3
    1+
    2
    t
    3
    ,Cos[
    τ
    4
    ]
    1-
    2
    t
    4
    1+
    2
    t
    4
    ,Sin[
    τ
    1
    ]
    2
    t
    1
    1+
    2
    t
    1
    ,Sin[
    τ
    2
    ]
    2
    t
    2
    1+
    2
    t
    2
    ,Sin[
    τ
    3
    ]
    2
    t
    3
    1+
    2
    t
    3
    ,Sin[
    τ
    4
    ]
    2
    t
    4
    1+
    2
    t
    4
    
    In[]:=
    commonDenom=Product[1+
    2
    t
    k
    ,{k,4}]
    Out[]=
    1+
    2
    t
    1
    1+
    2
    t
    2
    1+
    2
    t
    3
    1+
    2
    t
    4
    
    In[]:=
    Z
    cos
    =commonDenomFactor[cosSum/.algSubs]
    Out[]=
    5+3
    2
    t
    1
    -4
    t
    1
    t
    2
    +
    2
    t
    2
    +3
    2
    t
    1
    2
    t
    2
    +
    2
    t
    3
    -
    2
    t
    1
    2
    t
    3
    -4
    t
    1
    t
    2
    2
    t
    3
    -3
    2
    t
    2
    2
    t
    3
    -
    2
    t
    1
    2
    t
    2
    2
    t
    3
    -4
    t
    3
    t
    4
    -4
    2
    t
    1
    t
    3
    t
    4
    -4
    2
    t
    2
    t
    3
    t
    4
    -4
    2
    t
    1
    2
    t
    2
    t
    3
    t
    4
    +3
    2
    t
    4
    +
    2
    t
    1
    2
    t
    4
    -4
    t
    1
    t
    2
    2
    t
    4
    -
    2
    t
    2
    2
    t
    4
    +
    2
    t
    1
    2
    t
    2
    2
    t
    4
    +3
    2
    t
    3
    2
    t
    4
    +
    2
    t
    1
    2
    t
    3
    2
    t
    4
    -4
    t
    1
    t
    2
    2
    t
    3
    2
    t
    4
    -
    2
    t
    2
    2
    t
    3
    2
    t
    4
    +
    2
    t
    1
    2
    t
    2
    2
    t
    3
    2
    t
    4
    In[]:=
    Z
    sin
    =commonDenomFactor[sinSum/.algSubs]
    Out[]=
    -2
    t
    1
    +2
    t
    2
    -
    t
    1
    2
    t
    2
    -2
    t
    3
    -2
    2
    t
    1
    t
    3
    -2
    2
    t
    2
    t
    3
    -2
    2
    t
    1
    2
    t
    2
    t
    3
    +
    t
    1
    2
    t
    3
    +2
    t
    2
    2
    t
    3
    -
    t
    1
    2
    t
    2
    2
    t
    3
    -
    t
    4
    -
    2
    t
    1
    t
    4
    -
    2
    t
    2
    t
    4
    -
    2
    t
    1
    2
    t
    2
    t
    4
    +
    2
    t
    3
    t
    4
    +
    2
    t
    1
    2
    t
    3
    t
    4
    +
    2
    t
    2
    2
    t
    3
    t
    4
    +
    2
    t
    1
    2
    t
    2
    2
    t
    3
    t
    4
    +
    t
    1
    2
    t
    4
    +2
    t
    2
    2
    t
    4
    -
    t
    1
    2
    t
    2
    2
    t
    4
    +
    t
    1
    2
    t
    3
    2
    t
    4
    +2
    t
    2
    2
    t
    3
    2
    t
    4
    -
    t
    1
    2
    t
    2
    2
    t
    3
    2
    t
    4
    
    We first use a Gröbner basis to eliminate
    t
    4
    . The second argument to GroebnerBasis is a list of variables in which the first is the most expensive, that is, the one that the basis should try the hardest to avoid.
    In[]:=
    gbOmit4=GroebnerBasis[{
    Z
    cos
    ,
    Z
    sin
    },{
    t
    4
    ,
    t
    1
    ,
    t
    2
    ,
    t
    3
    }][[1]]//Factor
    Out[]=
    2
    1+
    2
    t
    1
    
    2
    1+
    2
    t
    2
    
    2
    1+
    2
    t
    3
    
    -15-3
    2
    t
    1
    +16
    t
    1
    t
    2
    +
    2
    t
    2
    -3
    2
    t
    1
    2
    t
    2
    +8
    t
    1
    t
    3
    +16
    t
    2
    t
    3
    -8
    t
    1
    2
    t
    2
    t
    3
    -3
    2
    t
    3
    +
    2
    t
    1
    2
    t
    3
    -3
    2
    t
    2
    2
    t
    3
    +
    2
    t
    1
    2
    t
    2
    2
    t
    3
    
    The first polynomial in the resulting basis is indeed free of
    t
    4
    . We ignore the quadratic factors, which are positive:
    In[]:=
    relate123=-Last[gbOmit4]
    Out[]=
    15+3
    2
    t
    1
    -16
    t
    1
    t
    2
    -
    2
    t
    2
    +3
    2
    t
    1
    2
    t
    2
    -8
    t
    1
    t
    3
    -16
    t
    2
    t
    3
    +8
    t
    1
    2
    t
    2
    t
    3
    +3
    2
    t
    3
    -
    2
    t
    1
    2
    t
    3
    +3
    2
    t
    2
    2
    t
    3
    -
    2
    t
    1
    2
    t
    2
    2
    t
    3
    Eliminating
    t
    2
    is similar and has a slightly simpler result:
    In[]:=
    gbOmit2=GroebnerBasis[{
    Z
    cos
    ,
    Z
    sin
    },{
    t
    2
    ,
    t
    1
    ,
    t
    3
    ,
    t
    4
    }][[1]]//Factor
    Out[]=
    1+
    2
    t
    1
    
    2
    1+
    2
    t
    3
    
    2
    1+
    2
    t
    4
    
    5+9
    2
    t
    1
    -3
    2
    t
    3
    +
    2
    t
    1
    2
    t
    3
    -8
    t
    3
    t
    4
    -8
    2
    t
    1
    t
    3
    t
    4
    -3
    2
    t
    4
    +
    2
    t
    1
    2
    t
    4
    -3
    2
    t
    3
    2
    t
    4
    +
    2
    t
    1
    2
    t
    3
    2
    t
    4
    
    In[]:=
    relate134=Last[gbOmit2]
    Out[]=
    5+9
    2
    t
    1
    -3
    2
    t
    3
    +
    2
    t
    1
    2
    t
    3
    -8
    t
    3
    t
    4
    -8
    2
    t
    1
    t
    3
    t
    4
    -3
    2
    t
    4
    +
    2
    t
    1
    2
    t
    4
    -3
    2
    t
    3
    2
    t
    4
    +
    2
    t
    1
    2
    t
    3
    2
    t
    4

    [A4]: Compatible values for
    t
    2
    and
    t
    3

    We introduce the new names
    x
    and
    y
    for
    t
    1
    and
    t
    4
    . We then compute the two possible values for
    t
    3
    in terms of
    x
    and
    y
    .
    In[]:=
    useXY={
    t
    1
    x,
    t
    4
    y}
    Out[]=
    {
    t
    1
    x,
    t
    4
    y}
    In[]:=
    t3vals=Solve[relate1340,
    t
    3
    ]/.useXY//Simplify
    Out[]=
    
    t
    3
    
    4(1+
    2
    x
    )y-
    15+22
    2
    y
    -9
    4
    y
    -
    4
    x
    2
    (-3+
    2
    y
    )
    +
    2
    x
    (22+60
    2
    y
    +6
    4
    y
    )
    (-3+
    2
    x
    )(1+
    2
    y
    )
    ,
    t
    3
    
    4(1+
    2
    x
    )y+
    15+22
    2
    y
    -9
    4
    y
    -
    4
    x
    2
    (-3+
    2
    y
    )
    +
    2
    x
    (22+60
    2
    y
    +6
    4
    y
    )
    (-3+
    2
    x
    )(1+
    2
    y
    )
    
    The radicand here is the discriminant of relate134 as a quadratic polynomial in
    t
    3
    , except that Simplify has pulled out a factor of
    Sqrt[4]2
    :
    In[]:=
    radQ=Discriminant[relate134,
    t
    3
    ]/4/.useXY
    Out[]=
    15+22
    2
    x
    -9
    4
    x
    +22
    2
    y
    +60
    2
    x
    2
    y
    +6
    4
    x
    2
    y
    -9
    4
    y
    +6
    2
    x
    4
    y
    -
    4
    x
    4
    y
    The possible values for
    t
    2
    are symmetric:
    In[]:=
    t2vals=t3vals/.{
    t
    3
    ->
    t
    2
    ,xy,yx}
    Out[]=
    
    t
    2
    
    4x(1+
    2
    y
    )-
    15+22
    2
    x
    -9
    4
    x
    +(22+60
    2
    x
    +6
    4
    x
    )
    2
    y
    -
    2
    (-3+
    2
    x
    )
    4
    y
    (1+
    2
    x
    )(-3+
    2
    y
    )
    ,
    t
    2
    
    4x(1+
    2
    y
    )+
    15+22
    2
    x
    -9
    4
    x
    +(22+60
    2
    x
    +6
    4
    x
    )
    2
    y
    -
    2
    (-3+
    2
    x
    )
    4
    y
    (1+
    2
    x
    )(-3+
    2
    y
    )
    
    To see which signs for the radicals in
    t
    2
    and
    t
    3
    are compatible, we substitute our two solutions for each into both
    Z
    cos
    and
    Z
    sin
    . It turns out that
    Z
    cos
    doesn’t care about the signs, but
    Z
    sin
    insists that the signs in the formulas for
    t
    2
    and
    t
    3
    be the same:
    In[]:=
    MatrixForm[Table[SameQ[0,Simplify[
    Z
    cos
    /.useXY/.t2vals[[i]]/.t3vals[[j]]]],{i,2},{j,2}]]
    Out[]//MatrixForm=
    
    True
    True
    True
    True
    
    In[]:=
    MatrixForm[Table[SameQ[0,Simplify[
    Z
    sin
    /.useXY/.t2vals[[i]]/.t3vals[[j]]]],{i,2},{j,2}]]
    Out[]//MatrixForm=
    True
    False
    False
    True

    [A5]: The first fundamental form of
    (
    E
    5
    ,
    e
    5
    )

    We denote the “bag-of-edges first fundamental form” as “bfff”.
    In[]:=
    halfAngleSubs={xTan[ξ/2],yTan[η/2]};
    In[]:=
    halfAngleBackSubs={Tan[ξ/2]x,Tan[η/2]->y,Cot[ξ/2]1/x,Cot[η/2]1/y};
    We compute “bfff” by plugging those phases into the definition of the bag-of-edges metric in (1). That computation itself is straightforward; but we don’t print the result, since it is lengthy and opaque.
    The following line then takes perhaps 15 seconds to do the crucial simplification:
    In addition, the formulas then lose some of their symmetry. For example, it is no longer obvious that swapping ξ with η swaps bE with bG — though that must still be the case:
    We compute the Gaussian curvature via the Brioschi formula (Alfred Gray [8], page 533):
    As in [A4], we replace the secants that differentiation has introduced with tangents and then FullSimplify without trigonometry.
    The curvature is intriguingly constant along the ditri-pentagon loops. We here consider three of them:
    Since FullSimplify didn’t manage to simplify that third one unaided, we give a hint:

    [A7]: The basics of solving the ODE in (9)

    We here compute the isothermal coordinates depicted in the bottom images of Figures 12 and 14. This involves solving the differential equation in (9) using machine precision — that is, IEEE double-precision floating point.
    Section 10.2 explains that, in some regions, we need a version of Sqrt in (9) that tracks the motion of the solution process along the two-sheeted Riemann surface of the square-root function, rather than just using the negative real axis as a fixed branch cut. To achieve this, the function memSqrt stores the last square root that it computed in the variable lastSqrt, which is dynamically bound by a Block in which the call to NDSolveValue occurs. For Figure 12, we have the following code:
    We could instead have computed bEr, bFr, and bGr as linear combinations of the unrotated bE, bF, and bG. Since the rotated variables μ and ν are related to the unrotated variables ξ and η by
    we have
    And indeed, corresponding entries of the matrices on the left and the right are equal:
    Given the rotated first fundamental form, the code for Figure 14 is no surprise:

    Solving the ODE in (9) accurately

    Along with IEEE doubles, Mathematica provides an arbitrary-precision floating point whose basic operations use significance arithmetic, a cheaper cousin of interval arithmetic, to bound the round-off errors that accumulate during easy stuff, such as arithmetic and the standard transcendental functions [18]. For sophisticated stuff, however, such as numerical integration or the numerical solution of differential equations, you specify a PrecisionGoal and an AccuracyGoal, goals for which you would like Mathematica to strive. You also specify the WorkingPrecision to be used during the calculation.
    Mathematica uses Accuracy to mean absolute accuracy, while Precision means relative accuracy. Given goals for both, Mathematica strives to achieve the easier of the two goals. We set our PrecisionGoal to Infinity, which tells Mathematica to focus on our AccuracyGoal. Mathematica seems able to solve Gauss’s ODE in (9) to whatever accuracy we like, but convincing it to do so is a bit complicated. (The following code strives for high accuracy without worrying about the computation time.)
    We don’t need the coordinates of the vertex v2, but here is a pair of valid formulas — one of many choices:

    Appealing to nonlinear least squares

    I evaluated polyApprox[deg] for each deg in the following sequence:
    Since the last few of those calls took many hours to compute, we here simply read in the results from the data file ApproximatingPolynomials.m, where I saved them. (We assume that file is available in the same directory from which this notebook file, “EquilateralPentagonMap.nb”, was read.) After reading that data file, the value poly[deg] will be the association that I computed via the call polyApprox[deg].
    The call polyApprox[250] took about 95 hours on my Windows PC, by the way. Computing the sample points and doing the nonlinear least squares took the first 30 of those, while computing the errors at the test points took the rest. The call polyApprox[225] took just 40 hours, and the times went down quickly from there.
    So the eventual convergence seems to be geometric. Fitting a line to the last five data points gives us:
    So a unit increase in the degree reduces the radial error by a factor of about 2.56. When the degree is small, a unit increase does a bit better than that; and there is a clear period-3 oscillation in the results:
    The radial-error functions pass through zero about deg times overall, when we combine the short and long legs; and they come reasonably close to equioscillating. In the following sample plots, the radial errors for the short leg appear above those for the long leg, each plotted from its left to its right:

    [A9]: The procedure PoincareCoordinates

    Other output regions

    The code

    We now give the code in the file “EquilateralPentagonMap.m”. It implements the routine PoincareCoordinates, distilling the above into a small package — well, small if we don’t charge for the data file ApproximatingPolynomials.m.
    We convert to rotated coordinates using the linear combinations discussed in [A7], expressing the results in terms of the rotated variables μ and ν:
    We next solve Gauss’s ordinary differential equation (9) in those rotated coordinates. We copy in the code for isoGoal and isoCare from [A8], since we are striving for high accuracy:
    The second argument to PoincareCoordinates is the desired accuracy. We add 3 to that, to allow for rounding errors in the subsequent arithmetic and just to be on the safe side; and we compare the result to the computed max errors of our precomputed approximating polynomials, choosing the first of those whose accuracy is adequate. If none of those has sufficient accuracy, we issue a warning message.
    The code for PoincareCoordinates in the non-error cases comes next. If the accuracy parameter is not specified, we default it to 17 decimal digits, which is roughly what you get in IEEE double precision.
    Note that telling isoCare to solve the ODE in (9) with accuracy acc is adequate, since cPoly divides by about 6, which improves the accuracy of its output.

    A few tests

    The π/2 vertex provides a more interesting test, since it is a long way from the horizontal axis during the solution of the ODE in (9). Computing the Poincare coordinates of the transformed π/2 vertex to 40 digits takes a couple of seconds, while computing it to 100 digits takes more than a minute:

    [A10]: Acknowledging my debt to Harry Calkins

    Back in 2021, when I first tried to solve the ODE in (9) accurately, I found that setting my AccuracyGoal higher than about 30 led to lots of warning messages from NDSolveValue — and associated inaccurate results. Fortunately, Cliff Hastings from Wolfram put me in touch with Harry Calkins. And Harry, bless him, managed to “touch base with the engineer who is at the top of the list when it comes to working with NDSolve”. That engineer advised me to try the method “StiffnessSwitching” — and hurrah! Things got vastly better! Even with that method, I still got sporadic warning messages, but only a few of them; and I had no difficulty programming around those, as described in [A8]. So, by exploiting Harry’s hint, I became able to solve the ODE in (9) with my accuracy limited only by my patience.
    Of course, the method “StiffnessSwitching” is described in the documentation for NDSolve; but that documentation is voluminous, and I had not studied it in sufficient depth to discover “StiffnessSwitching” on my own. Indeed, I am still fuzzy about exactly what makes an ODE “stiff” and why the ODE in (9) has “small intervals of stiff behavior”.
    By the way, the sporadic warning messages that remain seem to stem from infelicities with the GNU Multiple Precision Arithmetic Library, which Mathematica uses for its high-precision calculations. (The problem might be in that library itself, or in how Mathematica invokes that library, or in how the routines of that library get compiled for my Intel processor, or whatever. I am no expert, and I haven’t investigated further.) When representing a high-precision floating-point number, that library stores 64 bits of the mantissa in each computer word. And my remaining warnings arise for values of WorkingPrecision that differ by integer multiples of 19.3, where 19.3 decimal digits are equivalent to 64 binary bits. Thus, the warnings probably arise at the boundaries where one additional word becomes necessary.