remarkRemark hypothesisHypothesis claimClaim Mathieu Functions: Historical
PerspectiveC. Brimacombe, R. M. Corless, and M. Zamir

# COMPUTATION AND APPLICATIONS OF MATHIEU FUNCTIONS: A HISTORICAL
PERSPECTIVE††thanks: Submitted to the editors DATE. This work was supported by
NSERC under grants numbered RGPIN-2020-06438 (RMC) and RGPIN-2019-04749 (MZ);
and, while RMC was visiting the Isaac Newton Institute during the programme
Complex Analysis: Tools, techniques, and applications, by EPSRC Grant #
EP/R014604/1.

Chris Brimacombe University of Toronto, Toronto, Canada ().
<EMAIL_ADDRESS>Robert M. Corless School of Mathematical
and Statistical Sciences, Western University, London, Canada (, ).
<EMAIL_ADDRESS><EMAIL_ADDRESS>Mair Zamir33footnotemark: 3

###### Abstract

Mathieu functions of period $\pi$ or $2\pi$, also called elliptic cylinder
functions, were introduced in 1868 by Émile Mathieu together with so-called
modified Mathieu functions, in order to help understand the vibrations of an
elastic membrane set in a fixed elliptical hoop. These functions still occur
frequently in applications today: our interest, for instance, was stimulated
by a problem of pulsatile blood flow in a blood vessel compressed into an
elliptical cross-section. This paper surveys and recapitulates the historical
development of the theory and methods of computation for Mathieu functions and
modified Mathieu functions and identifies some gaps in current software
capability, particularly to do with double eigenvalues of the Mathieu
equation. We demonstrate how to compute Puiseux expansions of the Mathieu
eigenvalues about such double eigenvalues, and give methods to compute the
generalized eigenfunctions that arise there. In examining Mathieu’s original
contribution, we bring out that his use of anti-secularity predates that of
Lindstedt. For interest, we also provide short biographies of some of the
major mathematical researchers involved in the history of the Mathieu
functions: Émile Mathieu, Sir Edmund Whittaker, Edward Ince, and Gertrude
Blanch.

###### keywords:

Mathieu functions; modified Mathieu functions; historical survey; computation
of Mathieu functions; double eigenvalues; Puiseux series

01-02, 33-02, 33F05

##  1 Introduction

What is the sound of an elliptic drum? In a memoir presented at the Sorbonne
in 1868, Émile Mathieu showed the way to find the answer, when he described
the solution of a mechanical vibration problem characterized by an elliptic
boundary. The memoir was groundbreaking in that it introduced a new
differential equation whose eigenvalues and corresponding periodic solutions
led to the definition of a new class of functions. In 1912 Whittaker named
these new functions in honour of their discoverer: the differential equation
is now known as the Mathieu equation and the $\pi$ and $2\pi$ periodic
solutions (and only these periodic solutions) are known as the Mathieu
functions. While work on theoretical and analytical aspects of Mathieu
functions has continued since the introduction of these functions in 1868, the
focus has shifted in recent years to work on numerical and computational
aspects. This development, driven by steep advances in digital technology, has
given rise to heavy reliance today on “packaged software” for the evaluation
of Mathieu functions. This practice comes with the risk of concealing as yet
unresolved (or at least not completely resolved) analytical and computational
issues involved in the use of Mathieu functions.

Modern approaches to solving such vibrational problems are more likely to be
direct. See [36] for an exemplar of this approach. Spectral expansions using
special functions do remain of interest for approximation, however; see [82]
for an exemplar of this approach. But regardless of application, the Mathieu
functions themselves remain of interest because the equation is simple enough
to occur in a wide variety of contexts, and there is a need to compute them
directly.

We review methods of computation of periodic Mathieu functions and note that
all the methods we review fall short at some “exceptional” or “double” points,
as noted both by [14] and by [45]. Although those objections were made more
than fifty and forty years ago respectively, and theoretical work on these was
completed in [60], we believe that there is still no fully satisfactory code
available, as we will detail. We will discuss what to do in practice at these
double points, where the eigenvalues merge and we lose one independent
eigenfunction and thus lose completeness of our set of eigenfunctions.

While we survey many results here, we cannot cover everything. Instead we list
several references that contain many details (and indeed many important
results) that we omit. The first and most important reference is to Chapter 20
of the classic [1], written by Gertrude Blanch, and its successor, Chapter 28
of the Digital Library of Mathematical Functions (henceforth DLMF,
https://dlmf.nist.gov/28.1) written by G. Wolf. The DLMF is intended to be an
updated replacement for [1] in this computer age. Both are (perfectly legally)
free online. The DLMF gives pointers to many of the alternative notations used
elsewhere in the rather substantial literature. The original Chapter 20 of [1]
has yet more references. We will mention others as we go.

One of our favourite reference works for methods of computation of special
functions, namely the very thorough and carefully laid out book [38],
unfortunately does not include Mathieu functions as examples (indeed, there
may be a reason for that, which we will discuss). However, many of the
techniques that are described in that book are potentially applicable for
Mathieu functions, and we will see variants of some of them used in this
present paper.

###  1.1 Organization of the paper

In section 2 we briefly outline three applications of Mathieu functions,
including the one that motivated us to perform this study. We then
recapitulate in section 3 some of the historical development of Mathieu
functions. We are not historians of mathematics or science, but we have done
our best. In section 4 we look at a method to compute double points and the
corresponding eigenvalues, and Puiseux series, about those points, for the
eigenvalues. In section 5 we look at algorithms for computing solutions of the
Mathieu equations, including Mathieu functions (eigenfunctions). We finish
that section with a discussion of how to compute generalized eigenfunctions
for double eigenvalues, which are needed for completeness. We provide
concluding remarks in section 6. In appendix A we provide more details of the
application that motivated us to undertake this work. In appendix C we discuss
confocal ellipses, and give a singular perturbation argument relating Mathieu
functions to Bessel functions in the limit as ellipses become circles. In
appendix D we compare Mathieu’s perturbative solution to the Mathieu equation
with that produced by a computer algebra system, and apart from minor errors
and typos in his paper confirm his results. Finally, in appendix E we provide,
as an homage to the great scientific table-makers, a table of Puiseux
expansions about double points. Of course, the computer-readable version at
https://github.com/rcorless/Puiseux-series-Mathieu-double-points (together
with the Maple code used to generate it) is much more useful nowadays; still,
there’s something to be said for looking at a collection of related numbers.
We remark that printing that table nicely with LaTeX—at least, nicely enough
compared to old-style high-quality production—made for an interesting
challenge, and some quirks remain in the table as “Easter Eggs” for any that
are willing to hunt for them.

##  2 Applications

###  2.1 Columns and Strings Under Periodically Varying Forces

One of the more surprising physics experiments in the undergraduate curriculum
is the stabilization of a vertical hinged column by up-and-down vibration at
the right frequency. There are many YouTube videos demonstrating this, and,
after reading this paper, the reader may choose to search some of those videos
out: a good set of keywords to search with is “stability of the inverted
pendulum”. This phenomenon is explained using topological terms in [55]. The
underlying mathematics has been known at least since 1928 [87], and is
concerned directly with the Mathieu equation, if the periodic force is a
simple sine or cosine.

There are similar, though more complicated, problems which still need the
Mathieu equation. We will briefly describe the model studied in [57], namely
the stability of a column under periodically-varying compression, or of a
string under tension with a periodically-varying tension force. These problems
have infinitely many degrees of freedom, as opposed to the simple inverted
pendulum above. Consider for instance the case of a column under compression.
If the force $F$ is greater than the Euler load, then it is well known that
the column can buckle; the question at issue here is if $F=P+H\cos\omega t$,
which consists of a steady part plus a periodically-varying part, can one
choose $H$ and $\omega$ so that, even if $P$ is larger than the Euler load,
the column remains stable?

The answer is a qualified yes; one qualification is that, at least part of the
time, the force must be less than the Euler load, which makes sense. More, the
qualifications require the study of the stability of the solutions to the
Mathieu equation, and therefore require a good knowledge of the periodic
solutions of the Mathieu equation. The paper [57] also studies the motion of a
string with a periodically-forced tension with a similar model; see figure 1.
The equations of motion they derive are

$EI\frac{\partial^{4}w}{\partial x^{4}}-F(t)\frac{\partial^{2}w}{\partial
x^{2}}+m\frac{\partial^{2}w}{\partial t^{2}}=0\>,$

where they take $F(t)=P+H\cos\omega t$ and the familiar Young’s modulus $E$
characterizes the relationship between stress and strain, and $I$ is the
moment of inertia of the body’s cross-section.

Figure 1: A flexible string with a time-varying tension $F(t)$. The figure
indicates that the deformation $w(t)$ is confined to the plane. We model this
figure on Figure 1 of [57].

“In all of these problems the Mathieu equation (more properly, a sequence of
Mathieu equations in the continuous systems) plays a central rôle, since the
decision as to stability depends on the character of the solutions to such
equations.” [From the introduction in [57].]

In fact the study of the stability of these systems requires more than we are
going to cover in this paper: it needs the Floquet theory, which we only
lightly touch on.

###  2.2 Pulsatile Flow in Tubes of Elliptic Cross Sections

Our motivation for this present paper originated from a problem in pulsatile
blood flow. Under normal circumstances, blood flow occurs in vessels of
circular cross sections, but under a number of important pathological
conditions the vessels are deformed by external forces to the effect that
their cross sections are no longer circular. In a separate study (currently in
progress) we are examining this phenomenon using an elliptic cross section as
a simple, mathematically tractable, departure from a circular cross section,
with the main focus being on the hemodynamic consequences of that departure.
In the present paper our focus is on the mathematical and computational
consequences.

Fluid flow in a tube (see figure 2) is in general governed by a simplified
form of the Navier-Stokes equations [92]. If the tube is straight and of
uniform cross section, the flow can be described in terms of a single velocity
component $U$ along the axis of the tube. If the flow is pulsatile, as in the
cardiovascular system, the velocity $U$ consists of a steady part $u_{0}$ plus
an oscillatory part $u(t)$ such that

(1) $U=u_{0}+u(t)\>,$

where $t$ is time. Only the oscillatory component of velocity, $u(t)$, is
relevant to the present discussion since the transition from Bessel to Mathieu
functions occurs in the governing equations of this component of the flow as
the cross section of the tube changes from circular to elliptic. We give a
brief sketch here, with more details of our motivating application in appendix
A. In that appendix the focus is on the transition from Bessel to Mathieu
equations in the case of pulsatile flow in a tube as the cross section of the
tube transitions from circular to elliptic geometry.

While pulsatile flow in tubes of circular cross sections is governed by Bessel
equations, with solutions in terms of Bessel functions, the corresponding
situation in tubes of elliptic cross sections involves the Mathieu equation
and Mathieu functions. From a mathematical perspective, this transition is not
only a matter of curiosity but also a matter of practical importance because
of the comparative difficulties involved in the numerical computation of these
two kinds of functions. The difficulties and pitfalls involved in the
computation of Mathieu functions are the focus of much of the present paper.

Figure 2: Flow in a tube of elliptic cross-section. The elliptic cross-section
is meant to model some pathologies where the blood vessels are deformed so
that they are no longer circular. This figure is modelled after Figure 3.1,
page 45, of [92].

###  2.3 Vibrating membrane bounded by an ellipse

Figure 3: Contours and nodal lines of a possible pure vibration mode of an
elliptic drum. The aspect ratio of the ellipse is $3:5$. At the time of this
snapshot, the rightmost loops in this mode have positive values for the
contours, and all adjacent cells have opposite signs. Notice the hyperbolic
nodal lines where the membrane does not move, separating the cells. We have
suppressed here the details of which Mathieu functions were used to produce
this figure (we give them in section 3.1.5), but we followed the method
outlined by Mathieu in his 1868 paper [58]. See [22] for complete details.

One of the simplest physical problems whose solutions involve the Mathieu
equation and the Mathieu functions is the sound made by an elliptical drum.
This in fact was the problem Mathieu himself studied, and in section 3.1 we
will look in detail at how he solved it. The physical problem being modelled
is, in fact, remarkably easy to define: imagine an elliptical hoop, fixed
immovably, and a thin homogeneous membrane stretched tight across the hoop,
making a drum. What are the natural modes of vibration of this drum, and how
could we describe mathematically its motion once struck?

One possible natural mode for one particular drum is pictured in figure 3.
There we see contours of vibration, including contours where there is no
motion (the so-called “nodal lines”). This figure was drawn using our own
software to compute the relevant Mathieu functions, but many software packages
exist which could do this.

##  3 Historical overview, introducing notions and notation

We introduce the Mathieu equation and the Mathieu functions in historical
order, by discussing the contributions of several of the main researchers
involved. The result is a tour of several aspects of late nineteenth-century
and early twentieth-century mathematics. We also give some biographical
details of these main figures. For concreteness in the discussion to follow,
here is the Mathieu equation in one common modern notation:

(2) $\frac{d^{2}y}{dx^{2}}+\left[a-2q\cos(2x)\right]y=0\>.$

The parameter $q$ is given by the physics or the geometry of the specific
problem at hand; the eigenvalue $a$ must be calculated in order to ensure
periodicity of $y$, given $q$ and a desired order. The so-called modified
Mathieu equation is related to equation (2) by the transformation $z=\pm ix$
(the sign makes no difference):

(3) $\frac{d^{2}y}{dz^{2}}-\left[a-2q\cosh(2z)\right]y=0\>.$

The even solutions are conventionally written as
$\mathrm{Ce}_{g}(z,q)=\mathrm{ce}_{g}(\pm iz,q)$. The eigenvalues for even
solutions are conventionally written $a_{g}(q)$. The odd solutions are written
similarly, as $\mathrm{se}_{g}(z,q)$ and $\mathrm{Se}_{g}(z,q)$, and the odd
eigenvalues are conventionally written $b_{g}(q)$. Here $g$ is a nonnegative
integer, and the solutions split into further classes if $g$ is itself even or
odd, as we will see. The use of the letter $g$ for an integer contradicts the
usual $I-N$ convention in use nowadays, from Fortran; but Mathieu used the
letter $g$ in this way and we find it convenient when the parity of $g$ is
unspecified.

If we write the general solution of the Mathieu equation (with no initial or
boundary conditions applied) with arbitrary constants $\alpha$ and $\beta$ as

(4) $\displaystyle y(x)=\alpha w_{I}(x;a,q)+\beta w_{II}(x;a,q)$

using the fundamental pair of solutions satisfying $w_{I}(0;a,q)=1$ with
$w_{I}^{\prime}(0;a,q)=0$ and $w_{II}(0;a,q)=0$ with
$w_{II}^{\prime}(0;a,q)=1$, (using the notation of the DLMF and where ′
denotes $d/dx$) then the general solution of the Modified Mathieu equation can
be written

(5) $\displaystyle y(z)=\alpha w_{I}(\pm iz;a,q)+\beta w_{II}(\pm iz;a,q)\>.$

Some software packages denote these functions $C$ and $S$ respectively. When
$a(q)$ or $b(q)$ is an eigenvalue, the periodic Mathieu functions must satisfy
(now allowing $x$ to be complex and renaming it $z$)

(6) $\displaystyle\mathrm{ce}_{m}(z,q)$
$\displaystyle=c_{m}w_{I}(z;a_{m}(q),q)$ (7)
$\displaystyle\mathrm{se}_{m}(z,q)$
$\displaystyle=s_{m}w_{II}(z;b_{m}(q),q)\>,$

for some normalization constants $c_{m}$ and $s_{m}$. In this paper we take
those normalization constants to be $1$.

In theory, the use of Mathieu functions and modified Mathieu functions in the
solution of the aforementioned physical problems is attractive because the
functions are analogous to harmonic functions, and expansions in terms of them
can be efficient in comparison with direct numerical solution of the PDE
model. In practice, there are annoying difficulties: the available software
might be restricted to real arguments, or use a different normalization than
the one desired (the authors of [32] state that there are at least three
normalizations in common use; they themselves use the same normalization that
we do here), or the software may fail to be accurate for “difficult” values of
the problem parameters, say for large values. A more serious problem is the
approximation properties of the expansion itself, followed by the numerical
stability of the expansion. Approximation properties are explored in [77] for
real $q$, with all the power of the Sturm–Liouville theory (which indeed
Mathieu himself used in his 1868 paper). Numerical stability, on the other
hand, has received less attention.

Now that we have sketched where we want to go with Mathieu functions, let us
recapitulate their development.

###  3.1 Émile Léonard Mathieu (1835–1890)

Mathieu began his 1868 discussion [58] of the vibrations of an elastic
membrane held fixed by a hoop in the shape of an ellipse by first considering
the simpler problem when the hoop is, in fact, circular111This memoir was
translated from its nineteenth century French for us by Dr. Robert H. C. Moir,
and the translation—which we believe may be of interest on its own—has been
made available [59].. Mathieu’s discussion of the circular case starts with
the PDE

$\frac{\partial^{2}w}{\partial t^{2}}=m^{2}\left(\frac{\partial^{2}w}{\partial
x^{2}}+\frac{\partial^{2}w}{\partial y^{2}}\right)$

(Mathieu used $d$ and not our modern Russian $\partial$ for partial
derivatives) and then transformed $x=r\cos\alpha$, $y=r\sin\alpha$ to polar
coordinates to get

$\frac{\partial^{2}w}{\partial t^{2}}=m^{2}\left(\frac{\partial^{2}w}{\partial
r^{2}}+\frac{1}{r}\frac{\partial w}{\partial
r}+\frac{1}{r^{2}}\frac{\partial^{2}w}{\partial\alpha^{2}}\right)\>.$

Thereafter Mathieu used what is now the standard method of separation of
variables for a pure oscillation $w=\sin(2\lambda mt)u(r,\alpha)$ and
$u(r,\alpha)=P(\alpha)Q(r)$ to give a harmonic equation for $P(\alpha)$ and is
equivalent to what we now call Bessel’s equation for $Q(r)$:

${r}^{2}{\frac{d^{2}}{d{r}^{2}}}Q\left(r\right)+r{\frac{d}{dr}}Q\left(r\right)-\left({n}^{2}-4\,{\lambda}^{2}{r}^{2}\right)Q\left(r\right)\>.$

[Bessel’s equation in standard form is
$z^{2}d^{2}y/dz^{2}+zdy/dz+(z^{2}-\nu^{2})y=0$; see DLMF 10.2.1.] One can then
write the solution to the original vibration problem as a linear combination
of products of these eigenfunctions, and determine the unknown coefficients by
matching to the boundary conditions using orthogonality.

Although Bessel (1784–1846) did his work first (and although it was actually
Daniel Bernoulli who first identified this equation, even earlier) Mathieu did
not call this “Bessel’s equation”, or give it a name at all, but merely solved
it in series in what must have been common practice at the time. He then
demonstrated that one must find the zeros of the various series for the Bessel
functions, and cited Bourget’s Memoirs [17] for a method for doing so, in
order to identify the eigenvalues and eigenmodes of the vibration of the
membrane in the circular case.

Following the same strategy in the elliptic case, but using a confocal
elliptic transformation222Mathieu initially spelled the hyperbolic functions
out explicitly, and later used “Pour simplifier” the notation
$E(\beta)=(e^{\beta}+e^{-\beta})/2$ for $\cosh(\beta)$ and
$\mathcal{E}(\beta)=(e^{\beta}-e^{-\beta})/2$ for $\sinh(\beta)$. Johann
Heinrich Lambert had already in the 18th century introduced the notation which
we use today; it is interesting that it was not universally used in the 19th.
$x=c\cosh(\beta)\cos(\alpha)$ and $y=c\sinh(\beta)\sin(\alpha)$ where $2c$ is
the distance between the foci of the ellipse333Notice that this coordinate
transformation is singular if $c=0$. Therefore the connection of Mathieu
functions to Bessel functions, while present, requires art to tease out. See
Appendix C., Mathieu arrived at an equation that on separation gives two
equations equivalent to those now known as the Mathieu equation and the
modified Mathieu equation:

$-\frac{1}{P}\frac{d^{2}P}{d\alpha^{2}}+4\lambda^{2}c^{2}\cos^{2}(\alpha)=+\frac{1}{Q}\frac{d^{2}Q}{d\beta^{2}}+4\lambda^{2}c^{2}\cosh^{2}(\beta)\>.$

“Comme le premier membre ne peut renfermer que $\alpha$, et le second que
$\beta$, ils sont égaux a une même constante $N$.” Readers may be pleased to
learn that separation of variables reads the same in the 21st century as it
did in the 19th, mutatis mutandis.

There is one remaining difference to the modern notation, and that is the use
of $\cos^{2}\alpha$ and $\cosh^{2}\beta$. Using double-angle identities,
Mathieu later in this same paper transformed these to something that we write
in modern notation as

(8) $\displaystyle\frac{d^{2}P}{d\alpha^{2}}+(a-2q\cos 2\alpha)P$
$\displaystyle=0$ (9) $\displaystyle\frac{d^{2}Q}{d\beta^{2}}-(a+2q\cosh
2\beta)Q$ $\displaystyle=0\>.$

Notice that these two equations can be transformed into each other by the
change of variable $\beta=i\alpha$. Here $q=\lambda^{2}c^{2}$ is a relabeling
of the parameter that contained the physics, in Mathieu’s case the elastic
constant, as well as the focal distance $c$, and $a$ (which Mathieu called
$R$) is the separation constant, adjusted from Mathieu’s earlier variable $N$
by changing from $\cos^{2}\alpha=(1+\cos 2\alpha)/2$ and
$\cosh^{2}\beta=(1+\cosh 2\beta)/2$ to the double-angle forms, for a reason
that will become apparent. Mathieu used $h^{2}$ where we have $q$ here, and
that notation is still occasionally used.

The boundary conditions of the original problem reflect the elliptic geometry.
The angular coordinate $\alpha$ runs from $0$ to $2\pi$ (or, of course, from
$-\pi$ to $\pi$), requiring periodicity. Often the problem is taken to have
$\pi$-symmetry, to reflect the symmetry between the top and the bottom of the
ellipse.

####  3.1.1 Basic properties: orthogonality, eigenvalues

Mathieu noted that the theory of Jacques Charles François Sturm444Sturm was
the successor of Poisson in the chair of mechanics in the Faculté de Sciences,
Paris. This fact has a certain poignancy when read together with the obituary
of Mathieu [28], where we learn that Mathieu had desired that chair.
(1803–1855) implied that, for real $q$, what we now call the Sturm–Liouville
form of the Mathieu equation could be written as

(10)
$\displaystyle\mathcal{L}(y):=\frac{d}{dx}\left(L(x)\frac{dy}{dx}\right)-G(q;x)y=Nr(x)y(x)\>,$

(here, trivially, $L(x)=1$, $r(x)=1$, and $G(q;x)=2q\cos 2x$, while the
eigenvalue $N$ has been moved to the right hand side) and that several useful
properties naturally followed. See Chapter 39 of [54], and in particular Lemma
39.4 there; note that a necessary element of the proof of the properties
described in that lemma is that $G(q;x)$ be real. See also [71] for numerical
treatment of Sturm–Liouville problems.

Strictly speaking, Mathieu had to extend the now-classical theory to the case
of periodic boundary conditions in order to establish what he needed. Indeed,
he also showed how to compute the eigenvalues by an interesting perturbative
argument which we will take up in section 3.1.3. Before that, he established
an inequality around each neighbourhood of $g^{2}$, the square of an integer;
to us, this seems a convincing argument for the existence of each eigenvalue
(and thus of the infinite collection) for small real $q$.

The Sturm–Liouville theory shows that given a real value of $q$ there are a
countable number of eigenvalues $N$ of the Mathieu equation, conventionally
written either $N=-a_{k}(q)$ or $N=-b_{k}(q)$ depending on the type, which if
$q>0$ can be arranged in the sequence
$a_{0}(q)<b_{1}(q)<a_{1}(q)<b_{2}(q)<a_{2}(q)<\cdots$, with the eigenvalues
tending to infinity. Next, to each eigenvalue there corresponds an
eigenfunction, unique up to normalization. These eigenfunctions are now called
the Mathieu functions, and they come in four classes. If they are even and of
period $\pi$ they are denoted by $\mathrm{ce}_{2k}(x)$, for $k\geq 0$. If they
are even and of period $2\pi$ they are denoted by $\mathrm{ce}_{2k+1}(x)$, for
$k\geq 0$. If they are odd and of period $\pi$ they are denoted by
$\mathrm{se}_{2k}(x)$ for $k\geq 1$. If they are odd and of period $2\pi$ they
are denoted by $\mathrm{se}_{2k+1}(x)$ for $k\geq 1$. See figure 4 for a
representative graph of a few low-frequency Mathieu functions, with $q=1.5$.
Mathieu established that these eigenfunctions are orthogonal with respect to
the bilinear form defined by

(11) $\displaystyle\left\langle
y_{k},y_{\ell}\right\rangle:=\int_{0}^{2\pi}y_{k}(x)y_{\ell}(x)\,dx=\mathrm{const}\cdot[\mathcal{C}\
{\bf\mathrm{and}}\ k=\ell]\>.$

We prove this in detail in appendix B. By the proposition $\mathcal{C}$ we
mean here a proposition that is true if $y_{k}$ and $y_{\ell}$ are in the same
class; that is, both are even and of period $\pi$, both are odd and of period
$\pi$, both are even and of period $2\pi$, or both are odd and of period
$2\pi$. We have used here the Iverson convention $[\cdot]$ to mean $1$ if the
proposition inside the brackets is true and $0$ otherwise, because we prefer
its generality over the more restrictive Kroenecker delta symbol, and do not
think it will be confused with citations [53]. In other words, eigenfunctions
of one class are orthogonal to eigenfunctions of any other class using this
bilinear form, which is an inner product if $q$ and $a$ are real. If
appropriate, one could integrate to $\pi$ instead of $2\pi$. The exact method
used for normalization is a matter of convention. Finally, these classes of
orthogonal eigenfunctions are complete: every reasonable555As with Fourier
series, one can prove convergence for a quite wide class of functions;
however, for utility and rapid convergence, one needs a large number of
continuous derivatives, and the more smoothness the better. even/odd
period-$p$ function $f(x)$ can be expanded in a convergent series

(12) $\displaystyle f(x)=\sum_{k\geq 0}\alpha_{k}y_{k}(x)\>.$

Periodic functions which are neither even nor odd can of course be written as
a sum of an even function and an odd function

$f(x)=\frac{f(x)+f(-x)}{2}+\frac{f(x)-f(-x)}{2}$

and thus will use both classes of the given period in its expansion. Since the
coefficients $\alpha_{k}$ can be determined by orthogonality, this series is
expected to be practical and to enable us to find computationally useful
solutions to the original PDE by matching the boundary condition at the edge
of the ellipse.

(a) Some $\mathrm{ce}_{2k}(\alpha)$ for $q=1.5$

(b) Some $\mathrm{se}_{2k}(\alpha)$ for $q=1.5$

(c) Some $\mathrm{ce}_{2k+1}(\alpha)$ for $q=1.5$

(d) Some $\mathrm{se}_{2k+1}(\alpha)$ for $q=1.5$

Figure 4: The first few Mathieu functions with $q=1.5$. The normalization
shown here has $y(0)=1$ in the case of the even functions and
$y^{\prime}(0)=1$ in the case of the odd functions. Period $\pi$ functions are
shown in the top row, period $2\pi$ functions in the bottom row.

We emphasize that Mathieu only established this for real $q$. The case of
complex $q$ is a different matter, as noted by [60] and by [14]. The Sturmian
theory fails there (see appendix B) because there may be (and in fact are)
double eigenvalues, and in that case Blanch points out that for the
eigenfunction associated with the double eigenvalue
$\left<y_{k},y_{k}\right>=0$ and thus normalization by making the bilinear
form equal to a nonzero constant is not possible. We will return to this
later.

####  3.1.2 Normalizations

Mathieu normalized the solutions of his equation in a way that might seem
curious to modern eyes. First, he noted that it is easy to see that the
solution $y(x)$ of any linear second order ordinary differential equation
(ODE) may be written as $y=P_{1}+P_{2}$, where $P_{1}$ is either a maximum or
a minimum at $x=a$ while $P_{2}$ is zero at $x=a$ ($a$ is arbitrary). A
moment’s reflection shows that Mathieu was correct, and that this is true for
any $a$ in the domain of definition of $y(x)$. What is curious is that this
separates a second order equation into two parts, each of which imposes only
one condition: $P_{1}^{\prime}(a)=0$, or $P_{2}(a)=0$. Clearly any multiple of
$P_{1}$ or of $P_{2}$ will satisfy the same condition, but only one
combination of such functions will equal $y(x)$. This left Mathieu free to
normalize his functions in any way convenient to him, and he took great
advantage of it, in particular in his perturbative solutions. To normalize his
functions, he chose to make the coefficient of $\cos g\alpha$ in its series
expansion to be unity (and similarly the coefficient of $\sin g\alpha$ in the
case of odd eigenfunctions); see appendix D for a comparison with one modern
normalization. As Ince pointed out later in [49], this is not always possible
because for some values of $q$, even for some real values of $q$, the
coefficient of $\cos g\alpha$ is actually zero; but it is at least almost
always possible in the modern sense; that is, except on a set of measure zero
in parameter space.

Other normalizations are in use today: we use the universally-possible
normalization discussed at the end of section 20.5 in [1], namely,
$y(x)=a\mathrm{ce}_{g}(x)+b\mathrm{se}_{g}(x)$ and we specify that
$\mathrm{ce}_{g}(x)$ satisfy not only $\mathrm{ce}_{g}^{\prime}(0)=0$ but also
$\mathrm{ce}_{g}(0)=1$, and similarly $\mathrm{se}_{g}(0)=0$ and
$\mathrm{se}_{g}^{\prime}(0)=1$. As Blanch states in the aforementioned
reference, conversion between normalizations is “rather easy,” but we wanted
one that would always work.

However, in many published papers and codes the “norm” using the bilinear form
(11) is nominally enforced to be $\pi$ (or $2\pi$), which can (again) almost
always be done, but not universally: at double eigenvalues, the “norm” must
vanish. This is a little better than Mathieu’s normalization in that all such
exceptional values of $q$ must be complex: the norm will always be nonzero for
real $q$. But for some complex $q$, which we look at carefully in this present
paper, the “norm” does vanish. This means that in modern terms the norm is an
“indefinite” norm [4], and requires some care in handling. To emphasize, in
this paper we do not normalize by the bilinear form, but instead choose to
enforce initial conditions as above; this is also done elsewhere in the
literature, but not commonly.

####  3.1.3 Perturbative solution: anticipating Lindstedt

In the 1868 paper under discussion, Mathieu developed series solutions for the
first few eigenvalues, $a_{k}(q)$ and $b_{k}(q)$ in modern notation; in some
cases to sixth order in $q$ (twelfth order in $h$). It is interesting to note
that to do so he essentially used what we now know as anti-secularity: he
chose series coefficients in the eigenvalue expansion in order to eliminate
secular terms in the expansion for the eigenfunction and thereby enforced
periodicity of the solution. This notion is typically introduced nowadays as
the Lindstedt–Poincaré method. There are alternatives, now, too: one can
instead use the method of multiple scales, or in an even more modern way, use
renormalization. See for instance [23] and the references therein for more
details of those methods.

When Mathieu published his memoir in 1868, Anders Lindstedt was in his early
teens and his work on perturbation [56] was fourteen years in the future.
Mathieu might have good grounds for a claim to priority, even though (perhaps)
Lindstedt’s work was somewhat more general666Lindstedt’s method applies to
weakly nonlinear equations, which are linear if the small parameter is set to
zero. Lindstedt suggested simultaneous expansion of the eigenvalue. This
generates a sequence of linear equations to solve for subsequent terms, and
the overall process is not much different from what Mathieu did.. Mathieu’s
use of anti-secularity is clear, however, once one tries to retrace his steps;
it seems very natural, although Mathieu does not comment on it explicitly.
Indeed, his section 11 which details the perturbation solution reads more like
an informal summary of notes of how to proceed, with many details left out.
Nonetheless, using anti-secularity to enforce periodicity is exactly what he
did. He also made several elegant uses of his freedom to normalize in the
problem in order to reduce the labour involved. We have implemented his
solution in a computer algebra system, to retrace his steps and fill in the
details; some of our computations are compared to his in Appendix D.

Mathieu’s first computation was to find the even period $2\pi$ solution of
equation (8) when $q=h^{2}$ was small and the eigenvalue $a$ approached
$g^{2}$, the square of an unspecified integer. The solution in his notation
and with his normalization777Ince pointed out in [50] that this fails for some
$q$; we discuss this later. and to fewer terms than he calculated to is:

$\displaystyle\mathrm{ce}_{g}(\alpha)=$ $\displaystyle\cos
g\alpha+\left({\frac{\cos\left(g-2\right)\alpha}{4(g-1)}}-{\frac{\cos\left(g+2\right)\alpha}{4(g+1)}}\right){h}^{2}$
$\displaystyle+\left({\frac{\cos\left(g-4\right)\alpha}{32(g-2)\left(g-1\right)}}+{\frac{\cos\left(g+4\right)\alpha}{32(g+2)\left(g+1\right)}}\right){h}^{4}$
$\displaystyle+\left({\frac{\cos\left(g-6\right)\alpha}{384(g-4)\left(g-2\right)\left(g-1\right)}}+{\frac{\left({g}^{2}-4\,g+7\right)\cos\left(g-2\right)\alpha}{128(g-2)\left(g+1\right)\left(g-1\right)^{3}}}\right.$
(13)
$\displaystyle\left.-{\frac{\left({g}^{2}+4\,g+7\right)\cos\left(g+2\right)\alpha}{128(g+2)\left(g+1\right)^{3}\left(g-1\right)}}-{\frac{\cos\left(g+6\right)\alpha}{384(g+2)\left(g+3\right)\left(g+1\right)}}\right){h}^{6}+O(h^{8})$

As Mathieu noted, this series is valid only for large enough integers $g$. He
also correctly computed the corresponding eigenvalue (he called it $R$ in this
part of his paper) as

(14)
$a={g}^{2}+{\frac{{h}^{4}}{2(g-1)\left(g+1\right)}}+{\frac{\left(5\,{g}^{2}+7\right){h}^{8}}{32\,\left(g-2\right)\left(g+2\right)\left(g-1\right)^{3}\left(g+1\right)^{3}}}+\cdots\>.$

Mathieu then goes on to show how to compute perturbation solutions for
specific, smaller, frequencies $g$. See appendix D for details.

The idea of a series expression for the Mathieu functions was, of course,
natural for the time. Whether the idea of enforcing periodicity by expanding
the eigenvalue in series was original to Mathieu, we do not know; but its
presence in his paper certainly predates Lindstedt’s work.

For Mathieu, $q=h^{2}$ was real, and small (if the interfocal distance $2c$
was small). In many modern applications, $q$ might be complex, or large, or
both. It took many years of further research by others to go beyond these
series.

####  3.1.4 D-finite, or ‘holonomic’, formulation

After finding these perturbation expansions, Mathieu took a different
approach: he changed variables, first with $\nu=\cos\alpha$, whereupon the
Mathieu equation becomes (equation 28.2.3 in the DLMF)

(15)
$(1-\nu^{2})\frac{d^{2}P}{d\nu^{2}}-\nu\frac{dP}{d\nu}+(a+2q(1-2\nu^{2}))P=0$

and alternatively by $\mu=\sin\alpha$, (Mathieu used $\nu^{\prime}$ rather
than $\mu$, which we originally kept; but a referee pointed out that this was
confusing, and on second thought we agreed) whereupon the Mathieu equation
becomes

(16)
$(1-\mu^{2})\frac{d^{2}P}{d\mu^{2}}-\mu\frac{dP}{d\mu}+(a-2q(1-2{\mu}^{2}))P=0\>.$

The DLMF gives yet another algebraic form, using the change of variables
$\zeta=\sin^{2}\alpha$. These are interesting for several reasons, and we will
mention in section 3.2 some of the further properties that can be deduced from
these equations. What Mathieu used them for first was to generate recurrence
relations for their Taylor series expansion, which can be used about any
point. This can also be done for the original formulation, of course, but an
important difference is seen between the two forms: in the original
formulation, the Taylor series recurrence depends on all previously computed
terms; for the algebraic formulations, the recurrence relation depends only on
a finite number of the previous terms. We give the recurrence relation for the
original formulation in equation (5.2.1) in section 5, while the recurrence
relations for the algebraic formulations are already given in [58]. For
instance, for equation (15) if

$P(\nu)=\sum_{k\geq 0}\rho_{k}(\nu-\nu_{0})^{k}$

then (after the first few terms which have to be separately investigated),

$\displaystyle\rho_{n+4}=$
$\displaystyle\frac{1}{(\nu_{0}^{2}-1)(n+3)(n+4)}\Bigg{(}4q\rho_{n}+8\nu_{0}q\rho_{n+1}$
$\displaystyle+\left(4{\nu_{0}}^{2}q-{n}^{2}+a-4n-2q-4\right)\rho_{n+2}-\left(2n+5\right)\left(n+3\right)\nu_{0}\rho_{n+3}\Bigg{)}\>.$

We computed this recurrence relation automatically from equation (15) by using
the gfun package in Maple, specifically its diffeqtorec command [74]. Mathieu
gave the simpler form at $\nu_{0}=0$, and separated out the even and odd
series so that each recurrence relation involved only two prior terms. This
means that the Mathieu functions are what is now called $D$-finite or
holonomic, and can therefore be computed to high precision with an
asymptotically fast algorithm. See [8, 61, 62, 86].

Mathieu then considered properties of the functions that could be deduced from
these power series, which could also be interpreted as series in powers of
$\cos\alpha$ or of $\sin\alpha$. In particular, he used them to count real
roots.

####  3.1.5 Modified Mathieu Functions

In order to solve the vibrating drum problem, Mathieu had also to solve
equation (9). He chose to do this in a way analogous to the series solution
for Bessel’s equation that he gave in his introduction, and discussed how to
find the real roots thereof, which are necessary for matching the fixed
boundary condition at the elliptical rim of the membrane. In our terms, the
modified Mathieu functions are simply the Mathieu functions with purely
imaginary argument: $\mathrm{Ce}_{g}(x;q)=\mathrm{ce}_{g}(ix;q)$ and
$\mathrm{Se}_{g}(x;q)=-i\mathrm{se}_{g}(ix;q)$. We show two such functions in
figure 5.

We may now discuss the details of figure 3. We chose a drum shape with aspect
ratio $5:3$. We also chose to look at an even mode corresponding to
$a_{3}(q)$, so this means our pure tone will be described by
$\mathrm{Ce}_{3}(\beta;q)\mathrm{ce}_{3}(\alpha;q)$. We decided that it should
be, for the introduction, a simple picture with no elliptical nodal lines,
only hyperbolic; this meant that we were looking for the first zero of
$\mathrm{Ce}_{3}(\beta;q)$. Since the coordinates are
$x=c\cosh\beta\cos\alpha$, $y=c\sinh\beta\sin\alpha$ (alternatively,
$x+iy=\cos(\alpha-i\beta)$) we will want $c\cosh\beta=5$ and $c\sinh\beta=3$;
this gives $\beta=\ln 2$ and $c=4$. Now we want the value of $q$ so that
$\mathrm{Ce}_{3}(\ln 2;q)=0$. By zerofinding on $q$, we find that $q\approx
8.5676$. We used simple bisection, because we had not at that time implemented
differentiation with respect to $q$. The physics of the membrane would then
give the frequency of oscillation $\sin 2\lambda mt$ via $q=\lambda^{2}c^{2}$,
with the membrane parameter $m$. This value of $q$ gives the eigenvalue
$a_{3}\approx 14.6695$. The hyperbolic nodal lines are at approximately
$\alpha=\pm 0.9857$, $\pm\pi/2$, and $\pm 2.156$. The contours plotted were at
levels $\pm[0,11,22,33,44]/40$, remembering that our normalization is so that
$\mathrm{ce}_{3}(0;q)=1=\mathrm{Ce}_{3}(0;q)$. More details and more figures
can be found in [22].

(a) $\mathrm{Ce}_{0}(\beta;2)=\mathrm{ce}_{0}(i\beta;2)$

(b) $\mathrm{Se}_{1}(\beta;2)=\Im(\mathrm{se}_{1}(i\beta;2))$

Figure 5: (Left) A graph of $\mathrm{Ce}_{0}(q,\beta)$ when $q=2$. As the
argument $\beta$ increases, the function becomes increasingly oscillatory.
This value of $q$ could be used for an elliptical drum whose vertical
dimension was such as to coincide with a zero of this function (units
depending on the locations of the foci at $\pm c$.) (Right) A graph of
$\mathrm{Se}_{1}(\beta;q)=\Im(\mathrm{se}_{1}(i\beta;q))$ when $q=2$.

####  3.1.6 A short biography of Mathieu

Émile Léonard Mathieu (15 May 1835—19 October 1890) attended the École
Polytechnique de Paris, taking the entrance examination in 1854. He defended
his doctoral thesis in pure mathematics in 1859, before a committee consisting
of Lamé, Liouville, and Serret. He was well-regarded by the community at the
time for some of his work, and indeed is still known today for what are called
“Mathieu groups”. Apparently, though, he did not receive the positions that he
truly wanted; he then turned to applied mathematics (specifically,
Mathematical Physics) to see if that would “more engage the interest of
scientific men” [28]. In spite of this change, and in spite of winning a Gold
Medal in 1867, he was repeatedly passed over, and in 1869 he left for a
position at Besançon, becoming Chair of Pure Mathematics there in 1871. He
remained there until 1873 when he took up the Chair in Pure Mathematics at
Nancy, where he remained until his death in 1890. His obituary is a very
interesting read. In it Duhem [28] praises Mathieu’s achievements, calling him
the natural successor of Poisson, and blaming “fashion” and “politics” for
passing him over (thus suggesting that the fashion and politics in science and
mathematics were as alive and well then as they are today!).

His being passed over may have been a result of changing fashion, as Duhem
contends, or may simply have been an artifact of the Golden Age (for
mathematics) that he lived in. For instance, one of the positions that he
wanted, a Chair at Sorbonne, was awarded to Picard instead, who was Hermite’s
son-in-law and Mathieu thought this was a scandal; to be fair, it could be
argued that Mathieu’s record was superior at the time. But in other cases it
is clearer now. Mathieu complained that he came second to another favourite of
Hermite for another Chair, in this case to Hermite’s student Henri Poincaré.
The modern view must be different from Mathieu’s: It would be very difficult
today to imagine choosing Mathieu over Poincaré for any Chair.

We find in [16] (which contains an interesting view of the tension between
Paris and the provinces, and passages from Mathieu’s correspondence) still
other reasons why Mathieu was perpetually not chosen, and it seems that the
judgement of the Establishment that others were more worthy was, in the end,
justifiable. Even more, it was a turbulent time in France generally: the coup
when Napoleon III took power happened in 1851, and the Franco-Prussian war in
which Napoleon III was captured ended when the Prussians took Paris in 1870,
just as one example of how the larger world may have intruded on academic
life. It is quite believable that in these turbulent social circumstances many
deserving people did not receive all the recognition that they had earned.

In spite of all the difficulties of the times, however, Mathieu left a very
significant body of mathematical advances for posterity. We have surveyed only
a small corner of his work in this present paper, and even from just this it
is clear today that he was one of the best mathematicians of his age, which
included some of the greatest ever known.

###  3.2 Sir Edmund Taylor Whittaker (1873–1956)

####  3.2.1 Whittaker & Watson on Mathieu functions

It was Whittaker, another of the giants of nineteenth- and early twentieth-
century mathematics, who bestowed the name the Mathieu equation on equation
(2) and the name the Mathieu functions on the even and odd periodic solutions
of the Mathieu equation, “and these only” [90]. According to Whittaker’s
obituary [80], the name was given in the paper in the fifth ICM [88].
Attributing a person’s name to an equation or a function is a significant
event in mathematics because people (even mathematicians) are social animals,
and we simply pay more attention when a person’s name is involved. Such
namings often get it wrong, of course: “Stigler’s Law of Eponymy” states that
no scientific discovery is named after its original discoverer. For instance,
Puiseux series are named after Victor Puiseux (1820–1883) but were in fact
discovered by Newton, and for another, Young’s modulus—written about in 1807
by Thomas Young—was described 25 years prior to that by Ricatti, but in any
case was also discovered by Newton. In this case however we think Whittaker
got it right, and Mathieu deserves the credit. The more descriptive “elliptic
cylinder functions” is also used on occasion, and was also used by Whittaker.
Nowadays Mathieu functions and Modified Mathieu functions are also called
Angular Mathieu functions and Radial Mathieu functions, in accordance with
their roles in the elliptic coordinate system: $\beta$ is more like a radius,
and $\alpha$ more like an angle.

A full chapter of “A Course of Modern Analysis” by Whittaker and Watson [90]
is devoted entirely to Mathieu functions. Any serious study of these functions
should begin with that book. For some reason, the authors rescaled the
equation and have $16q$ there instead of the $2q$ (or $2h^{2}$) that Mathieu
had. This is of no real consequence. More interestingly, they use the same
normalization convention that Mathieu did in his perturbation series
computations: they take the coefficient of $\cos gz$ in the expansion of
$\mathrm{ce}_{g}(z)$ to be unity, and similarly the coefficient of $\sin gz$
in the expansion of $\mathrm{se}_{g}(z)$ to be unity. This differs from modern
practice.

Several important theorems are established in that chapter: the orthogonality
under the bilinear form of equation (11) is proved (by appeal to results from
an earlier chapter), several integral equations are established, and the
Floquet theory of the solutions to periodically-forced ODEs888Achille Marie
Gaston Floquet (15 December 1847, Épinal – 7 October 1920, Nancy) is applied
to the non-periodic solutions of the Mathieu equation, and more generally to
Hill’s equation. The Floquet theory shows that solutions must exist in the
form $\exp(\mu z)\phi(z)$ where $\phi(z)$ is periodic with period $\pi$
(because the periodic forcing of the Mathieu equation has that period) and
$\mu$ (actually, in modern works starting in [1], $\nu$ where $\mu=i\pi\nu$)
is called the characteristic exponent. Regions where $\Re(\mu)>0$ indicate
that the solution is unstable. The Mathieu equation is an important but
special example for Floquet theory, because the Mathieu equation is also even
and therefore $\exp(-\mu z)\phi(-z)$ is also a solution; this implies that for
both solutions to be stable, we must have $\Re(\mu)=0$; moreover the regions
in $(a,q)$ parameter space where the solutions are stable are bounded by
characteristic lines containing periodic solutions. For more information about
Floquet theory, an excellent introduction can be found in [81].

Below is one of the integral equations established in [90]: if $G(\eta)$ is an
even Mathieu function, then there is a characteristic number $\lambda$ for
which (translating from the $16q$ convention of Whittaker to the $2q$
convention used in this paper)

(17)
$G(\eta)=\lambda\int_{-\pi}^{\pi}e^{2\sqrt{q}\cos\eta\cos\theta}G(\theta)\,d\theta\>.$

We have not tried using this as a method of computing even Mathieu functions,
although Whittaker and Watson claim that it “affords a simple manner of
constructing” them.

Whittaker and Watson attribute to Lindemann the change of variable
$\zeta=\cos^{2}z$ (the DLMF has $\sin^{2}z$ but this is the same) which turns
the Mathieu equation into an algebraic differential equation or ADE, even
though Mathieu had done the same using $\nu=\cos z$ and $\mu=\sin z$, as
previously discussed. An ADE is, loosely speaking, an ordinary differential
equation with rational functions—usually polynomials—as coefficients; this is
not the same as a differential algebraic equation or DAE, which notion is more
familiar to numerical analysts. Whittaker and Watson make further references
to works of Abel, Stieltjes, Sylvester and others in the study that results
from this observation. Whittaker and Watson make an asymptotic connection of
Mathieu functions to Bessel functions using this form of the equation: in this
formulation it is more natural, but still a bit involved.

Whittaker and Watson also show that the Fourier series for Mathieu functions
are well-defined, at least for small enough $q$, by exhibiting convergent
power series for the coefficients. This marks an important step.

####  3.2.2 A short biography of Whittaker

We take the following material from the remarkably well-written and well-
informed Wikipedia article on Whittaker, supplemented by the mathematical
obituary written by G. Temple [80]. Sir Edmund Taylor Whittaker (24 October
1873–24 March 1956) studied mathematics and physics at Trinity College,
Cambridge, starting in 1892. He was elected Fellow of the College in 1896 and
continued there until 1906. He was elected Fellow of the Royal Society of
London in 1905. He then became Royal Astronomer of Ireland, and Andrews
Professor of Astronomy at Trinity College Dublin. In 1911 he joined the
University of Edinburgh and in 1912 was elected Fellow of the Royal Society of
Edinburgh, later becoming President. He was knighted by George VI in 1945 “for
service to mathematics”.

His mathematical obituary previously cited runs $22$ pages, plus a facing
portrait and a further four pages listing Whittaker’s complete works. It is a
mathematical biography well worth reading. The topic headings are 1.—Algebra
2.—Interpolation (exhibiting significant work in statistics and in numerical
analysis, which the anonymous Wikipedia author also takes particular note of)
3.—Automorphic functions 4.—Astronomy 5.—Potential theory and special
functions [It is here that we learn that Whittaker’s 1912 paper established
that the Mathieu equation is “the simplest linear differential equation which
is not reducible to a particular or degenerate case of the hypergeometric
differential equation”; and it is also here that we learn that E. L. Ince was
a research student working with Whittaker in Edinburgh and it was from this
period that Ince gained his interest in Mathieu functions.] 6.—Dynamics
7.—Relativity and Electromagnetic Theory 8.—Quantum Theory 9.—Scientific Books
and Monographs and 10.—Historical and philosophical writings. The very final
section, “11.—Influence” ends on a slightly ironic note: Temple claims that
Whittaker’s influence was felt not just by his works, but also by his coinage
of mathematical terms, some of which are listed in the final paragraph. Sadly,
of all those terms listed, namely ‘isometric circle’, ‘adelphic integral’,
‘cotabular functions’, ‘cardinal function’, ‘congruent hypergeometric
function’, ‘Mathieu function’, and ‘calamoids’, few apart from the Mathieu
functions are widely known today. We confess to curiosity as to what he meant
by a “calamoid”: the term seems to have survived only in botany, and has to do
with palm leaves.

But the mountain of scientific achievement that Whittaker built still stands
on its own.

###  3.3 Edward Lindsay Ince (1891–1941)

####  3.3.1 Fourier series recurrence relations

Mathieu’s perturbative solution of his equation and the resulting series in
powers of $\cos\alpha$ or $\sin\alpha$ already suggest that it is natural to
think of using Fourier series to represent these periodic functions. It seems,
however, that it was Ince who first made practical use of Fourier series for
Mathieu functions.

The basic idea is simple: we expand our proposed periodic solution $y(x)$ in
Fourier series:

(18) $\displaystyle y(x)=\sum_{k\geq 0}A_{k}\cos kx+\sum_{k\geq 1}B_{k}\sin
kx\>.$

Inserting this series into the Mathieu equation (2), using the multiplication
identities

$\displaystyle\cos(2x)\cos(kx)$
$\displaystyle=\tfrac{1}{2}\Bigl{(}\cos(k+2)x+\cos(k-2)x\Bigr{)}$
$\displaystyle\cos(2x)\sin(kx)$
$\displaystyle=\tfrac{1}{2}\Bigl{(}\sin(k+2)x+\sin(k-2)x\Bigr{)}\>,$

and then equating coefficients of $\cos kx$ and of $\sin kx$, gives us a
collection of recurrence relations. By circumstance (which Mathieu made
simpler by converting from $\cos^{2}x$ to the double-angle form), the $A_{k}$
coefficients only involve other $A_{k}$ coefficients, and moreover only those
that differ by $2$ in index; similarly for the $B_{k}$ coefficients. The edge
conditions (those recurrences specialize when $k=0$ and $k=1$ to slightly
different forms) produce a set of equations that can be written as follows:

(19) $\displaystyle aA_{0}\phantom{XXX}-q$ $\displaystyle A_{2}$
$\displaystyle=0$ (20) $\displaystyle-2qA_{0}+(a-4)$ $\displaystyle
A_{2}\phantom{XXX}-qA_{4}$ $\displaystyle=0$

and, for all $k\geq 2$,

(21) $\displaystyle-qA_{2k-2}+(a-(2k)^{2})A_{2k}-qA_{2k+2}=0\>.$

The odd cosine coefficients must instead satisfy

(22) $\displaystyle(a-1-q)A_{1}-q$ $\displaystyle A_{3}$ $\displaystyle=0$

and, for all $k\geq 1$,

(23) $\displaystyle-qA_{2k-1}+(a-(2k+1)^{2})A_{2k}-qA_{2k+3}=0\>.$

The even sine coefficients give

(24) $\displaystyle(a-4)$ $\displaystyle B_{2}-qB_{4}$ $\displaystyle=0$

and, for all $k\geq 2$,

(25) $\displaystyle-qB_{2k-2}+(a-(2k)^{2})B_{2k}-qB_{2k+2}=0\>.$

Finally, the odd sine coefficients give

(26) $\displaystyle(a-1+q)B_{1}-q$ $\displaystyle B_{3}$ $\displaystyle=0$

and, for all $k\geq 1$,

(27) $\displaystyle-qB_{2k-1}+(a-(2k+1)^{2})B_{2k+1}-qB_{2k+3}=0\>.$

All four of these sets of recurrence relations give rise to infinite
tridiagonal matrices, each one symmetric except the first. By an artifice,
namely replacing $A_{0}$ by $A_{0}\sqrt{2}$, we may symmetrize even the first
one.

The infinite eigenvalue problem then becomes

(28) $\left[\begin{array}[]{cccccc}0&\sqrt{2}q&0&0&0&\cdots\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr\sqrt{2}q&4&q&0&0&\cdots\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr 0&q&16&q&0&\cdots\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&q&36&q&\cdots\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&0&0&q&64&\ddots\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\vdots&\vdots&\vdots&\vdots&\ddots&\ddots\end{array}\right]\left[\begin{array}[]{c}\sqrt{2}A_{0}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr A_{2}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr A_{4}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr A_{6}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr A_{8}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\vdots\end{array}\right]=a\left[\begin{array}[]{c}\sqrt{2}A_{0}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr A_{2}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr A_{4}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr A_{6}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr A_{8}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\vdots\end{array}\right]\>.$

The eigenvalues of this matrix are denoted $a_{0}(q)$, $a_{2}(q)$, $a_{4}(q)$,
$\ldots$ and indeed for real $q$ these occur in increasing order:
$a_{0}(q)<a_{2}(q)<a_{4}(q)<\cdots$. For complex $q$ the situation is more
complicated, but the set of eigenvalues remains discrete and countable,
essentially because the diagonal of the matrix above contains entries that
grow rapidly enough with the row index.

The other three symmetric tridiagonal matrices are constructed analogously and
have analogous properties. For completeness, we include them here:

(29) $\left[\begin{array}[]{ccccc}1+q&q&0&0&\cdots\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr q&9&q&0&\cdots\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&q&25&q&\cdots\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&q&49&\ddots\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\vdots&\vdots&\vdots&\ddots&\ddots\end{array}\right]\left[\begin{array}[]{c}A_{1}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr A_{3}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr A_{5}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr A_{7}\\\ \vskip 6.0pt
plus 2.0pt minus
2.0pt\cr\vdots\end{array}\right]=a\left[\begin{array}[]{c}A_{1}\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr A_{3}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr A_{5}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr A_{7}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr\vdots\end{array}\right]$

The eigenvalues of equation (29) are denoted $a_{2k+1}(q)$.

(30) $\left[\begin{array}[]{ccccc}4&q&0&0&\cdots\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr q&16&q&0&\cdots\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&q&36&q&\ddots\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&q&64&\ddots\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\vdots&\vdots&\vdots&\ddots&\ddots\end{array}\right]\left[\begin{array}[]{c}B_{2}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{4}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr B_{6}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{8}\\\ \vskip 6.0pt
plus 2.0pt minus
2.0pt\cr\vdots\end{array}\right]=a\left[\begin{array}[]{c}B_{2}\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr B_{4}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr B_{6}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{8}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr\vdots\end{array}\right]$

The eigenvalues of equation (30) are denoted $b_{2k}(q)$.

(31) $\left[\begin{array}[]{ccccc}1-q&q&0&0&\cdots\\\ \vskip 6.0pt plus 2.0pt
minus 2.0pt\cr q&9&q&0&\cdots\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
0&q&25&q&\cdots\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&0&q&49&\ddots\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\vdots&\vdots&\vdots&\ddots&\ddots\end{array}\right]\left[\begin{array}[]{c}B_{1}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{3}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr B_{5}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{7}\\\ \vskip 6.0pt
plus 2.0pt minus
2.0pt\cr\vdots\end{array}\right]=a\left[\begin{array}[]{c}B_{1}\\\ \vskip
6.0pt plus 2.0pt minus 2.0pt\cr B_{3}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr B_{5}\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr B_{7}\\\ \vskip 6.0pt
plus 2.0pt minus 2.0pt\cr\vdots\end{array}\right]$

The eigenvalues of equation (31) are denoted $b_{2k+1}(q)$.

Remark. If $q$ is real, then these are real symmetric eigenvalue problems,
which have special properties connected to orthogonal polynomials [91]. If $q$
is complex, then the matrices are complex symmetric, not Hermitian. This has
some important consequences that we will pursue in section 4.

Since matrices are nearly universally familiar nowadays, it may be simplest to
think of these infinite complex symmetric tridiagonal matrix eigenproblems as
determining the eigenvalues and corresponding eigenfunctions of the Mathieu
equation. Truncating these infinite matrices causes no problems, as previously
stated. This, or at least in terms of determinants, seems to have been the way
that Ince thought of the problem [49], although he did not have direct
numerical methods to solve the eigenproblem and so he went on to solve the
problems perturbatively, recovering and extending Mathieu’s power series in
$q$ for the first few eigenvalues and for the coefficients $A_{k}$ and $B_{k}$
of the corresponding eigenfunctions. These series and their computation remain
of interest, with (to quote Blanch [13]) “an algorithm for generating the
successive [series coefficients], suitable for computers” being published in
[73] 999This last paper is interesting: the author, a researcher at “TRG
Incorporated in Melville, NY” appeared to be disappointed at their main
result, which was a set of recurrence relations which could be used
numerically to generate the numerical series coefficients of both
eigenfunctions and eigenvalues: it seemed that they were really searching for
a good test problem for early symbolic computation programs! We were
intrigued, but unable to find out much about TRG Incorporated except that
staff there carried out research on computational fluid simulations and lasers
in the 1960’s. As a further remark on this quotation, we are not sure if
Blanch is referring to human computers or to machines, although we incline to
the latter interpretation because she does use two kinds of IBM computers for
that article.. Subsequent papers give similar recurrence relations, and, if
the order is large enough and not too many terms are needed, explicit symbolic
formulae for the perturbative coefficients. Such a formula is termed “generic”
in [34].

Ince also derived a continued fraction101010The authors of [20] attribute the
derivation of the continued fraction form to Heine in [44] but it seems likely
that Ince’s derivation was independent. from the above recurrence relations.
We will continue to use tridiagonal matrices for the moment, as is done in
many numerical methods today such as [93]. The paper [20] shows especially
good results, and the authors advocate strongly for the conceptual and
computational use of these matrices.

When $q=0$ the matrix eigenvalues are simply the squares of the whole numbers:
$0$, $1$, $4$, $9$, and so on. They are, technically, double eigenvalues
(taking even and odd functions together), but in this case they retain two
independent eigenfunctions $\cos kx$ and $\sin kx$, except if $k=0$. For real
$q$ the eigenvalues can be sorted in increasing order. See Figure 6.

The evenness of the Mathieu equation, and its invariance under the
transformation $z\to z\pm\pi/2$ and $q\to-q$, mean that the eigenvalues have
the following symmetries: $a_{2n}(-q)=a_{2n}(q)$, $a_{2n+1}(q)=b_{2n+1}(q)$,
and $b_{2n+2}(-q)=b_{2n+2}(q)$. These are equations 28.2.26–28.2.28 in the
DLMF. There is also the conjugate symmetry $a(\overline{q})=\overline{a(q)}$
and similarly for $b(q)$.

Figure 6: The first few eigenvalues $a_{g}(q)$ and $b_{g}(q)$, that is, the
eigenvalues of the infinite tridiagonal matrices in equations (28)–(31), as
functions of the real variable $q$. This graph is modelled after the similar
one in Fig 28.1 of the DLMF.

Obviously if $y(x)$ is an eigenfunction corresponding to a given eigenvalue
$a_{k}(q)$ or $b_{k}(q)$, then so is any multiple of $y(x)$. This can be seen
in the eigenvalue/eigenvector formulation by multiplying the vector of Fourier
coefficients by any nonzero constant. We therefore need to choose a
normalization in order to choose a definite eigenfunction. Sadly, this is done
in several different ways in the literature. When using a particular software
package, one has to pay attention to the choices that the authors have made.

The eigenvalues of truncations of these infinite matrices converge quite
quickly to the desired eigenvalues of the Mathieu equation although the rate
of convergence depends on the eigenvalue. See figure 7. There are convergent
series containing eigenvalues that can be used as numerical checks on the
accuracy, but for this graph we simply computed high-precision values of the
eigenvalues by a continued fraction method and compared the results to the
truncated matrix eigenvalues. As stated, the matrices are real symmetric
tridiagonal for real $q$ (or can be made to be). As such, there are fast
algorithms for their computation, some based on Rayleigh iteration [69]; in
LAPACK there are SSTEV and DSTEV; see the Netlib repository. If $q$ is
complex, then the matrices are complex symmetric tridiagonal, which are harder
to solve, but for which there are still interesting algorithms [4]. There are
specialized algorithms to find multiple eigenvalues in the infinite dimension
case [63].

We have not felt the need to resort to fast methods: an ordinary slow
($O(n^{3})$) QR eigenvalue algorithm is perfectly fine, because the matrix
dimensions are so small in modern terms. After all, the Fourier series
converges spectrally to the Mathieu equation [77], and thus relatively few
eigenvalues are needed. However, in an application where there were very many
modes needed for many different values of $q$, fast methods would be quite
valuable.

The connection to orthogonal polynomials is in the case of complex $q$ more
strained, and may instead be more akin to the skew-symmetric eigenvalue
problem: see [51] for some surprises there.

(a) Convergence when $q=2$

(b) Convergence when $q=15+4\,i$

Figure 7: What dimension matrix is needed to get an accurate value of
eigenvalue $a_{2k}$? Here we plot the errors in $a_{2k}(q)$ for $3\leq k\leq
30$ by using a matrix of dimension $k+\Delta_{n}$ where $\Delta_{n}$ is as
indicated on the figure: either $1\leq\Delta_{n}\leq 6$ as on the left, or
twice that as on the right. That is, for $a_{6}(2)$ (top symbol in each column
of points in the left hand graph) we plot the errors in using matrices of
dimension $3+1:6=4$, $5$, $6$, $7$, $8$, and $9$. This tells us that to get
double precision we need dimension $9$ (dimension $8$ almost works). In
contrast, for $a_{6}(15+4\,i)$, (top symbol in each column of points in the
right-hand graph) we consider matrices of dimension $3+2:2:12=5$, $7$,
$\ldots$, $15$. In this case we need dimension $13$ to get double precision
accuracy in this eigenvalue. Higher-order eigenvalues always need a matrix big
enough to contain the desired eigenvalue and all smaller ones; it may be
surprising to see that they get more accuracy from any extra dimensions than
the smaller ones do. When $k=30$ (lowest points in each column) the accuracy
is significantly better; the matrix dimension for the low point on the right
is $42$. All computations done in $60$ Digits.

####  3.3.2 Asymptotics

There is some asymptotic work for Mathieu functions in [90], and indeed even a
little in Mathieu’s original work; but the first serious studies of the
asymptotics of large eigenvalues was [49]. In this paper Ince first
recapitulates some of Mathieu’s zero-counting work (it is not clear that Ince
knew that Mathieu had done this as well) using, as Mathieu had, Sturm’s
theory. Next, Ince established what are now the first two terms of DLMF 28.8.1
(we print the first four from there, below): if we denote $h=\sqrt{q}$ and
$s=2m+1$, then as $h\to\infty$ with $m=0$, $1$, $2$, $\ldots$,

$\left.\begin{matrix}a_{m}\left(h^{2}\right)\\\
b_{m+1}\left(h^{2}\right)\end{matrix}\right\\}\sim-2h^{2}+2sh-\frac{1}{8}(s^{2}+1)-\frac{1}{2^{7}h}(s^{3}+3s)-\frac{1}{2^{12}h^{2}}(5s^{4}+34s^{2}+9)-\cdots\>.$

This was created as a purely real result, but works for at least some complex
values of $q$ as well. These have been implemented to arbitrary order in some
computer algebra systems, for instance Maple.

####  3.3.3 Nonperiodicity of the other independent solution

In [47], Ince proved that the Mathieu equation (and similarly the Hill
equation) “can admit but one solution of period $\pi$ or $2\pi$”. This is
generally taken to mean that the other linearly-independent solution cannot be
periodic, and indeed this is true for periods $\pi$ and $2\pi$. See figure 8
where one such secularly-growing solution is graphed. In [48] Ince established
for small $q$ that if, for instance, the periodic solution was
$\mathrm{ce}_{p}(z)$, then the second, necessarily non-periodic, solution
would be

$y=\mathrm{se}_{p}(z)+Kz\mathrm{ce}_{p}(z)+\phi(z)$

where both $K$ and $\phi(z)$ would be $O(q^{p})$ in size and $\phi(z)$ would
be periodic. As we can see in figure 8(a) this is true even if $q$ is not
small (in fact we took $q=1.4688\,i$ approximately, and its eigenvalue
$a\approx 2.0886$ for this figure; the periodic solution for this value is
plotted later, in figure 11(a)). A general theorem to this effect can be found
for instance in [30].

However, in a throwaway remark in [13], Gertrude Blanch states “In contrast to
this, it is known that there exist other eigenvalues that give rise to
solutions of period $k\pi$, where $k$ is an arbitrary integer greater than
$2$. For these eigenvalues, all solutions are periodic, if one is.” She gives
the reference [68] in Chapter 20 of [1]. This is also covered in [30].

(a) A non-periodic solution

(b) Unstable solutions

Figure 8: Some non-periodic solutions of the Mathieu equation. On the left, we
have the nonperiodic solution for the Mulholland-Goldstein double point
$q\approx 1.4688\,i$ showing secular growth as Ince proved. Black is the real
part of the solution, red is the imaginary part. On the right we show absolute
values of solutions for $a=4$ and $q=10$, well into an unstable region in
$(a,q)$ space by the Floquet theory, showing exponential growth. Indeed the
solutions are of the predicted form $\exp(\mu z)\phi(z)$ where $Re(\mu)>0$ and
$\phi(z)$ is periodic; we can see that this particular periodic function has
two zeros by the “spikes” in the graph.

####  3.3.4 A short biography of Ince

Edward Lindsay Ince (30 November 1891 – 16 March 1941) was, as previously
mentioned, a research student of Whittaker and had spent the years 1909–1913
at the University of Edinburgh. He had a short but colourful life and career,
which included what would have then been an exotic episode at the newly-
founded Egyptian University of Cairo111111Founded 1908, it is now Cairo
University. from 1926 to 1931. Returning to England for reasons of health and
family, he took various positions, ultimately becoming Head of Technical
Mathematics at the University of Edinburgh. He was elected a Fellow of the
Royal Society of Edinburgh in 1923. He died of leukemia in Edinburgh, aged
just 49. Just before his death he was awarded the Makdougall-Brisbane prize
for his work on Lamé functions, but he died before he learned of the award.

There are two obituaries of Ince that we know of: one by Whittaker [89], and
another by Aitken [2]. Both make interesting reading, although there is more
mathematics in the one by Whittaker. The one by Aitken is perhaps more
personal, although it does contain the passage “Ince firmly believed that
theoretical solutions of problems, however abstractly elegant, were incomplete
unless the mathematician either tabulated the solving functions himself or
rendered them tabulable. Perusal of his papers will show that in his chosen
field of research he achieved both of these objects”. As an aside, Aitken’s
biography is itself worth reading; in addition to his mathematical
accomplishments he was elected to the Royal Society of Literature for his
memoir “Gallipoli to the Somme: Reflections of a New Zealand Infantryman”,
having served in both battles. In contrast, Ince was not permitted active duty
in the First World War owing to ill-health, although he did do a term of
National Service.

Returning to mathematics and to Ince, Whittaker tells us that Ince started out
reading mathematics under Professor George Chrystal, who we principally know
nowadays as the person who wrote the text that inspired Ramanujan; after
Chrystal died in 1911, “Ince continued under his successor”, namely Whittaker
himself. Whittaker then carefully describes Ince’s achievements with the
Mathieu functions, including Ince’s proof [47] in 1922 that there could not be
two independent solutions with period $\pi$ or $2\pi$ to the same Mathieu
equation. Whittaker echoes Aitken’s remark about Ince’s sensibility with
regard to computation, and amplifies it: “Ince held that an important part of
a pure mathematician’s duty is to provide tables for the use of physicists and
astronomers, and he was well aware that the possibility of constructing such
tables without a colossal expenditure of time and energy depends on the
progress of theoretical analysis.” Whittaker remarks that Ince’s 1932 tables
of the Mathieu functions their zeros and turning points was “A splendid piece
of work, performed single-handed save for some help by an Egyptian assistant.”
Ince gave the name of his helper, when he acknowledged Mansy Shehata, who was
then an Assistant in Pure Mathematics at the Egyptian University in Cairo.
Ince also acknowledged grant support in purchasing calculating machines, which
seemed to be of significant use.

Ince’s work on the Mathieu functions was important in making computation of
them practical; Gertrude Blanch, as we shall see next, relied on many of his
results.

###  3.4 Gertrude Blanch (1897–1996)

####  3.4.1 Numerical computation via continued fractions

Blanch’s first publication on Mathieu functions was [10]. This was reprinted
in 1950 by the National Bureau of Standards “with the permission of the
editors of the Journal of Mathematics and Physics to meet a continuing
demand”. This refined and improved the continued fraction method that Ince
used and in particular allowed the error in the Fourier coefficients to be
estimated. She wrote a paper in the Transactions of the American Mathematical
Society on the asymptotics of the odd periodic Mathieu functions [11],
extending and correcting the works of others. Blanch wrote an influential
paper in SIAM Review on continued fractions, namely [12]. This paper contained
a detailed linearized rounding error analysis for continued fractions, and
used both Bessel functions and eigenvalues of Mathieu functions as examples.
She argued that continued fractions gave the preferred method for computing
Mathieu function eigenvalues in [13]. Then Blanch & Clemm [14], still later
improved by [45], went on to systematically compute the double
eigenvalues121212Blanch 1966 quotes using an “I.B.M. 1620 and 7094”. We merely
note the printed periods in the abbreviation “I.B.M.” used there, as opposed
to the simpler trademark IBM used nowadays, and refrain from comparing the
speed and memory of the early computers with what is available today (or what
is coming)..

Blanch’s use of continued fractions is different from the use of Gauss-type
continued fractions for evaluation of other special functions, such as those
of hypergeometric type. Those methods are explained in detail in, for example,
Gil et al. [38]; who also cite Blanch for her study of the numerical stability
of evaluation of continued fractions over the complex plane (as well as give
more modern results, of which there are many). But since the Mathieu equation
is not of hypergeometric type, continued fractions cannot be applied directly
in this fashion. Instead, they are used to define a class of nonlinear
equations whose roots are the eigenvalues of the Mathieu function; and then
the partial quotients give the elements of the corresponding eigenvector.

We believe that Blanch’s choice of this particular method of continued
fractions was at least partly because of the kind of computing she typically
did, which in her own words was “experimental in nature” and which by habit
were continually checked as the computation went on. [Our computations for
this paper were carried out in a similar style.] In the first part of her
career the computations were, like Ince’s, carried out by hand and by desk
calculating machine. A significant difference is that while Ince had only one
assistant when he was in Cairo, Blanch organized a group of as many as 450
assistants, and later used digital computers. But her computations were in
many ways what we might call “artesenal” with a great deal of human
involvement.

Today people generally prefer software packages that can run without care on
the part of the user: in Blanch’s words, in a “robot-like” manner [12]. Most
people want to call a subroutine and be certain that it would never give error
messages, would always be fast, and would always return the correct answer. If
one is writing general-purpose software for evaluation of the Mathieu
functions, therefore, one might not choose continued fractions (as we will
see) because their occasional instabilities in the complex plane are somewhat
variable; on the real line Blanch had analyzed this behaviour and given a
sound rule for routine computation, but not for general complex $q$. For our
purposes in this paper, however, the continued fraction algorithm is perfect:
when carefully monitored, it is fast, flexible, generalizable, and any
instability can be controlled simply by using extra precision.

Blanch’s variant of the continued fraction algorithm can be described as
follows. This treatment looks a little different than the standard connection
of three-term recurrence relations, tridiagonal matrices, and orthogonal
polynomials [91], but in essence it is similar. We suppose first that $a$ and
$q$ are given (later we will perform iterations, looking for values of $a$ and
sometimes also of $q$ that make the continued fraction equal to zero). We will
need the auxiliary quantities

(32) $V_{m}=\frac{a-m^{2}}{q}$

for $m=0$, $1$, $2$, $\ldots$. These are not defined if $q=0$, but if $q=0$
then we already know that the eigenvalues are $g^{2}$ for integers $g$.
Equation (32) is equation (1.05) in [13]. We use this notation in the
recurrence relations (19)–(27), and further with the auxiliary quantities
$c_{m}$ where $c_{0}=2$ and all other $c_{k}=1$, in order to look after the
first edge case. We rewrite the recurrence relations (apart from the base
cases) as

(33) $A_{m+2}+c_{m-2}A_{m-2}-V_{m}A_{m}=0\>.$

Putting $G_{m}=A_{m}/A_{m-2}$ and dividing equation (33) by $A_{m}$ (Blanch
was very careful about what happened when any $A_{m}$ was zero, but here we
rely on IEEE arithmetic with signed zeros and infinities to get everything
right) we get

(34) $G_{m+2}+\frac{c_{m-2}}{G_{m}}-V_{m}=0\>.$

This recurrence can be written either as

(35) $G_{m+2}=V_{m}-\frac{c_{m-2}}{G_{m}}$

or

(36) $G_{m}=\frac{c_{m-2}}{V_{m}-G_{m+2}}$

Of course, these recurrences must be started correctly by using the
appropriate equations (19)–(27). Running equation (36) until the denominator
$Q$ is so large that $1/|Q^{2}|\leq\varepsilon$ where $\varepsilon$ is our
tolerance, starting from some $M>2$ so that $c_{k}=1$ for all $k\geq M-2$,

(37) $G_{M}=\frac{1}{V_{M}-\frac{1}{V_{M+2}-\frac{1}{V_{M+4}-\ddots}}}$

and because the $V_{m}$s grow like $m^{2}$ this continued fraction converges
for all $a$ and all $q\neq 0$. Call the $G_{M}$ computed in this way
$G_{M,\mathrm{tail}}$.

Now we use equation (35) with increasing $m$ starting from our known edge
cases (depending on which class of Mathieu eigenvalue we wish to compute), and
if and only if $a$ is an eigenvalue then the two values of $G_{M}$ in the
middle will agree. Call the $G_{M}$ computed in this way
$G_{M,\mathrm{head}}$. Let

(38) $T(a,q)=G_{M,\mathrm{head}}(a,q)-G_{M,\mathrm{tail}}(a,q)\>.$

$T(a,q)$ must be zero for $a$ to be an eigenvalue corresponding to $q$. The
edge cases for $G_{M,\mathrm{head}}$ determine whether this is an $a_{2k}$,
$a_{2k+1}$, $b_{2k}$, or $b_{2k+1}$ eigenvalue. Just which integer $k$
depends, as a rule, on whether there is an unambiguous continuous path in $q$
back to the eigenvalue with that index when $q=0$. Blanch gave a rule, as
previously stated, for choosing the $M$ in the middle so as to minimise
numerical instability for real $q$.

There are many methods one could use to find zeros of equation (38), but since
differentiation of $T(a,q)$ with respect to $a$ is simply a matter of
differentiating equations (35) and (36) (convergence is uniform in compact
neighbourhoods and so differentiation is permissible), we choose Newton’s
method. For a given $q$ and starting with an initial estimate $a^{(0)}$ for
the eigenvalue131313In order to find a good initial estimate, people usually
use continuation from $q$ near $0$: one solves for $q=q_{n}$, and then uses
that eigenvalue as an initial estimate for the eigenvalue at $q=q_{n}+\Delta
q$., the iteration is

$a^{(k+1)}=a^{(k)}-\frac{T(a^{(k)},q)}{T_{a}(a^{(k)},q)}\>.$

Blanch used eigenvalues with slightly different values of $q$ as initial
estimates for the eigenvalues she required, and this worked very well, unless,
of course, the eigenvalue was not simple (i.e. did not have multiplicity $1$),
in which case even more derivatives of $T$ turn out to be useful.

The eigenvalues of the Mathieu equation are as previously stated almost always
simple, but for isolated complex values of $q$ may have multiplicity $2$. In
particular, if $q=is$ where $s$ is real and $i^{2}=-1$, that is if $q$ is
purely imaginary, then as $s$ increases from zero we will necessarily
encounter double points: first at approximately $s=1.4688$ and then at
approximately $s=6.9289$ (see section 28.7 of the DLMF). The double point near
$s=1.4688$ was first studied in [65] and we will refer to it as the
Mulholland-Goldstein double point (we use it as an example, frequently). There
are a countable infinity of double points. It is proved in [60] that there are
no triple points of the Mathieu equation, or simultaneous double points—that
is, a value of $q$ for which two (or more) pairs of eigenvalues merge.

One might be tempted to dismiss double points because they occur “with
probability zero” if one chooses the parameter $q$ “at random”. But in
applications requiring complex $q$ the parameter will not usually be chosen at
random, and indeed the problem being modelled may call for a continuum of
values of the parameter. In the application that motivated us to study the
Mathieu functions, we were interested in all imaginary values of $q$, and this
necessarily included some double points. In some sense this means that the
“probability zero” events actually occur with “probability one”. This reversal
of expectation is a common occurrence in bifurcation phenomena, for instance
[40].

####  3.4.2 Double points

Blanch and Clemm used a variation on Newton’s method in their systematic
search for double points [14]. Their method was not as efficient as two-
dimensional Newton iteration, but it worked and it was the first method to do
so. Instead of just using one derivative with respect to $a$, they used two,
and expanding

$T(a,q)=T(a^{(k)},q)+T_{a}(a^{(k)},q)(a-a^{(k)})+\frac{1}{2}T_{aa}(a^{(k)},q)(a-a^{(k)})^{2}+\cdots$

they set this to zero and solved the resulting quadratic for the update to
$a^{(k)}$, choosing the smaller magnitude square root. This gives a
superlinearly-converging iteration for double roots (even faster for simple
roots), and together with interpolation enabled them to compute the smallest
$72$ double eigenvalues to what we would now call double precision. We
confirmed their work and recomputed all their roots. We used their tabulated
values as initial estimates, and have plotted the results in figure 9. Blanch
and Clemm carried out their computations on digital computers (we believe
models IBM 1620 and IBM 7024).

####  3.4.3 A short biography of Blanch

As previously stated, Gertrude Blanch wrote the chapter on Mathieu functions
in the classic handbook [1]. We have cited many of her papers on the Mathieu
functions and related matters. She apparently wrote in 1943 what might be the
first modern textbook on numerical analysis, and an updated version in 1982,
according to [41]. Unfortunately, we have not seen a copy of either edition.

We draw material for this section from several sources, including a transcript
of an interview with her in 1973 [83], an extensive biography by Grier [41],
and the collection of her papers at the Charles Babbage Institute at the
University of Minnesota. Gertrude Blanch’s extremely interesting life was
well-documented. We have in this short biography concentrated on mathematical
aspects of her life and left out important non-mathematical aspects. The above
resources are well worth consulting for a fuller picture.

Gertrude Blanch was born Gittel Kaimowitz in 1898 in Kolno, then part of
Russia. She came to the US in 1907 and attended high school in Brooklyn,
graduating in 1914. She changed her name, after her father died, to an
Anglicized version of her mother’s family name, Blanc. She became an American
citizen in 1921. She worked for fourteen years to get enough money to attend
university; her employer paid her tuition for New York University, where she
graduated summa cum laude in 1932. Apparently following the advice of one of
her professors there, Fay Farnum, she then went to Cornell, receiving her PhD
in algebraic geometry in 1935 under the guidance of Virgil Snyder and Wallie
A. Hurwitz141414Blanch credits them both, but the Mathematics Genealogy
Project only reports Snyder as her advisor. Snyder also seems to have advised
Farnum and ten (by our count) other women for their PhDs, the earliest—Anna
van Benschoten—in 1908.. After a short stint teaching at Hunter College for
someone on sabbatical leave, she took an office clerical administration and
accounting job. This administrative job was to prove important for her later
work with the Mathematical Tables Project, as detailed in [41]. In order to
keep her mathematical interests alive, she took a night course in relativity
at Washington Square College given by Arnold Lowan. When Lowan was asked to
create the Mathematical Tables Project under the New Deal Works Progress
Administration, he asked Blanch to join. She became Technical Director,
eventually organizing a group of 450 (human) computers. According to Grier,
she deserves much of the credit for the success of the project, and part of
that credit is due to her prior business-oriented administrative experience.
She published several papers during this time, including one with Hans Bethe
[15].

“During her time at the Mathematical Tables Project, she particularly enjoyed
working with the Mathieu functions, and these functions would be central to
the rest of her career.” [41, p.23].

When asked how she first got interested in Mathieu functions, she responded:
“Morse151515Philip Morse, at MIT, the first author of [64]. This monumental
work also describes Mathieu functions. was interested, for example, in Mathieu
functions. I got started on Mathieu functions because Morse needed them and
there were any number of small things and some special integrals that he came
across within his field.”

She remained interested in Mathieu functions for the rest of her career and
indeed after her retirement. Her ultimate academic appointment was as Head of
Mathematical Research in the Aerospace Research Laboratory at Wright Paterson
Air Force Base in Dayton, Ohio, where she wrote many of her papers. She became
a Fellow of the American Association for Advancement in Science in 1962. She
received the Federal Women’s Award from President Lyndon Johnson in 1964. She
retired in 1967, and died in 1997, aged 99.

###  3.5 Selected other works by other contributors

There is an excellent concise treatment of Mathieu functions, including the
Floquet theory, in Chapter XVI of Volume 3 of [30], “The Bateman Manuscript”.
In particular, it is there that we learn two interesting facts about
$\mu=i\pi\nu$, the Floquet characteristic exponent for the Mathieu equation:
first, that Poincaré found a way to compute it from the two basic solutions
$w_{I}(z)$ and $w_{II}(z)$ using the evenness of the Mathieu equation: at
least one of

$w_{I}(z)=\frac{e^{\mu z}\phi(z)+e^{-\mu z}\phi(-z)}{2\phi(0)}$

or

$w_{II}(z)=\frac{e^{\mu z}\phi(z)-e^{-\mu
z}\phi(-z)}{2(\phi^{\prime}(0)+\mu\phi(0))}$

will have nonzero denominator; now differentiate $w_{II}$ and compute
$w_{I}(\pi)$ and $w_{II}^{\prime}(\pi)$. Since $\phi(0)=\phi(\pm\pi)$ and
$\phi^{\prime}(0)=\phi^{\prime}(\pm\pi)$, we have that

(39) $\cosh\pi\mu=w_{I}(\pi)=w_{II}^{\prime}(\pi)\>.$

In the DLMF this equation (using cosine and not hyperbolic cosine, because
$\nu$ is used instead, where $\mu=i\nu$) is called the characteristic
equation, number 28.2.16. Blanch uses another convention, namely
$\mu=i\pi\nu$, defining $\nu$ differently, in [1]. The second interesting fact
is that the periodic solutions, that is the Mathieu functions, correspond to
the case $\mu=i\,n$ where $n$ is an integer; if $n$ is an even integer the
solution is periodic with period $\pi$ and if $n$ is odd the solution is
periodic with period $2\pi$. The case when $\mu=i\,r$ where $r$ is a rational
integer generates those mysterious solutions alluded to earlier which are
periodic with greater periods. See section 20.3 in [1] for more discussion of
this fact.

The earlier book [79] by M. J. O. Strutt was critically reviewed in [39],
which states “It is a useful book, both for pure mathematicians interested in
the theory of special functions, and for applied mathematicians compelled to
use the functions in their researches. It is worthy of consultation by both
classes, but it is rather superficial$\ldots$”.

The most extensive and thorough treatment of Mathieu functions is that of
Meixner, Schäfke, and Wolf [60], which is based on an earlier book by Meixner
and Schäfke that we have been unable to get hold of. The later edition handles
significant generalizations of the theory of Mathieu functions, and as
previously stated deals completely with the theory of double points. Nearly
every reference after this work cites it; indeed the Bateman Manuscript itself
cites (the 1950 book) as “forthcoming”. The additional author, G. Wolf, of the
second edition is the author of the DLMF Chapter 28 on Mathieu functions,
which updates Blanch’s Chapter 20 of [1].

There are many tables of the Mathieu functions. We believe that Blanch’s work
remains the most extensive, but we can also mention Bickley [9]. There are
tables of integrals and series for Mathieu functions in Prudnikov et al. [70].
The value of the numerical tables nowadays is of course lessened by the
availability of software; the tables of integrals and series may also have
been supplanted now by computer algebra systems, at least to some extent.

There are a large number of relatively modern books, papers, and software
packages for the computation of Mathieu functions. Analytical work on the
Mathieu function is extensive; to cite a paper almost at random, consider
[66]. Some of the works ignore most of the others, and hence rediscover
things, but others are very insightful. We particularly recommend [20] and
[43] for exceptional visualizations and clarity of exposition. On the purely
computational side, there are three sets of ACM Transactions on Mathematical
Software papers discussing implementations, the latest being [31]. There are
Python (scipy) and third-party Matlab implementations of (real) Mathieu
functions. The Mathematica implementation of complex Mathieu functions may be
very good (unfortunately, we have only had limited access to Mathematica so we
cannot be sure, but it may even have facilities for computing double
eigenvalues) and is described in [84]. The Maple implementation, with which we
are most familiar, encodes a substantial amount of analytical work including
both $q$-series and asymptotic series.

##  4 Double points

As previously noted, at certain isolated points in the complex $q$-plane, for
instance at the Mulholland-Goldstein point $q^{*}\approx 1.468768613785142\,i$
(reporting $16$ digits161616This double point was the first found: studied in
[65] and later computed by [18] to $3$ digits and then to double precision in
[14].), we have a double eigenvalue: $a_{0}=a_{2}\approx 2.088698902749695$.
Several interesting things happen at double points. First, the eigenfunctions
coalesce: here, $\mathrm{ce}_{0}(q^{*},\alpha)=\mathrm{ce}_{2}(q^{*},\alpha)$,
leaving a gap in the completeness of the set of eigenfunctions. This means
that we will have to supplement the set of eigenfunctions in some way in order
to expand arbitrary functions in series containing Mathieu functions.

Second, near to these double points, the ordering of the eigenvalues becomes
ambiguous. For $\Im(q)<\Im(q^{*})$ we have
$a_{0}(q)<a_{2}(q)<a_{4}(q)<\cdots$, but at $q^{*}$ equality occurs. For
$\Im(q)>\Im(q^{*})$, both $a_{0}(q)$ and $a_{2}(q)$ are complex, and ordering
is a matter of convention. The DLMF adopts the convention in this case that
$a_{0}(q)$ continues as $\Im(q)$ increases by choosing the branch with
negative imaginary part, while $a_{2}(q)$ takes the conjugate. See the
visualization in section 28.7 of the DLMF: they have paths for $a_{0}(q)$ and
$a_{2}(q)$ coming in and merging at the double point, and then emerging at
right angles: which way the paths are numbered on emergence is a convention.

The convention given in the DLMF makes sense provided one thinks of $q$
varying along the imaginary axis, and increasing. Approaching $q^{*}$ along
some other path in the $q$-plane may cause puzzlement: the ordering is
conventional, no more.

The discussion in [45] links the numbering of eigenvalues to the paths used
for continuation in $q$, and they describe branch cuts ending at each double
point which can be used to disallow paths that would disrupt the conventional
numbering. This analysis becomes more involved the larger the double points
get.

Also as $q$ approaches $q^{*}$, the coefficients in the Mathieu series
expansion for a given function will usually become singular. For example,
consider

(40) $\cos
2\alpha=c_{0}(q)\mathrm{ce}_{0}(q,\alpha)+c_{2}(q)\mathrm{ce}_{2}(q,\alpha)+\cdots\>.$

By numerical experimentation at high precision171717We are slightly
embarrassed to admit to how many figures we took these computations to: we
worked at $250$ decimal Digits in Maple, and solved the Mathieu equation with
a tolerance of $10^{-120}$. We then worked with $5$, $8$, $13$, $21$, $34$,
and $55$ Digit truncations of $q^{*}$ and calculated the corresponding $a_{0}$
and $a_{2}$ to $120$ Digit accuracy; this enabled us to identify the constants
in this section to $50$ Digits or more. We only report the double precision
values. Later, we confirmed these by simply computing the Puiseux expansion of
this point. we find that

(41) $\displaystyle c_{0}+c_{2}$
$\displaystyle=1.009185957186356-0.1210349964877181\,i$ (42) $\displaystyle
c_{0}-c_{2}$
$\displaystyle=\frac{-1.023431886611575+0.2551095295356106\,i}{\sqrt{q-q^{*}}}$

when we use the normalization convention that $\mathrm{ce}_{2k}(q,0)=1$, and
moreover that

(43) $\displaystyle a_{0}=$ $\displaystyle
a^{*}+d\cdot(q-q^{*})^{1/2}+O(q-q^{*})$ (44) $\displaystyle a_{2}=$
$\displaystyle a^{*}-d\cdot(q-q^{*})^{1/2}+O(q-q^{*})$

where $d\approx
1.659487804320256+1.659487804320256\,i=1.659487804320256(1+i)$. These values
were found by using orthogonality, which holds if $q\neq q^{*}$, and by high-
precision computation of

$\int_{x=0}^{2\pi}\mathrm{ce}_{2k}^{2}(q,x)\,dx=O(q-q^{*})^{1/2}\>.$

Then if $y^{*}=\mathrm{ce}_{0,2}(q^{*},\alpha)$ is the coalesced solution to
Mathieu’s equation, we can examine the most important contribution to the
perturbation by considering the solution to the perturbed equation

(45) $y^{\prime\prime}+(a^{*}\pm d\cdot(q-q^{*})^{1/2}-2q\cos(2\alpha))y=0\>.$

If $y=y^{*}+u(q^{*},\alpha)\cdot d\cdot(q-q^{*})^{1/2}+\cdots$, then a short
calculation dropping terms of $O(q-q^{*})$ and higher gives

(46) $u^{\prime\prime}+(a^{*}-2q^{*}\cos 2\alpha)u+y^{*}=0\>.$

The function $u$ must be periodic with the same period as $y^{*}$. It can
therefore be computed numerically alongside $y^{*}$ (by solving a boundary-
value problem) or alternatively can be expressed as an integral of Mathieu
functions against $y^{*}$. This argument is extended and formalized in a short
section in [60] starting on p. $82$. Here and with just this simple example,
we see that $u$ is essentially $\partial y/\partial a$ (found by solving a
variational equation). By combining the equations above, we can see that

(47) $\cos
2\alpha=(c_{0}+c_{1})\mathrm{ce}_{0,2}^{*}(q^{*},\alpha)+d(c_{0}-c_{2})\sqrt{q-q^{*}}\cdot
u(q^{*},\alpha)+O(q-q^{*})^{1/2}\>.$

This arrangement is continuous as $q\to q^{*}$ because
$(c_{0}-c_{2})\sqrt{q-q^{*}}$ is $O(1)$ in that limit. This shows explicitly
that by adding $\partial\mathrm{ce}_{0,2}(q^{*},\alpha)/\partial a$ to the
collection of Mathieu functions we are expanding with we ensure completeness
of the set and the possibility of the expansion. We know of no freely-
available software package for Mathieu functions that provides for the
computation of these extra functions. We discuss methods to compute
$u(\alpha)$ in section 5.4.

###  4.1 Computing double eigenvalues

We here discuss an effective method for computing the double eigenvalues. When
computing $T(a,q)$, we will need to also compute the derivatives $T_{a}(a,q)$,
$T_{q}(a,q)$, $T_{aa}(a,q)$ and $T_{aq}(a,q)$. Then we will be able to carry
out a two-dimensional Newton iteration for solving the two equations
$T(a,q)=0$ and $T_{a}(a,q)=0$ simultaneously. The equations for the iteration
looks like this:

(48) $\begin{bmatrix}T_{a}&T_{q}\\\
T_{aa}&T_{aq}\end{bmatrix}\begin{bmatrix}\Delta a\\\ \Delta
q\end{bmatrix}=\begin{bmatrix}-T\\\ -T_{a}\end{bmatrix}$

where all function evaluations and derivative evaluations occur at the current
estimates $(a^{(k)},q^{(k)})$. Then as usual $a^{(k+1)}=a^{(k)}+\Delta a$ and
$q^{(k+1)}=q^{(k)}+\Delta q$. Given sufficiently good initial estimates, this
iteration converges quadratically to double points $(a^{*},q^{*})$. This
method was first used by [45], and as previously noted is different from the
method of Blanch and Clemm. Initial estimates are typically obtained by
numerical continuation in $q$.

######  Remark 4.1.

Even if one has $q^{*}$ to double precision, one cannot naively compute
$a^{*}$ to the same precision: after all, being Hölder continuous but not
Lipschitz continuous at the double point it will be sensitive to errors on the
order $\sqrt{q-q^{*}}$. In ordinary numerical parlance, the double eigenvalue
is infinitely ill-conditioned, although the Hölder continuity constrains that
to some extent. With only double precision computation, this means one can
only expect to know $a^{*}$ to about single precision, and the double
eigenvalue will ordinarily have spuriously split into two simple (but close)
eigenvalues. Their mean, as is well-known, will be a double-precision
approximation to the true double eigenvalue (we will see that in adding the
two Puiseux series expansions the leading error terms cancel). This numerical
split, however, might be used to advantage in computation of independent
eigenfunctions: the norms of these erroneous approximate eigenfunctions will
be $O(q-q^{*})^{1/2}$ and therefore the associated Fourier coefficients in the
eigenfunctions will be large, but perhaps this is tolerable. The resulting
spectral expansion of the function may well be accurate enough for one’s
purposes. Perhaps this is the real reason no-one has developed the tools to
deal precisely with these extra eigenfunctions in practice; of course in
theory this has been understood since at least [60].

###  4.2 Initial estimates by continuation

As previously mentioned, to compute $a_{g}(q)$ or $b_{g}(q)$ by Newton’s
method, one typically needs an initial estimate for the eigenvalue. The
standard way to do this is to choose a path in the complex $q$ plane from
$q=0$ to the desired $q$, and then one knows $a_{g}(0)=g^{2}$ and
$b_{g}(0)=g^{2}$; one increments $q$ by a small amount, and then uses the
previous value of $a_{g}(q)$ or $b_{g}(q)$ as the initial estimate for a
Newton iteration for $a_{g}(q+\Delta q)$. This works well except in the
neighbourhood of double points, at which it is sensible to switch to two-
dimensional Newton iteration, locate the double point carefully, and then step
from there in the direction of the desired $q$.

The asymptotic formulae of [45] can be used instead to directly estimate the
location of large double eigenvalues for subsequent refinement by two-
dimensional Newton iteration. We have not made a systematic study of these
formulae, although we have tested using just their simplest approximations
from equations 4.14a–4.14d in that paper which give essentially an empirical
estimate of the approximate values of $q$ at double eigenvalues; other
equations in the paper give methods for numerically computing asymptotically
correct values of the double eigenvalues, based on WKBJ approximations and the
elliptic integral

(49) $I=\int_{0}^{\pi/2}\sqrt{a-2q\cos 2\alpha}\,d\alpha\>.$

This must be equal to $m\pi/2$ if $a=a_{m}$ or $a=b_{m}$ is an eigenvalue.
They dryly point out that “no simple expression for this integral is
available” and indeed Maple’s current exact expression for this value contains
$7,287$ characters and takes up more than two screens. We have not
investigated if that can be usefully simplified. Hunter and Guerrieri use
numerical evaluation of that integral together with rootfinding in order to
compute asymptotic estimates of the location of double points.

###  4.3 Computing local series for the eigenvalues by Newton’s method

We may use Newton’s method not just to compute a simple eigenvalue, or two-
dimensional Newton’s method to compute a double eigenvalue, but we can also
use (one-dimensional) Newton’s method to compute a series for an eigenvalue by
working in the ring of formal power series. In the case of a simple eigenvalue
at, say, $q=q_{s}$, we may compute a power series in $(q-q_{s})$ for the
eigenvalue $a_{g}(q)$, and simultaneously if we wish for the associated
eigenfunction. If instead we want a series expansion around a double
eigenvalue at $q=q^{*}$, then we may compute a Puiseux series for the
eigenvalues at nearby $q$, again by Newton’s method. In that case, we will
need an initial estimate for $a(q)$ correct to $O(q-q^{*})$. This may seem
surprising, if one has only ever used Newton’s method purely numerically, but
the technique is well-known to the symbolic computation community. Either of
these series can be used for continuation: one computes a series about a given
$q$, then uses that series to predict the value of the eigenvalue for a nearby
$q+\Delta q$, which can then be corrected by Newton’s method at the new point.
This may allow larger $\Delta q$, although the danger of branch switching is
always present with too-large a $\Delta q$, and a certain degree of caution is
encouraged.

What allows this series computation to work is that Blanch’s version of the
continued fraction algorithm can be carried out in series. One simply uses
series arithmetic when adding, multiplying, or dividing. This automatically
allows the computation of all derivatives needed. The convergence test only
needs to consider the constant term. More, this allows computation of both
local Taylor series for the eigenvalues, that is

$a(q)=\sum_{k\geq 0}\alpha_{k}(q-q_{0})^{k}\>,$

by carrying out the Newton iteration in series with $q=q_{0}+x$ where $x$ is
the series variable. We are solving

$T(a(x),q_{0}+x)=0$

by iterating

$a^{(k+1)}=a^{(k)}-\frac{T(a^{(k)},q_{0}+x)}{T_{a}(a^{(k)},q_{0}+x)}$

in series; because $x=q-q_{0}$ we get the desired power series. In this case,
we start with the initial estimate $a^{(0)}=\alpha_{0}$, and a single Newton
iteration gets us $\alpha_{0}+\alpha_{1}x$ (plus higher order terms that are
incorrect and we may ignore), and another iteration gets us
$\alpha_{0}+\alpha_{1}x+\alpha_{2}x^{2}+\alpha_{3}x^{3}$ (plus higher order
terms that are incorrect and we may ignore), and so on. The initial estimate
has error $O(x)$; the first iterate has better error $O(x^{2})$; the second
has even better error $O(x^{4})$, and so on, showing a familiar quadratic
convergence; yet somehow nicer than numerical convergence, because more
predictable than numerical Newton’s method in that after $n$ steps we have
error $O(x^{2^{n}})$, and all lower degree terms are (apart from rounding
errors) exactly correct. See [37] for a proof that this method converges “in
series” if the second derivative exists, and for a discussion of the linearly
convergent iteration using the constant derivative $T_{a}(\alpha_{0},q_{0})$
in the denominator instead; that alternative iteration takes more iterations
of course but the series arithmetic is cheaper.

A little thought shows that we may also compute Puiseux series

$a(q)=a^{*}+\sum_{k\geq 1}\alpha_{k}(q-q^{*})^{k/2}\>$

for the eigenvalue about double points, again by carrying out Newton iteration
in series, this time with $q=q^{*}+x^{2}$ where $x$ is the series variable.

For a survey of Puiseux series, see [5]. For a rigorous algorithmic treatment
of expansion of solutions of systems of differential equations in Puiseux
series about singular points, see [19]. In the treatment of Puiseux series in
this paper, which we keep informal so as to maintain readability, we only show
how to compute the first few terms of the series, use them as asymptotic
series only, and do not demonstrate convergence of the resulting series. One
expects, however, by standard results in analysis that the resulting series,
if taken to an infinite number of terms, would in fact converge in a disk
$|q-q^{*}|\leq\rho$ where $\rho$ was strictly less than the distance to the
nearest other double point.

Returning to the problem at hand, we need the initial estimate to be more
accurate than we needed for simple Taylor series: we need the first two terms
correct, namely $a(q)=a^{*}+\alpha_{1}x$ where $\alpha_{1}$ is found by
setting the coefficient of $x^{2}$ to zero in the following series expansion:
$0=T(a(x),q^{*}+x^{2})=$

(50)
$T(a^{*},q^{*})+T_{a}(a^{*},q^{*})(\alpha_{1}x+\cdots)+T_{q}(a^{*},q^{*})x^{2}+\frac{1}{2}T_{a,a}(a^{*},q^{*})(\alpha_{1}x)^{2}+\cdots\>.$

The constant coefficient $T(a^{*},q^{*})$ and the linear coefficient
$T_{a}(a^{*},q^{*})$ are both zero at a double point. The coefficient of
$x^{2}$ is $\alpha_{1}^{2}T_{a,a}(a^{*},q^{*})/2+T_{q}(a^{*},q^{*})$ and so
will be zero if and only if

(51)
$\alpha_{1}=\pm\left(\frac{-2T_{q}(a^{*},q^{*})}{T_{a,a}(a^{*},q^{*})}\right)^{1/2}\>.$

Since the Mathieu equation has only isolated double points, neither
$T_{q}(a^{*},q^{*})$ nor $T_{a,a}(a^{*},q^{*})$ is ever zero181818Certainly
$T_{a,a}$ is never zero because there are only double roots, not triple roots.
If however $T_{q}$ were zero then there would still only be two roots, but in
this case $\alpha_{1}=0$ and $a=a^{*}+\alpha_{2}x^{2}+\cdots$ where
$\alpha_{2}$ is one of two nonzero roots of a quadratic equation. However, we
believe that the theorem of [60] guarantees that $T_{q}$ is never zero so this
should never happen, and indeed we never saw it happen., so $\alpha_{1}$ is
finite and nonzero. These distinct choices for $\alpha_{1}$ lead to distinct
series expansions; together these two series describe the eigenvalues that
merge as $q\to q^{*}$.

With the initial estimate $a^{(0)}=a^{*}+\alpha_{1}x$ we may again use Newton
iteration, even though this time $T_{a}(a(x),q_{0}+x^{2})$ will be $O(x)$
because that derivative is zero when $x=0$. This means that even if $a^{(k)}$
is correct up to $O(x^{m})$, so that the residual $T(a(x),q_{0}+x^{2})$ will
be $O(x^{m})$, we will lose one power of $x$ from the Newton correction and so
$a^{(k+1)}$ will “only” be correct up to $O(x^{2m-1})$. Starting with $m=1$
(i.e. just with $a^{*}$) is therefore not accurate enough; we must have $m=2$
(i.e. start with $a^{*}+\alpha_{1}x+O(x^{2})$) to get off the ground, and then
$2m-1=3$ is higher order, and the next step will have $2m-1=5$, and then $9$,
and so on. This gives a kind of quadratic convergence—still approximately
doubling the number of terms correct with each iteration and after $m$
iterations we will have the series for $a(x)$ correct to $O(x^{2^{m}+1})$—in
computation of the Puiseux series.

See Algorithm 1, which covers both Taylor series and Puiseux series. This
algorithm has been implemented as a Maple procedure and is publically
available at Rob Corless’s GitHub repository. That Newton’s method
automatically converges in formal power series (including Puiseux series) may
be surprising, but it is really the same behaviour as in the numerical world:
the initial estimate has to be close enough, i.e. has to have enough correct
terms in the series, for convergence to start. Once it does, convergence is
rapid. The “asymptotic constant” which complicates analysis in the numerical
world is hidden under the $O(x^{m})$ symbol in the formal power series
analysis, which makes it simpler.

Algorithm 1 Solving $T(a,q)=0$ in series. This algorithm uses Blanch’s
algorithm for evaluating $T(a,q)$, except that operations may be carried out
in series. This allows Taylor series solution near regular points, or Puiseux
series solution at double points.

0: If Taylor series desired, $q_{0}$ and a simple eigenvalue $a_{0}=a(q_{0})$
computed by (say) one-dimensional Newton iteration

0: If Puiseux series desired, a double eigenvalue pair $(a^{*},q^{*})$
computed by two-dimensional Newton iteration, and $T_{a,a}(a^{*},q^{*})$ and
$T_{q}(a^{*},q^{*})$ to compute $\alpha_{1}=\pm 2T_{q}/T_{a,a}$ as in the
text. Choose a sign for $\alpha_{1}$.

0: Positive integer $N$ for the desired number of terms in the series for
$a(x)=a_{0}+a_{1}x+\cdots+a_{N}x^{N}$.

If Taylor series, put $q\leftarrow q_{0}+x$ and $a\leftarrow a_{0}$ and
$n\leftarrow 1$

If Puiseux series, put $q\leftarrow q^{*}+x^{2}$ and $a\leftarrow
a_{0}+\alpha_{1}x$ and $n\leftarrow 2$

while $n<N$ do

$R\leftarrow T(a,q)$ (Trimming leading coefficients $[x^{k}]$ for $k<n$ b/c
rounding errors )

If Taylor series, $n\leftarrow\min(2n,N)$

If Puiseux series, $n\leftarrow\min(2n-1,N)$

$a\leftarrow a-R/T_{a}(q,a)$ to $O(x^{n})$

end while

######  Remark 4.2.

Rounding errors can complicate matters here. In exact arithmetic, the residual
$T(a^{(k)},q(x))$ would be $O(x^{m})$ exactly, for some integer $m$. In
practice, the coefficients of the terms $r_{0}+r_{1}x+\cdots+r_{m-1}x^{m-1}$
are contaminated by rounding errors and while small are typically nonzero.
Especially for the Puiseux series computation, where the derivative starts
with a zero constant term and is $O(x)$, this would mean that the change to
$a^{(k+1)}$ would have spurious nonzero terms of order $1/x$, $1$, $x$,
$\ldots$, $x^{m-1}$. This can rapidly invalidate the results. To make the
algorithm work, then, one must recognize the rounding errors in the
coefficients of the residuals, or simply avoid using terms that one knows
ought to be zero. This is not usually difficult. In our practice, we used
ultra-high precision to check, sometimes working in $100$ or more Digits, that
terms that ought to be zero but looked nonzero were really the result of
rounding errors and not blunders in programming. This allowed us to clearly
distinguish the effects of rounding errors in our experiments.

###  4.4 Examples of Puiseux series about double points

For the double eigenvalue $a^{*}=2.088698902749695\ldots$ corresponding to the
Mulholland-Goldstein double point $q=q^{*}=1.46876861378514\ldots\,i$, we have

(52)
$a=a^{*}+\alpha_{1}\sqrt{q-q^{*}}+\alpha_{2}(q-q^{*})+\alpha_{3}(q-q^{*})^{3/2}+\cdots\>.$

Computation according to the method of the previous section gives that

$\displaystyle\alpha_{1}\approx$ $\displaystyle\qquad\pm
1.65948780432026\ldots(1+\,i)$ $\displaystyle\alpha_{2}\approx$
$\displaystyle\qquad-0.119150377434444\,i$ $\displaystyle\alpha_{3}\approx$
$\displaystyle\qquad\alpha_{1}\cdot(-0.177731786327682\,i)$
$\displaystyle\alpha_{4}\approx$ $\displaystyle\qquad-0.0383269616582290$
$\displaystyle\alpha_{5}\approx$
$\displaystyle\qquad\alpha_{1}\cdot(0.0107135404169547)$
$\displaystyle\alpha_{6}\approx$ $\displaystyle\qquad-0.00154061238466389\,i$
$\displaystyle\alpha_{7}\approx$
$\displaystyle\qquad\alpha_{1}\cdot(0.00273004721440515\,i)$
$\displaystyle\alpha_{8}\approx$ $\displaystyle\qquad 0.000276547402694740$
(53) $\displaystyle\alpha_{9}\approx$
$\displaystyle\qquad\alpha_{1}\cdot(0.000563051707888754)\>.$

Puiseux series can be computed about every double point by the method
suggested in the last section. We do not believe that such series have been
reported in the literature. We see a certain amount of unexplained regularity
in this particular series. Note that this particular $\alpha_{1}$ is a
multiple of $(1+i)$, which is a consequence of the purely imaginary character
of this first double point because $\sqrt{i}=(1+i)/\sqrt{2}$. In table 1 we
tabulate the first few coefficients of some representative series, essentially
as an homage to all the great tabulators of Mathieu functions.

###  4.5 Confirming Blanch & Clemm

Consider figure 9 where the smallest seventy-two double points191919We took
the tabulated values in [14], used Optical Character Recognition to convert
them to computer-readable form, and ran our algorithm (which is the same as
that of [45] and different from that of [14]) from section 4 to confirm them.
For an interesting article connecting OCR with [1], see [76]. Here, after
correcting several amusing OCR errors including the near-inexplicable
occurrence of Russian characters masquerading as numerals, we found (as did
[45]) that all printed decimals in the tables of [14] were correct. Blanch and
Clemm did not plot their double eigenvalues but only tabulated them, but
Hunter and Guerrieri did plot some of theirs. More, they computed farther into
the complex plane than did Blanch and Clemm, and found asymptotic formulae. In
the figure, we see arcs of apparent square-root like curves spreading to
positive infinity; we also see oval arcs of similar beads coming from the
imaginary axis to the real axis, forming a kind of peacock’s tail of double
points., are plotted: as we said, there are a lot of them. Morever, as we saw
above, for values of $q$ close to double points, the vanishing of the norm of
the affected Mathieu functions means that the expansion coefficients must
become singular, and thus numerically troublesome: at the very least, there
will be cancellation error entailed by the subtraction of large nearly equal
quantities. We know of no discussion of this feature of the Mathieu functions
anywhere in the literature. Of course, this is a familiar phenomenon from
other areas of mathematics, such as elementary linear algebra: consider the
analogous problem of finding the eigenvectors of the following matrix, and
expanding another vector as a linear combination thereof:

$\mathbf{A}=\begin{bmatrix}a&1\\\ t&a\end{bmatrix}$

which if $t\neq 0$ has eigenvalues $a\pm\sqrt{t}$, and linearly independent
eigenvectors $\mathbf{v}_{0}=[1,\sqrt{t}]^{T}$ and
$\mathbf{v}_{1}=[1,-\sqrt{t}]^{T}$. Expanding (say)
$[1,1]^{T}=c_{0}\mathbf{v}_{0}+c_{1}\mathbf{v}_{1}$ requires $c_{0}+c_{1}=1$
and $c_{0}\sqrt{t}-c_{1}\sqrt{t}=1$ or $c_{0}-c_{1}=1/\sqrt{t}$. This is
obviously analogous to the situation above. It is even more analogous when one
considers the generalized eigenvector that arises at $t=0$:
$\mathbf{A}\mathbf{u}=a\mathbf{u}+[1,0]^{T}$. Exactly as in the Mathieu
function case above, the numerical difficulties in expanding as a linear
combination of eigenvectors show up for small nonzero $t$, but these are
alleviated on adding the generalized eigenvector to the mix, and writing
instead

(54) $\begin{bmatrix}1\\\ 1\end{bmatrix}=c_{0}\begin{bmatrix}1\\\
\sqrt{t}\end{bmatrix}+c_{1}\begin{bmatrix}1\\\
-\sqrt{t}\end{bmatrix}+c_{3}\begin{bmatrix}0\\\ 1\end{bmatrix}$

Now, of course, the set is not linearly independent, and one has to choose the
coefficients in a sensible way.

(a) Blanch and Clemm double points

(b) double points in all quadrants

Figure 9: In [14] some forty double eigenvalues $a_{k}$ and thirty-two double
eigenvalues $b_{k}$ were tabulated. We plot their results here in figure 9(a).
In the first quadrant, eigenvalues $a_{2k}$ merging with eigenvalues
$a_{2k+2}$ are plotted as blue diamonds. Eigenvalues $a_{2k+1}$ merging with
eigenvalues $a_{2k-1}$ are plotted as black circles. Eigenvalues $b_{2k}$
merging with $b_{2k+2}$ are plotted as red diamonds. Eigenvalues $b_{2k+1}$
merging with $b_{2k-1}$ are plotted as fuchsia circles. We see arcs of
apparent square-root like curves spreading to positive infinity; we also see
oval arcs of similar beads coming from the imaginary axis to the real axis,
forming a kind of peacock’s tail of double points. Extension to all four
quadrants in figure 9(b) follows by conjugate symmetry and by the symmetries
$a_{2k}(-q)=a_{2k}(q)$, $b_{2k}(-q)=b_{2k}(q)$, and
$a_{2k+1}(-q)=b_{2k+1}(q)$, which last implies that the fuchsia circles and
black circles exchange meaning in the left half plane: merging $a_{2k+1}$ and
their conjugates are fuchsia circles in the left half plane, black circles in
the right half plane; merging $b_{2k+1}$ and their conjugates are black
circles in the left half plane, fuchsia in the right half plane. The real axis
is special: only at $q=0$ do $a$ and $b$ eigenvalues merge to $k^{2}$, and
they do so while keeping their independent eigenfunctions, which become $\cos
k\alpha$ and $\sin k\alpha$ at $q=0$.

##  5 Algorithms for Mathieu functions and Modified Mathieu functions

Once one has the eigenvalue, one needs a way to construct the eigenfunction.
Now, our original motivation for studying this problem required accurate
eigenvalues and eigenfunctions as a means to an end: namely the study of blood
flow in a vessel of elliptic cross-section. We consulted the literature and
the available software and were not confident in the selection: there was a
wide variety available, of varying quality and applicability. Most frequently,
the Mathieu functions were computable for real $q$ only, and we needed them
for complex values of $q$. For example, the Mathieu functions in SciPy fall
into this category. We found no software that promised to deal with double
eigenvalues (except possibly Mathematica, but we do not have a license for
Mathematica).

We decided to thoroughly investigate the available software and algorithms and
see what they could do. To that end we decided to work from the ground up: we
would investigate all known methods to compute Mathieu functions. Our goal was
to understand what was going on, not necessarily invent our own method.
However, as a comparison to existing software, we did construct an arbitrary-
precision method that is not too inefficient, which we could assess for
accuracy internally. This gave us a good standard with which to measure other
methods’ accuracy.

We started with the simplest method.

###  5.1 An impractical algorithm

The most straightforward method for computing an entire function is to use its
Taylor series at a convenient point, say the origin. We are guaranteed,
because the series for an entire function always converges, that by taking
enough terms and using enough precision in our arithmetic we can get an
accurate answer. For instance, the basic Mathieu function $w_{I}(a,q,z)$,
called MathieuC in Maple, has the Taylor series beginning

(55)
$w_{I}(a,q,z)=1+\frac{a-2q}{2!}z^{2}+\frac{(a-2q)^{2}-8q}{4!}z^{4}+\cdots$

and since the function is entire, this series converges for all $z$. However,
the series is impractical, (as is well-known to numerical analysts) which we
demonstrate explicitly now.

Taking (say) the modest values $a=29/20$ and $q=3/5$, we still use several cpu
minutes (on a 2017 vintage Windows tablet computer) to compute all the terms
in the above series, using exact arithmetic, up to $O(z^{800})$. This expense
might be considered as pre-computation cost, and ignored, and so we do. This
truncated series can then be used to compute, for instance

(56)
$w_{I}(a,q,4.0i)=-0.13260687218535639758416285435423{\color[rgb]{1,0,0}86957557534433}\>,$

using 46 Digits of precision in the computation; the final nonzero term of the
truncated series is about $-1.4\cdot 10^{-47}$ and because all subsequent
terms are smaller202020The series is not an alternating series because the
sign pattern is more complicated. This is not a rigorous error estimate,
therefore., we would hope that the computed value is accurate. But the final
12 Digits, coloured red, are incorrect, because some intermediate terms in the
series are about $10^{12}$ in size and rounding errors in those terms are
revealed by the cancellation of the large terms. This is a well-known effect,
of course, discussed in many textbooks and educational papers; see e.g. [25],
which discusses it as an effect of the ill-conditioning of the truncated
Taylor series.

The situation gets very much worse as $|z|$ increases. Already by $z=4.36i$
the best accuracy we can achieve with the $800$ term series is double
precision, because although in exact arithmetic the truncation error is about
$10^{-17}$, the condition number is about $10^{20}$ and so we have to carry
$20$ extra digits; by $z=4.45i$ we can only achieve single precision, and that
by using $22$ extra digits; and by $4.51i$ only half precision by using $24$
extra digits. By $z=5.0i$ it is already true that $800$ terms are not enough
in the series, no matter what precision we use because the truncation error is
too big. Moreover, even if we had enough terms, because the condition number
is about $10^{40}$, the computation would need about $40$ extra decimal digits
of precision anyway.

The final conclusion is that this method is unaffordable. It can be rescued by
using analytic continuation instead of a single series (which of course is the
underlying basis for most numerical methods to integrate ODE), and we look at
that method in section 5.2. Most researchers instead choose to jump straight
to a spectral method for the periodic Mathieu functions, expanding in Fourier
series; and then by a trick these can be used also for the modified Mathieu
functions. We will discuss this method in section 5.3.

###  5.2 Analytic continuation (marching)

One straightforward thing to try, valid for both the Mathieu functions and for
the modified Mathieu functions, is numerical solution of the Mathieu equation.
We could in principle use any method that allows integration along a complex
path, such as any Runge-Kutta method or multistep method. Because these
solutions are to be later used as functions, requiring evaluation at any point
in their domain, we will need good interpolants and modern codes automatically
supply these.

However, because the equation is linear, and because its variable coefficient
has a known Taylor series algorithm, it is a straightforward exercise to
develop a specialized numerical solution by marching Taylor series, a method
known as analytic continuation if infinite Taylor series are used and as a
Taylor series method if a finite order is used at each step. Taylor series
methods are often introduced in numerical analysis textbooks and then
dismissed practically in the same breath as being too costly or insufficiently
general, but these methods do not actually suffer from those faults when
properly implemented: see for example [67] which shows how to use them to
solve DAE but also gives many references to papers using them to solve
initial-value problems (IVP) for ODE. See also [38] for a discussion of their
use for computation of special functions. Taylor series methods are especially
attractive when the differential equations are very simple, but they work very
well on most smooth problems. For such problems, Taylor series methods are of
cost polynomial in the number of bits of accuracy required, on finite
intervals, as the number of bits of accuracy required goes to infinity [46];
in contrast, fixed-order methods such as any given Runge–Kutta method are of
cost exponential in the number of bits of accuracy required. So Taylor series
methods are actually a quite natural method to try when high accuracy is
wanted, which might surprise the casual user of mathematical software for
solving ordinary differential equations.

One practical advantage of such methods is that they come with free
interpolants, and another advantage is they come with an inexpensive estimate
for the residual. A third advantage is that the residual can be accurately
measured afterwards (as opposed to estimated beforehand) to validate the
solution. These methods can be implemented as variable stepsize and as
variable order. In the case of stiff problems they can be implemented as
implicit methods [6]. Indeed, these methods are interesting for other
problems, not just the Mathieu equation.

So, that’s what we did. In fact, we have implemented (in Maple) what is called
a Hermite-Obreschkoff method, which is implicit and uses Taylor series at both
ends of the marching step, and is more numerically stable than explicit Taylor
series methods. Implicit methods are more effective for stiff problems, as is
well-known [78].

The Mathieu equation is not usually stiff, per se, but rather is frequently
oscillatory (especially along the imaginary $z$-axis, i.e. for modified
Mathieu functions). This causes a pure Taylor series method to suffer some
instability. In contrast, the Hermite-Obreschkoff method, which is implicit
(something like a generalization of the implicit midpoint method), performs
more satisfactorily.

We implemented this, not because we thought it would be the best method to
evaluate Mathieu functions, but because we thought the implementation would be
useful in our quest for understanding of Mathieu functions, and moreover
provide an independent check on other methods that we explore. The fact that
such methods are possibly interesting for other problems is merely a bonus.

####  5.2.1 Specifics of the Hermite-Obreschkoff method we use

This description is quite short because the basic technology is widely
understood and our implementation is not the main focus of this paper. We
provide details only for reproducibility. The key piece is the recurrence
relations for the Taylor coefficients of the expansion of $y(t_{n}+\Delta
t)=\sum_{k\geq 0}w_{k}\Delta t^{k}$. Simultaneously, we need the Taylor
coefficients for $\cos 2t$, which necessitates the Taylor coefficients for
$\sin 2t$. This gives

$\displaystyle C_{k+1}=$ $\displaystyle-\frac{2}{k+1}S_{k}$ $\displaystyle
S_{k+1}=$ $\displaystyle\frac{2}{k+1}C_{k}$ (57) $\displaystyle w_{k+2}=$
$\displaystyle-\left(aw_{k}-\frac{2q}{(k+1)(k+2)}\sum_{j=0}^{k}C_{j}w_{k-j}\right)\>.$

$C_{0}=\cos(2t_{n})$ and $S_{0}=\sin(2t_{n})$, while $w_{0}=y(t_{n})$ and
$w_{1}=y^{\prime}(t_{n})$. This recurrence needs to be scaled if we are
integrating on a complex path.

We interpolate over the interval $t_{n}\leq t\leq t_{n+1}$ by using what we
call a “blend”: this is nothing more than Hermite interpolation using the
Taylor coefficients at each end of the step [24]. A string of such blends put
together as a piecewise polynomial is called a “string of blends” because it
is tied together at “knots”. Such interpolants are quite remarkably stable
numerically, even at very high order, and can be evaluated in cost linear in
the degree of the interpolant. The high order means that they are potentially