Computation and applications of Mathieu functions: A historical perspective

Mathieu functions of period $\pi$ or $2\pi$, also called elliptic cylinder functions, were introduced in 1868 by \'Emile 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: \'Emile Mathieu, Sir Edmund Whittaker, Edward Ince, and Gertrude Blanch.


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 π and 2π 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 L A T E X-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.

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 ωt, which consists of a steady part plus a periodically-varying part, can one choose H and ω 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 where they take F (t) = P +H cos ω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. "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.

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 Fig. 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].
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 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.
2.3. Vibrating membrane bounded by an ellipse. 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- 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.
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: 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 (3.1) by the transformation z = ±ix (the sign makes no difference): The even solutions are conventionally written as Ce g (z, q) = ce g (±iz, q). The eigenvalues for even solutions are conventionally written a g (q). The odd solutions are written similarly, as se g (z, q) and 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 α and β as using the fundamental pair of solutions satisfying w I (0; a, q) = 1 with w I (0; a, q) = 0 and w II (0; a, q) = 0 with w II (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 y(z) = αw I (±iz; a, q) + βw II (±iz; a, q) . (3.4) 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) 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. Mathieu (1835Mathieu ( -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, circular 1 . Mathieu's discussion of the circular case starts with the PDE ∂ 2 w ∂t 2 = m 2 ∂ 2 w ∂x 2 + ∂ 2 w ∂y 2 (Mathieu used d and not our modern Russian ∂ for partial derivatives) and then transformed x = r cos α, y = r sin α to polar coordinates to get

3.1.Émile Léonard
Thereafter Mathieu used what is now the standard method of separation of variables for a pure oscillation w = sin(2λmt)u(r, α) and u(r, α) = P (α)Q(r) to give a harmonic equation for P (α) and is equivalent to what we now call Bessel's equation for Q(r): [Bessel's equation in standard form is z 2 d 2 y/dz 2 + zdy/dz + (z 2 − ν 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 transformation 2 x = c cosh(β) cos(α) and y = c sinh(β) sin(α) where 2c is the distance between the foci of the ellipse 3 , 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: "Comme le premier membre ne peut renfermer que α, et le second que β, ils sont egaux 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 α and cosh 2 β. Using double-angle identities, Mathieu later in this same paper transformed these to something that we write in modern notation as Notice that these two equations can be transformed into each other by the change of variable β = iα. Here q = λ 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 α = (1 + cos 2α)/2 and cosh 2 β = (1 + cosh 2β)/2 2 Mathieu initially spelled the hyperbolic functions out explicitly, and later used "Pour simplifier" the notation E(β) = (e β + e −β )/2 for cosh(β) and E(β) = (e β − e −β )/2 for sinh(β). 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. 3 Notice 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.
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 α runs from 0 to 2π (or, of course, from −π to π), requiring periodicity. Often the problem is taken to have π-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 Sturm 4 (1803-1855) implied that, for real q, what we now call the Sturm-Liouville form of the Mathieu equation could be written as (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) < · · · , 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 π they are denoted by ce 2k (x), for k ≥ 0. If they are even and of period 2π they are denoted by ce 2k+1 (x), for k ≥ 0. If they are odd and of period π they are denoted by se 2k (x) for k ≥ 1. If they are odd and of period 2π they are denoted by se 2k+1 (x) for k ≥ 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 We prove this in detail in appendix B. By the proposition C we mean here a proposition that is true if y k and y are in the same class; that is, both are even and of period π, both are odd and of period π, both are even and of period 2π, or both are odd and of period 2π. We have used here the Iverson convention [·] 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 π instead of 2π. The exact method used for normalization is a matter of convention. Finally, these classes of orthogonal eigenfunctions are complete: every reasonable 5 even/odd period-p function f (x) can be expanded in a convergent series Periodic functions which are neither even nor odd can of course be written as a sum of an even function and an odd function and thus will use both classes of the given period in its expansion. Since the coefficients α 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. 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 y k , y k = 0 and thus normalization by making the bilinear form equal to a nonzero constant is not possible. We will return to this later.

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 (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α in its series expansion to be unity (and similarly the coefficient of sin gα 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α 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) = ace g (x) + bse g (x) and (a) Some ce 2k (α) for q = 1.5 (b) Some se 2k (α) for q = 1.5 (c) Some ce 2k+1 (α) for q = 1.5 (d) Some se 2k+1 (α) for q = 1.5 we specify that ce g (x) satisfy not only ce g (0) = 0 but also ce g (0) = 1, and similarly se g (0) = 0 and se g (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 (3.10) is nominally enforced to be π (or 2π), 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 general 6 . 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π solution of equation (3.7) 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 normalization 7 and to fewer terms than he calculated to is: 3 − g 2 + 4 g + 7 cos (g + 2) α 128(g + 2) (g + 1) 3 (g − 1) − cos (g + 6) α 384(g + 2) (g + 3) (g + 1) 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 (3.13) a = g 2 + h 4 2(g − 1) (g + 1) 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. ( and alternatively by µ = sin α, (Mathieu used ν rather than µ, which we originally kept; but a referee pointed out that this was confusing, and on second thought we agreed) whereupon the Mathieu equation becomes The DLMF gives yet another algebraic form, using the change of variables ζ = sin 2 α. 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.3) in section 5, while the recurrence relations for the algebraic formulations are already given in [58]. For instance, for equation (3.14) if then (after the first few terms which have to be separately investigated), We computed this recurrence relation automatically from equation (3.14) by using the gfun package in Maple, specifically its diffeqtorec command [74]. Mathieu gave the simpler form at ν 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 α or of sin α. In particular, he used them to count real roots.

Modified Mathieu Functions.
In order to solve the vibrating drum problem, Mathieu had also to solve equation (3.8). 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: Ce g (x; q) = ce g (ix; q) and Se g (x; q) = −ise 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 Ce 3 (β; q)ce 3 (α; 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 Ce 3 (β; q). Since the coordinates are x = c cosh β cos α, y = c sinh β sin α (alternatively, x + iy = cos(α − iβ)) we will want c cosh β = 5 and c sinh β = 3; this gives β = ln 2 and c = 4. Now we want the value of q so that Ce 3 (ln 2; q) = 0. By zerofinding on q, we find that q ≈ 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λmt via q = λ 2 c 2 , with the membrane parameter m. This value of q gives the eigenvalue a 3 ≈ 14.6695. The hyperbolic nodal lines are at approximately α = ±0.9857, ±π/2, and ±2.156. The contours plotted were at levels ±[0, 11,22,33,44]/40, remembering that our normalization is so that ce 3 (0; q) = 1 = Ce 3 (0; q). More details and more figures can be found in [22].
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 wellregarded 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 con- As the argument β 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 ±c.) (Right) A graph of Se 1 (β; q) = (se 1 (iβ; q)) when q = 2.
tends, 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.  [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: β is more like a radius, and α 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 ce g (z) to be unity, and similarly the coefficient of sin gz in the expansion of 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 (3.10) 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 ODEs 8 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(µz)φ(z) where φ(z) is periodic with period π (because the periodic forcing of the Mathieu equation has that period) and µ (actually, in modern works starting in [1], ν where µ = iπν) is called the characteristic exponent. Regions where (µ) > 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(−µz)φ(−z) is also a solution; this implies that for both solutions to be stable, we must have (µ) = 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(η) is an even Mathieu function, then there is a characteristic number λ for which (translating from the 16q convention of Whittaker to the 2q convention used in this paper) 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 ζ = 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 ν = cos z and µ = 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.

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

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 α or sin α 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: Inserting this series into the Mathieu equation (3.1), using the multiplication identities 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: and, for all k ≥ 2, The odd cosine coefficients must instead satisfy and, for all k ≥ 1, The even sine coefficients give and, for all k ≥ 2, Finally, the odd sine coefficients give and, for all k ≥ 1, 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 √ 2, we may symmetrize even the first one.
The infinite eigenvalue problem then becomes The eigenvalues of this matrix are denoted a 0 (q), a 2 (q), a 4 (q), . . . and indeed for real q these occur in increasing order: a 0 (q) < a 2 (q) < a 4 (q) < · · · . 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: The eigenvalues of equation (3.28) are denoted a 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] 9 . 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 fraction 10 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.
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].  . 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 ≤ k ≤ 30 by using a matrix of dimension k + ∆n where ∆n is as indicated on the figure: either 1 ≤ ∆n ≤ 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, . . ., 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.
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.

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 = √ q and s = 2m + 1, then as h → ∞ with m = 0, 1, 2, . . ., 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  where Re(µ) > 0 and φ(z) is periodic; we can see that this particular periodic function has two zeros by the "spikes" in the graph. algebra systems, for instance Maple.

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 π or 2π". This is generally taken to mean that the other linearly-independent solution cannot be periodic, and indeed this is true for periods π and 2π. 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 ce p (z), then the second, necessarily non-periodic, solution would be where both K and φ(z) would be O(q p ) in size and φ(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 ≈ 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π, 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]. 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 π or 2π 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.

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 eigenvalues 12 .
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 for m = 0, 1, 2, . . .. 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 (3.31) is equation (1.05) in [13].
We use this notation in the recurrence relations (3.18)- (3.26), 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 Putting G m = A m /A m−2 and dividing equation (3.32) 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 This recurrence can be written either as Of course, these recurrences must be started correctly by using the appropriate equations (3.18)-(3.26). Running equation (3.35) until the denominator Q is so large that 1/|Q 2 | ≤ ε where ε is our tolerance, starting from some M > 2 so that c k = 1 for all k ≥ M − 2,  T (a, q) must be zero for a to be an eigenvalue corresponding to q. The edge cases for G M,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 (3.37), but since differentiation of T (a, q) with respect to a is simply a matter of differentiating equations (3.34) and (3.35) (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 eigenvalue 13 , the iteration is 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].

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) ) + 1 2 T aa (a (k) , q)(a − a (k) ) 2 + · · · 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. Hurwitz 14 . 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].

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 µ = iπν, 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) = e µz φ(z) + e −µz φ(−z) 2φ(0) 14 Blanch 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. 15 Philip Morse, at MIT, the first author of [64]. This monumental work also describes Mathieu functions. or w II (z) = e µz φ(z) − e −µz φ(−z) 2(φ (0) + µφ(0)) will have nonzero denominator; now differentiate w II and compute w I (π) and w II (π).
In the DLMF this equation (using cosine and not hyperbolic cosine, because ν is used instead, where µ = iν) is called the characteristic equation, number 28.2.16. Blanch uses another convention, namely µ = iπν, defining ν differently, in [1]. The second interesting fact is that the periodic solutions, that is the Mathieu functions, correspond to the case µ = i n where n is an integer; if n is an even integer the solution is periodic with period π and if n is odd the solution is periodic with period 2π. The case when µ = 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. . .".
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 * ≈ 1.468768613785142 i (reporting 16 digits 16 ), we have a double eigenvalue: a 0 = a 2 ≈ 2.088698902749695. Several interesting things happen at double points. First, the eigenfunctions coalesce: here, ce 0 (q * , α) = ce 2 (q * , α), 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 (q) < (q * ) we have a 0 (q) < a 2 (q) < a 4 (q) < · · · , but at q * equality occurs. For (q) > (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 (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.
By numerical experimentation at high precision 17 we find that when we use the normalization convention that ce 2k (q, 0) = 1, and moreover that where d ≈ 1.659487804320256 + 1.659487804320256 i = 1.659487804320256(1 + i). These values were found by using orthogonality, which holds if q = q * , and by high- 16 This double point was the first found: studied in [65] and later computed by [18] to 3 digits and then to double precision in [14]. 17 We 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.
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 ∂y/∂a (found by solving a variational equation). By combining the equations above, we can see that

This arrangement is continuous as
in that limit. This shows explicitly that by adding ∂ce 0,2 (q * , α)/∂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(α) in section 5.4.

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: where all function evaluations and derivative evaluations occur at the current estimates (a (k) , q (k) ). Then as usual a (k+1) = a (k) + ∆a and q (k+1) = q (k) + ∆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 √ 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+∆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 This must be equal to mπ/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.

Computing local series for the eigenvalues by Newton's method.
We may use Newton's method not just to compute a simple eigenvalue, or twodimensional 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 + ∆q, which can then be corrected by Newton's method at the new point. This may allow larger ∆q, although the danger of branch switching is always present with too-large a ∆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 by carrying out the Newton iteration in series with q = q 0 + x where x is the series variable. We are solving in series; because x = q − q 0 we get the desired power series. In this case, we start with the initial estimate a (0) = α 0 , and a single Newton iteration gets us α 0 + α 1 x (plus higher order terms that are incorrect and we may ignore), and another iteration gets us α 0 +α 1 x+α 2 x 2 +α 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 (α 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 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 * | ≤ ρ where ρ 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 * + α 1 x where α 1 is found by setting the coefficient of x 2 to zero in the following series expansion: 0 = T (a(x), q * + x 2 ) = (4.11) T (a * , q * ) + T a (a * , q * )(α 1 x + · · · ) + T q (a * , q * )x 2 + 1 2 T a,a (a * , q * )(α 1 x) 2 + · · · .
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 α 2 1 T a,a (a * , q * )/2 + T q (a * , q * ) and so will be zero if and only if (4.12) Since the Mathieu equation has only isolated double points, neither T q (a * , q * ) nor T a,a (a * , q * ) is ever zero 18 , so α 1 is finite and nonzero. These distinct choices for α 1 lead to distinct series expansions; together these two series describe the eigenvalues that merge as q → q * .
With the initial estimate a (0) = a * + α 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 * + α 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 4.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.
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+· · ·+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, . . ., 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.
Algorithm 4.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. Require: If Taylor series desired, q 0 and a simple eigenvalue a 0 = a(q 0 ) computed by (say) one-dimensional Newton iteration Require: 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 α 1 = ±2T q /T a,a as in the text. Choose a sign for α 1 . Require: Positive integer N for the desired number of terms in the series for a(x) = a 0 + a 1 x + · · · + a N x N . If Taylor series, put q ← q 0 + x and a ← a 0 and n ← 1 If Puiseux series, put q ← q * + x 2 and a ← a 0 + α 1 x and n ← 2 while n < N do R ← T (a, q) (Trimming leading coefficients [x k ] for k < n b/c rounding errors ) If Taylor series, n ← min(2n, N ) 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.
Computation according to the method of the previous section gives that 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 α 1 is a multiple of (1 + i), which is a consequence of the purely imaginary character of this first double point because √ i = (1 + i)/ √ 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. figure 9 where the smallest seventy-two double points 19 , 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:

Confirming Blanch & Clemm. Consider
which if t = 0 has eigenvalues a ± √ t, and linearly independent This is obviously analogous to the situation above. It is even more analogous when one considers the generalized eigenvector that arises at t = 0: Au = au + [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 Now, of course, the set is not linearly independent, and one has to choose the coefficients in a sensible way.

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 19 We 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. Fig. 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 squareroot 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α and sin kα at q = 0. 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.

(a) Blanch and Clemm double points (b) double points in all quadrants
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 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 (5.2) w I (a, q, 4.0i) = −0.1326068721853563975841628543542386957557534433 , using 46 Digits of precision in the computation; the final nonzero term of the truncated series is about −1.4 · 10 −47 and because all subsequent terms are smaller 20 , 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.

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 + ∆t) = k≥0 w k ∆t k . Simultaneously, we need the Taylor coefficients for cos 2t, which necessitates the Taylor coefficients for sin 2t. This gives C 0 = cos(2t n ) and S 0 = sin(2t n ), while w 0 = y(t n ) and w 1 = y (t n ). This recurrence needs to be scaled if we are integrating on a complex path.
We interpolate over the interval t n ≤ t ≤ 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 competitive with spectral methods in efficiency for smooth problems, when high accuracy is desired [46]. They are simple to integrate and differentiate, and reasonably simple to multiply together, and this gives inexpensive ways to evaluate integrals containing these numerical solutions, e.g. π 0 f (z)ce 3 (q, z) dz. There are also simple methods based on companion matrices to find zeros of blends; alternatively, there is an iterative scheme of higher order than Newton's method (1 + √ 3 ≈ 2.732 instead of 2) for the same cost per iteration [21]. One can imagine implementing an analogue of Chebfun [7] using blends, and this might be interesting because the approximation properties are different (mostly not as good, but in some niche circumstances when many derivatives are easily computed at isolated points in the complex plane, possibly competitive). An analogue of Chebfun using sinc functions [72] has some interesting properties; perhaps an analogue of Chebfun using blends would be worth investigating. A final advantage for this method is that it works directly on the Mathieu equation and the modified Mathieu equation both, so that both Mathieu functions and modified Mathieu functions can be computed with the same code merely by changing the initial conditions. The Hermite interpolant ("blend") essentially doubles the order of the method for the same number of Taylor coefficients taken at each step, and moreover greatly improves stability: when the degrees of the Taylor coefficients are the same at each end (as we use) this is quite stable [67]. Indeed on the Dahlquist test problem, a balanced implicit Taylor series method is A-stable (this was known as early as 1991 [52]). Because the problem is linear, the method winds up being only linearly implicit.
To take a step, we predict the stepsize by the PI step control suggested in [42], which uses data from the previous two steps (N.B.: at the beginning we use a predictor based on Taylor series alone). Then we generate two separate Taylor series at t n + ∆t, one (say y c (t)) with initial conditions y(t n + ∆t) = 1 and y (t n + ∆t) = 0 and another (say y s (t)) with y(t n + ∆t) = 0 and y (t n + ∆t) = 1; we then use collocation at the points τ 1 = t n + ∆t/4 and τ 2 = t n + 3∆t/4 (these are Chebyshev-Lobatto points): we blend the Taylor series at t n (which we knew already) with a linear combination z = αy c (t) + βy s (t) of the new ones at t n + ∆t with the linear combination chosen to make the residual r(t) = z − (a − 2q cos(t))z equal to zero at τ 1 and τ 2 . This gives us two linear equations in the two unknowns α and β. Solution of the 2-by-2 system is hand-coded in. The second derivative z (t) of the blend z(t) is computed by program differentiation (actually hand-coded, so not technically "automatic" differentiation: call it "semi-automatic differentiation").
We then measure the residual of the resulting blend at s = t n + ∆t/2. This is asymptotically (as ∆t → 0) the location of the maximal residual of this blend. If this residual is small enough, smaller than the specified tolerance, we accept the step and use this estimate in the PI control to predict the size of the next step.
This therefore gives us an independent, reliable, and reasonably efficient specialpurpose numerical solver for the Mathieu equation. In our experience it performs well, especially because we can measure the residual r(t) separately at many points, to provide reassurance that we have obtained the exact solution to a differential equation very near to the Mathieu equation. Trivially, we have solved y + (a − 2q cos 2t)y = r(t) denoting the residual by r(t), and we have chosen stepsizes to ensure that |r(t)| is small. By the Alexeev-Gröbner nonlinear variation of constants formula this means that the global error will also be small: there exists a function G so that (with the same initial conditions for the reference solution y(t) and the computed solution z(t) solving the above) Indeed because the Mathieu equation is linear we may write G(t, τ ) explicitly as a Green's function: We will use this again later. The Wronskian of the Mathieu equation is 1. We did not plot G(t, τ ) but instead merely estimated its effects by integration at different tolerances. Outside of regions of exponential growth, i.e. where the real part of the characteristic exponent is zero, the Mathieu equation is well-conditioned. Even when (µ) > 0 this method allows us to keep the global error relatively small over short enough integration paths.
In short, we wrote a special-purpose numerical ODE solver to evaluate the Mathieu functions. By varying the initial conditions, we may compute either of the basic solutions, w I (t) or w II (t); once one has an eigenvalue given q this already allows computation of any of the Mathieu functions. We may compute the other independent (nonperiodic) solution as we did for figure 8(a). By integrating on a complex path, we may compute any of the modified Mathieu functions. The modified Mathieu functions, which are also needed already for the solution of the vibrating elliptic membrane, are of course not periodic. We can compute Floquet solutions as we did for figure 8(b), and estimate the Floquet exponent numerically [33] and compare it with Poincaré's formula for µ (equation (3.38)), see figure 10. It is also possible to use this method with the Green's function in equation (5.4) to compute the generalized eigenfunctions needed for double eigenvalues (see the discussion in section 5.4).
In spite of these capabilities, we are not recommending our implementation of this method as a general-purpose Mathieu solver, both because we have not made it bulletproof and because we believe that for most purposes it would not be as fast as a spectral solver based on Fourier series-after all, with a Fourier method all the Fourier coefficients are essentially immediately available once the eigenvectors have been computed-but for our purposes, which included an independent check on existing software, it was adequate. By varying the precision at which we computed and by taking tolerances to be extremely small we were able to buy accurate solutions by paying for computing cycles. That said, it is a variable order method that can tolerate quite high order 21 without numerical difficulty, and so is quite fast. We have not made a detailed speed comparison with a spectral method. However, for the Mathieu functions on which we concentrate this paper-which are periodic-one expects that the spectral method of expanding in Fourier series to be more efficient than any possible marching method. For the modified Mathieu functions, the issue is not so clear-cut. Again, though, we wanted an independent and reliable method  Here we track (µ) where the Floquet solutions are exp(µz)φ(z) and exp(−µz)φ(−z) with φ(z) periodic with period π. Regions near the a axis are stable as q is increased in a purely real fashion, except only neutrally so at q = 0 and a = g 2 for integers g; regions near the a axis are stable as s is increased where q = is, except as before only neutrally so at a = g 2 . Near those difficult-to-contour regions we supplemented the graph with Taylor series expansions of the characteristic curves, in red. Double eigenvalues are indicated with blue dots. as a comparator for other methods, and in this the Hermite-Obreschkoff code we implemented is perfectly satisfactory.

Complexity and other issues.
One interesting complication is that when using the Hermite-Obreschkoff or Taylor method (or, indeed, when using any method), accurate solution becomes very expensive for large |z|. This is because solutions to the modified Mathieu equation become highly oscillatory along the real axis (so the solutions to the Mathieu equation become highly oscillatory along the imaginary axis). Indeed the asymptotics of the solution to the modified Mathieu equation contain a term exp(i 2 √ q cosh(z)) (see Formula 28.25.1 in the DLMF; we ignore the denominator, which grows only more slowly). This is confirmed by inspection of the formulae in Chapter 28.23 of the DLMF (https://dlmf.nist.gov/28.23) where we see the arguments 2 √ q cosh(z) and 2 √ q sinh(z) appearing inside the Bessel functions used in the expansions for the modified Mathieu functions. Resolving any solution for graphical purposes requires computing a fixed number of points in each cycle, say 5 or 10. One therefore sees that the number of computed points required to resolve the solution grows exponentially with the real part of z, because the number of cycles grows exponentially with the width of the interval. We therefore see that the cost to directly compute-in this sense-the solution accurately must grow exponentially with the argument. This in part was why the simple Taylor series at the origin failed already by z = 5.0i, in section 5.1.
This does not contradict the theorem in [46] because that only holds on a fixed finite interval, in the limit as the required accuracy goes to infinity.
Off the real or imaginary axes, if z = x+iy has nonzero x and y, then the solution not only oscillates but grows doubly exponentially as x grows, containing the terms exp(±2 √ q sinh(x) sin(y)). This fact does not seem to have been explicitly remarked on elsewhere, perhaps because the information is readily visible in the formula, and the user is expected to understand this without being told. That fact suggests that most applications of the Mathieu functions must have only small |z| or else either purely real z or purely imaginary z; that is, either the real Mathieu functions or the real modified Mathieu functions will be the relevant functions. Asymptotics cannot recover the phase information, of course-so if you want an accurate solution for large z, you will have to pay a significant price for it (and know z to exponential accuracy, in a certain sense). This conclusion may be surprising, so take the simpler function S(x) = cos(cosh x) (which occurs when setting q = 1/4 in formula 28.25.1 of the DLMF and taking only one part of the complex exponential and ignoring all slower-varying terms). If x is about 1000 or so, how many bits of x are necessary to know before you can extract one bit of information, namely the sign, of S(x)? Since cos θ = 0 when θ = (2k + 1)π/2 for some integer k it follows that we must know x so well that we can detect that it is between x k with cosh x k = (2k + 1)π/2 and x k+1 with cosh x k+1 = (2k + 3)π/2. Now x k+1 − x k = inv cosh((2k + 3)π/2) − inv cosh((2k + 1)π/2) ∼ 1/k − 1/k 2 + . . . as k → ∞. But how large is k? If x k is about 1000, then (2k + 1)π/2 is about cosh 1000, or k is about 2 1440 . So 1/k is about 2 −1440 . To detect the difference between x k+1 and x k , then, we need to carry twice that, or 2880 bits. That corresponds to about 180 hexadecimal digits, or about 240 decimal digits. For x = 10 4 instead, we need to know about 2400 decimal digits of x; all this just to get the sign of S(x) correct. Yet we know the asymptotics of this function perfectly well-it's just a cosine with exponentially increasing frequency. We remark that the cost of arithmetic with big floats grows at least linearly with the length of the float; FFT multiplication is a logarithmic factor more costly (and naive multiplication grows quadratically in cost).
This discussion suggests that any numerical method that we use will only be helpful and efficient for |z| bounded by a modest constant, say 2π. For larger values of z one must use the asymptotic formula instead, and give up on phase information.
Remark 5.1. The highly oscillatory nature of the solutions to the modified Mathieu equation also induce an instability for large |z| in the Taylor series method, akin to the instability induced for stiff problems [23,78]. If the stepsize is not aggressively reduced as |z| increases, then eventually this instability takes over and the numerical solution becomes meaningless and rapidly overflows. One would hope that instead using an implicit Taylor series method would alleviate this problem, and to a certain extent it does, but even with an implicit method the stepsize is forced to be so small for accuracy that failure is assured in practice when |z| is at all large.
Remark 5.2. The Bessel function expansions in https://dlmf.nist.gov/28.23, and similarly those of [32] which we will discuss presently, seem to offer a way around the difficulty, at least when a is an eigenvalue of the problem and the solution is periodic. Accurate computation of the series coefficients, together with accurate computation of the appropriate Bessel functions, seem to offer an inexpensive way to accurately evaluate a modified Mathieu function. Inspection of the derivative, however, shows that the large argument z = x + iy also needs to be known to exponential accuracy in order to find an accurate value of the Mathieu function in question, especially near a zero. The difficulty seems to be intrinsic. Since hardly anyone complains about this in the literature, we conclude that most people only want the values of Mathieu functions and modified Mathieu functions for small or at most moderately large values of z, or magnitude information and not phase information when z is large.
Mathieu himself thought there might be issues for large z and introduced the change of variable ν = cos(z), which leads to the following algebraic differential equation (ADE) (equation 28.2.3 in the DLMF). Recall that an ADE is generally different from a DAE.
This algebraic differential equation (and indeed also the similar one in equation 28.2.2 in the DLMF that arises on ζ = sin 2 z) has some interesting computational properties: for one, they are what is called D-finite or holonomic, meaning the recurrence relation for the Taylor series terms is of fixed order, and which means accurate computation can be done asymptotically quickly [86]. See also [61,62,8].
Indeed near the imaginary axis this DE is also numerically easier-one may take more nearly equal integration steps in the new variable, instead of the ever-decreasing ones necessary in the original variable z-but somehow this is just "sweeping the problem under the rug" because the value of ν itself becomes very large; if the cost of the stepsize selection is low (and the control is effective) one should get very nearly the same numerical performance in the original variable, except that the recurrence relation for the Taylor series terms in the original variable depends on all previous terms, not just a fixed number as for D-finite functions and so each step is also more expensive. However, the presence of the singularities ν = ±1 complicates mattersindeed the initial condition in α at α = 0 corresponds to a singular point ν = 1 and so something special must be done to get the integration started, and moreover the branching structure of the nonlinear transformation also makes its presence felt because the inverse transformation is multivalued. Assessment of the resulting solution by the method of residuals also needs an extra step. So although all these obstructions can be dealt with, it is not actually clear which of these two approaches is numerically best. Experiments seem necessary. Since our purposes were served by the Hermite-Obreschkoff method in the original variable, we did not pursue this.

Spectral methods.
However, the method of choice, at least for most implementers, for computation of Mathieu functions is the use of what is effectively a spectral method. This is prima facie valid only for the case when a is an eigenvalue and the solution is periodic. One computes the eigenvalues by the matrix method as in section 3.3.1, and then the resulting eigenvector gives the coefficients in the Fourier series expansion of the Mathieu function (if using the continued fraction instead, the recurrence relations can be used, although one has to take numerical care). The Fourier coefficients decay extremely rapidly, as is usual for Fourier expansion of smooth functions. To give an example to show just how rapidly, here is Equation 28.4.24 from the DLMF (here the superscripts refer to which eigenfunction and the subscripts refer to which coefficient in its Fourier expansion): .
This holds as m → ∞, for fixed n. This states that the mth Fourier coefficient ultimately decays like (q/4) m /(m!) 2 , much faster than exponentially. The other three types of Fourier coefficients decay similarly rapidly. This very rapid decay means that good approximations can be made with only a few terms of the Fourier series.
For the corresponding modified Mathieu functions these Fourier series also converge, but now only slowly. As an alternative, one can use the same Fourier coefficients in a Bessel function expansion, for instance equation 28.23.2 of the DLMF: where the modified Mathieu function on the left can be approximated by the series on the right; the notation therein is different from the notation in our paper, but the Cs are Bessel functions and the c 2n s are the Fourier coefficients for the ordinary even period-π Mathieu functions. These are not the only series one might use. In [32] (and in [1] and the DLMF) we find the following expansions, which the authors claim are rapidly convergent: Here s = √ q exp(x) and t = √ q exp(−x). We find these series to be preposterous: the A k and the B k are the Fourier coefficients of the corresponding Mathieu functions, namely the same coefficients as in the Fourier series for the corresponding Mathieu function. Now they are to be used in a completely different, almost alien-looking, series? But these preposterous formulae are both correct and useful, and go back at least to [26]. Erdélyi thought the 'coincidence' of these series coefficients to be significant [29]. Here we have translated to the notation of this paper-although the normalization used in the above does not agree with the normalization here, even though the authors of [32] claim to use the same normalization that we do; instead the formulas are the same as those in [30]. We tried this, and it worked well. Indeed, such a Fourier-Bessel method is similar to those advocated by almost everyone, from [3] to [93]. Against this Fourier-Bessel method it is not clear that a straightforward numerical solution such as the one we have implemented would be competitive, but for instance the authors of [75] seem to think that something like it might be. Several authors indeed claim that numerical evaluation of Bessel functions for large arguments and high order is difficult or inefficient, but we do not believe this statement: in our experience the standard recurrence method performs well. There do seem to be numerical instabilities in some Bessel series that can cause difficulty, although these can be mitigated by careful choices among the expansions [85]. We have not done a detailed comparison of methods.
In general one has to be a bit careful with rounding error if the recurrence relations are used to compute the Fourier coefficients: the relations can be unstable, as noted by many authors. A practical solution is to use forward recurrence for the first few and backward recurrence for the rest. This seems to have been first advocated by Blanch [10]. Similar considerations apply if numerical eigenvectors are used (after all, the matrix is simply an arrangement of the recurrence relations). A rule of thumb is that if one is working to d digits with the forward recurrence, then the Fourier coefficients will decay down to about 10 −d/2 times the magnitude of the largest coefficient and then rounding error will start to impact the results thereafter. In our experiments we simply used ultra-high precision and didn't worry about rounding errors at all.

Computing a generalized eigenfunction.
At least four ways suggest themselves to compute the generalized eigenfunction u = ∂y/∂a at a double point q * with double eigenvalue a * , which is a solution of equation (4.7) that satisfies periodic boundary conditions. We duplicate that equation here for convenience (recall indicates d/dα): u + (a * − 2q * cos 2α)u + y * = 0 .
The first and simplest way is undoubtedly what people actually use: one pretends that the eigenvalue is not actually a double one-typically because of rounding error it would have split anyway into a * + d √ q − q * + · · · and a * − d √ q − q * + · · · where q is a floating-point approximation to q * anyway-and then use the computed eigenfunctions from the matrix method, each with norm O((q − q * ) 1/2 ) and simply live with the errors. That does not sound like professional practice, but if it is done knowingly then we suspect that it will usually give perfectly reasonable answers. If done unknowingly then we disapprove, but the criminals will likely get away with it.
The second way is to compute a generalized eigenvector of the infinite tridiagonal matrix A for the problem. This is scarcely harder than the crude approach above: first, one averages the two computed approximate eigenvectors for the double eigenvalue split pair, and averages the eigenvalues, to get a more accurate double eigenvalue a * and eigenvector v * . For convenience in the exposition below, suppose that the eigenvectors of A are numbered v 1 (corresponding to eigenvalue a 1 ), v 2 (corresponding to eigenvalue a 2 ), and so on. Suppose further that we have reordered the eigenvalue and eigenvector pairs so that the pair of eigenvalues that arose on splitting the double eigenvalue numerically are put in positions 1 and 2, and similarly put their associated eigenvectors in columns 1 and 2 of the matrix of eigenvectors, and number them v 1 and v 2 . Put a * = (a 1 + a 2 )/2 and v * = (v 1 + v 2 )/2. Averaging is well-known to produce good estimates of double eigenvalues and eigenvectors. If, of course, by some miracle the eigenvalue routine doesn't split the double eigenvalue and instead produces a single, accurate, double eigenvalue and only one corresponding eigenvector, then we use that.
This system is singular, but v * is in its range and this is not difficult; one could use the SVD, for instance 22 . Since the matrix is of low dimension, the expense of the SVD is no obstacle. Then one uses the entries of u to construct the generalized eigenfunction u(z) in the same manner one constructs the eigenfunction v 1 (z) from the eigenvector v * .
To expand a given function f (z) as a sum of these eigenfunctions, notice that the "norm" of v 1 is zero, but the bilinear form of v 1 (z) with u(z) is nonzero. The generalized eigenfunction and v 1 (z) are each orthogonal to all other eigenfunctions, however. Put Then the γ k are easily found by orthogonality as usual: We remind you that this bilinear form does not use the conjugate. It is not an inner product. The "norm" of a nonzero function can be negative, complex, or indeed zero. Such "norms" are called indefinite norms, in for instance [4].
To find α and β we use the fact that while the norm of v 1 (z) is zero, the norm of the generalized eigenfunction u(z) is not zero, and also p 0 v 1 (z)u(z) dz = 0 .

Thus the two equations
give us a triangular two-by-two system (indeed with constant diagonal) to solve for the unknown coefficients.
For example, consider q = 1.468768613785142 i, the Mulholland-Goldstein double point again, and its associated eigenvalue a = 2.08869890274970 (computed this time by averaging the computed eigenvalues of the N by N matrix, where we took N = 25: its split eigenvalues were 2.088698902749696 72 27 ± 8.31667446021810974 · 10 −8 i), where we have displayed the distinct real digits in red. We averaged the corresponding eigenvectors to get v * = [c 1 , c 2 , . . . , c N ], and put v 1 (z) = c 1 √ 2 + 20 k=2 c k cos(2(k − 1)z) . 22 It is straightforward to set up the recurrence relations for this generalized eigenvector: they are merely forced versions of the recurrences in equations 3.18-3.30; but as mentioned these recurrences are known to be unstable sometimes, while in contrast the SVD has the virtue of answering the question in a manner that relieves all doubts.

The
√ 2 is needed because the symmetrizing trick for the matrix of equation (3.27) gives an extra √ 2 in the 0th coefficient. We then enforced v 1 (0) = 1 by scaling. We index from 1 in the above equation because that is usual for matrices. The real and imaginary parts of this are plotted in figure 11(a); v 1 (z) is, of course, an approximation for ce 0 (q, z) = ce 2 (q, z), the coalesced eigenfunction. By examining its residual v 1 + (a − 2q cos 2z)v 1 we see that it is accurate to 10 −14 (plot not shown). We verified that v 1 (z) has numerically zero norm, as well: We then use the SVD to solve the singular system (a * I − A)u = −v * , and form We then removed a multiple of v 1 (z) so that u(0) = 0. The real and imaginary parts of this generalized eigenfunction are plotted in figure 11(b). We now use these eigenfunctions to expand a smooth function of period π. We took as an example (5.12) f (z) = e cos 2z cos 6z ≈ αv 1 (z) + βu(z) + 20 k=3 γ k v k (z) .
As described above, we computed the coefficients of eigenfunctions v k (z) for k = 3, 4, . . ., 20 by orthogonality. The final five eigenvalues and eigenvectors were not needed.
Then because we may identify β ≈ 0.3152 + 0.1086 i . Now because where both of the integrals on the right are nonzero and we now know β, this gives α ≈ 0.1536 − 0.08560 i. As you can see from the graph of the magnitudes of the computed γ k for 3 ≤ k ≤ 25 in figure 12(a), these are appreciable. Notice also the exponential decay of the coefficients; this is typical for a spectral expansion of a smooth function. We plot the error f (z)−αv 1 (z)−βu(z)− k≥3 γ k v k (z) in figure 12(b) where we see that it is less than 5 · 10 −14 .
The third way that we were thinking that one could compute these eigenfunctions seems a little more complicated: we solve the initial-value problem sweeping forward for v 1 (z) using the Hermite-Obreschkoff method described earlier, which chooses the mesh; we record the local Taylor series for the two local functions satisfying y(z j ) = 1, y (z j ) = 0 and y(z j ) = 0, y (z j ) = 1. We then solve the boundary value problem for u(z) on that interval by imposing periodic boundary conditions and using collocation at two points in each interval. If there were M subintervals, this gives an almost block diagonal matrix 23 of 2M equations in the 2M unknowns, namely the coefficients α k . Top row left: Real and imaginary parts of the coalesced eigenfunctions v 1 (z) = ce 0 (q, z) = ce 2 (q, z) corresponding to the Mulholland-Goldstein double point q ≈ 1.4688 i (real part in black, imaginary part in red). On the right, we have the corresponding generalized eigenfunction obtained by solving y + (a − 2q cos 2z)y + v 1 = 0. Bottom row left: Real and imaginary parts of the coalesced eigenfunctions v 1 (z) = se 2 (q, z) = se 4 (q, z) corresponding to the next-largest pure imaginary double point q = 6.92895 . . . i with eigenvalue approximately 11.1905. On the right, we have the corresponding generalized eigenfunction obtained by solving y + (a − 2q cos 2z)y + v 1 = 0. and β k of the linear combination of the two solutions at each interior node, together with α 0 = α M and β 0 = β M . This sounds involved, but in fact it is straightforward. We suspect that few people will actually implement this method, however; we haven't, yet, either.
A fourth way, similar but perhaps even simpler, which we actually did do, is to use the Green's function from equation (5.4); we already have methods for integrating strings of blends and for multiplying strings of blends, so all this requires is the ability to transfer w I and w II onto the same string of blends. We found it simplest just to compute them at the same time. The general solution of equation (  the dependence of w I (z; a, q) on a and q and similarly w II for brevity) If we are solving for a generalized eigenfunction for ce 2g (q * , z) then y * = ce 2g (q * , z) = w I (z; q * , a * ) and this is already periodic so α = 0 and therefore u(0) = 0; notice that w II (0) = 0; moreover the integral to π for w I (ζ)y * (ζ) is also zero, so that at z = π 0 = u(π) = βw II (π) − w I (π) π 0 w II (ζ)y * (ζ) dζ .
Since w I (π) = 0 because it, being the eigenfunction in question in this example, is periodic, we see that β = 0 as well. Thus all we need to identify the eigenfunction is the integral against the Green's function. We tried this, and it confirmed the results in the top row of figure 11. The bottom row in that figure was computed this way.
6. Concluding remarks. We have presented a historical survey of the computation of the Mathieu functions, which are defined as the period π and period 2π solutions of the Mathieu differential equation (3.1). Our original motivation for this undertaking was to solve a problem of pulsatile blood flow in a vessel of elliptic cross-section. To actually do that we used the spectral method, but with the brute force of multiple precision to power through the double-eigenvalue difficulty. That was inelegant, so we investigated the alternatives.
To carry out this investigation, we implemented our own procedures (in Maple) for the computation of Mathieu functions (and generalized eigenfunctions) in order to explore some of the difficulties involved. Our code is at present, like we imagine Blanch's to have been, "artesanal" and was intended only for use with careful human supervision. For experimentation, of course, this is a feature, not a bug. The task of constructing fully general, bulletproof code for the computation of Mathieu functions is one that calls for dedicated effort and analysis. We hope that this present paper encourages a team to undertake the task. We know of no such code in existence currently.
In 1957 G. Temple called Mathieu functions "indispensable but intractable instruments of mathematical physics". The word "intractable" has a technical meaning nowadays, and certainly computing Mathieu function is not intractable in that sense. But they are somewhat involved, and the major remaining question is how best to compute the sometimes-needed generalized eigenfunctions. Of course, we must not forget about the advances in direct numerical solution of the underlying PDE: the methods of [36] may be preferable in many applications over expansion in Mathieu functions.
Mathieu's work on the nodal lines of an elliptic drum occurred a little after Ernst Chladni (1756-1827) demonstrated the existence of nodal lines on freely vibrating plates, apparently following work of Hooke; for a lovely discussion of these see [35].
Perhaps the most interesting results of our readings of the literature to us were, first, the discovery that Mathieu had anticipated Lindstedt's anti-secularity perturbation method by over a decade 24 ; and that Blanch's algorithm for computing eigenvalues could be implemented (easily!) in series, thereby allowing one to compute series for eigenvalues which would allow greater surety in continuation methods (Newton iteration starting with estimates from nearby q is usually used) or to compute Puiseux series about double points. We believe that the computation of these Puiseux series is new to this paper. The notion of generalized eigenfunctions for the Mathieu equation is not new to this paper (it is in [60], as previously stated) but we believe that the details of their computation are presented for the first time here.
where A and B are arbitrary constants and J 0 and Y 0 are Bessel functions of order zero and of the first and second kind, respectively, satisfying the standard Bessel equations The new variable ζ is related to the radial coordinates by where Ω is a frequency parameter related to the nondimensional frequency parameter In the case of a tube of elliptic cross section the boundary conditions suggest a transformation to elliptic coordinates Using the translation The foregoing analysis demonstrates clearly that the boundary and boundary conditions play a key role in the transition from Bessel equations to Mathieu equations in the case of pulsatile flow in a tube. However, while this explains the mathematical aspects and origin of the transition from Bessel to Mathieu equations in this particular application, it does not provide a hint as to the corresponding origin of this transition in the context of the physics of the flow. One of the aims of our study was to understand the link between the mathematical and the physical aspects and origin of this transformation.
While, at first, properties of the flow in a slightly elliptic tube appear as a small perturbation of the flow in a circular tube (and can indeed be analyzed in this manner [27]), a closer look at the shear stress on the boundary of the tube shows a more significant change that occurs as the cross section of the tube changes from circular to elliptic as shown in Figure 13. The figure shows that no matter how small the change from circular to elliptic cross section is, the two vertices 25 of the ellipse lead to periodicity in the distribution of shear stress on the boundary. This periodicity is the origin of the transition from Bessel to Mathieu equations in the description of the flow. Fig. 13. Shear stress at the boundary of a tube in pulsatile flow. The curves represent a range of tubes with different cross sections ranging from circular to increasingly elliptic. The numbers on the right identify the tubes in terms of the ratio of their minor to major axes of their cross section. The y-axis on the left of the figure represents that ratio of the shear stress normalized in terms of the shear stress in a tube of circular cross section.
We start as usual with the supposition that we have two separate eigenfunctions, y m (x) and y n (x), of Mathieu's equation, corresponding to different eigenvalues λ m and λ n : y m (x) − q cos(2x)y m (x) = λ m y m (x) y n (x) − q cos(2x)y n (x) = λ n y n (x) .

(B.2)
Here either of the eigenvalues λ m and λ n can be an a k (q) or b k (q), so long as they are different from each other. Multiply the first equation by y n (x) and the second by y m (x) and subtract to get Since the eigenvalues are distinct, this ensures that the bilinear form is zero.
If we instead tried to use the genuine inner product u, v := 2π 0 u(x)v(x) dx and norm u 2 = u, u 1/2 , instead of the bilinear form and indefinite "norm," we fail. The Mathieu functions, when q is not real, are not orthogonal with respect to this inner product, as can be verified by a straightforward computation with some non-real q, say q = 1.0i. Computing the eigenvalues a 2 (1.0i) and a 4 (1.0i), which indeed have distinctly different numerical values, and numerically computing the integral 2π 0 ce 4 (x, 1.0i)ce 2 (x, 1.0i) dx , we get a complex number of magnitude about 0.5138. That inner product is not zero. In contrast, we find that the bilinear form (without the complex conjugate) is zero to numerical accuracy, as it should be.
Where does the proof fail, when q is not real? In trying the method of the proof, we would first have to conjugate one of the equations to start with; say the second one.
y m (x) − q cos(2x)y m (x) = λ m y m (x) y n (x) − q cos(2x)y n (x) = λ n y n (x) . (B.5) Since we will integrate from x = 0 to x = 2π again, x can be taken as real. Now multiply the first equation by y n (x) and the second by y m (x) and subtract: but now, the linear term does not cancel-unless, of course, q is real: (B.6) y n (x)y m (x)−y m (x)y n (x)+cos 2x (q − q) y m (x)y n (x) = (λ m −λ n )y m (x)y n (x).
Without that cancellation, the integrals of the next step tell us nothing useful. And indeed as we saw above by direct computation for a specific non-real q, the naively desired result-namely, orthogonality under a true inner product instead of a bilinear form-isn't true, and so no possible proof could exist.

Appendix C. Confocal Ellipses verses Fixed Aspect Ratio Ellipses and
Bessel functions as a limit of Mathieu functions. To help people visualize the difference between a family of confocal ellipses, all of which have the same set of foci, and a family of fixed aspect ratio ellipses, which is likely what most people imagine when they think of a family of ellipses, we made figure 14. The point is to demonstrate how quickly confocal ellipses become circular-looking.
If ε is small, we recognize this as a small (admittedly singular) perturbation of a form of Bessel's equation, so the outer solution, away from ξ = 0, will be y(ξ) = C 1 I M (ξ) + C 2 K M (ξ) + O(ε 2 ) ,  . The circle (red dashed line) was chosen to have radius r = c exp(2)/2, which as the mean of the semi-major and semi-minor axes of the largest ellipse, corresponds well. We see that the largest confocal ellipse shown here is appreciably circular. Larger β would generate even more circular-looking ellipses. In contrast, the fixed aspect-ratio ellipses all have the same aspect ratio as the smallest nonsingular confocal ellipse shown on the left, namely a/b = tanh(2/3) ≈ 0.58.
where I and K denote Bessel functions. The next term, which is O(ε 2 ), in the perturbation expansion is somewhat complicated, though small in size when we computed it, but in fact we don't need it. The result is plotted in figure 15 and compared with the relevant Mathieu function. Now, as ε → 0, for ξ to remain O(1) we must have sin(z) → ∞; that is, this matching will only be valid for large imaginary z. This is the double limit mentioned earlier. Owing to the exponential growth of sin, however, it won't have to be that large. We already saw in figure 14(a) that with c = 1/3 the confocal ellipse is pretty circular already for z = 2; and indeed already by q = 1/10 and for ce 0 (q, iη) (that is, we use a 0 = −0.00495 . . .) we see a marked resemblance to (I M (ξ)) in figure 15. Indeed if we choose a point (η = 2.5) and choose a normalization factor C = 0.58698 so that the asymptotic approximation agrees at that point then the error ce 0 (q, iη) − C (I M (ε sin(iη))) gets very small very quickly as shown in figure 15(b). Notice that here ε ≈ 0.632, which isn't even all that small. The various complex numbers in this example have shuffled the usual behaviour amongst the Bessel functions, and we have simply taken a real part for comparison; but the real part of the solution to a real linear differential equation, in this case the modified Mathieu equation, is also a solution.
Appendix D. Comparing Mathieu's perturbation series to Maple's. In order to compare the q-series produced by Maple's series command with the hand computation of Mathieu, we have to ensure that the same normalization is enforced. We also had to check our work carefully to ensure that the Maple results were correct. Since Mathieu's normalization was to make the coefficient of cos gα equal to unity, and Maple's normalization is to make p 0 ce 2 g (α) dα = π except when g = 0 when it is 2π, we compute Maple's series and then divide by the coefficient of cos gα. Similarly and M ≈ 0.4527i. Here z = iη is purely imaginary. We have ignored the O(ε 2 ) term of the outer solution. Even so, the match seems good.
for the se g (α) functions. By independently computing residuals, we assure ourselves of the correctness of those computations. The Maple code for this appendix can be found at Rob Corless's GitHub Repository. For ce g (q, t) with symbolic g we get, on dividing by the coefficient of cos gt, the following terms in the series. The coefficient of q is, provided g > 1, (D. 1) cos (g − 2) α 4 (g − 1) − cos (g + 2) α 4 (g + 1) .
Finally, he has an extra factor (g 2 − 4) in the denominator of the q 6 (h 12 ) term of the eigenvalue and an apparently incorrect numerator, 9 g 5 + 22 g 4 − 203 g 2 − 116. But if you replace the 5th power with a 6th (and, really, reading the PDF of this manuscript, it's hard to tell whether it should be a 6 anyway), this factors into the correct form: 9 g 6 + 22 g 4 − 203 g 2 − 116 = (g 2 − 4)(9 g 4 + 58 g 2 + 29). Actually, on the line above Mathieu's final form for R, the power is more clearly a 6: this isn't even a typo, just something hard to read given the printing and transcribing process.
The end of the story is that Mathieu was able to give the correct generic perturbation series to ce g (h 2 , α)/F up to and including terms of order q 5 (h 10 ), on the understanding that the factor F was chosen to make all coefficients of cos gα equal to zero apart from the first one. Special series. The generic series is good only for "large enough" g. For specific small g, and indeed for any fixed g if one wants to compute enough terms, special computations have to be done. We here give the results that Mathieu attempted, and comment on some minor errors in his paper.
which is nearly the same as Mathieu's, (D.12) R = 16 + 1 30 h 4 + 433 864000 h 8 − 189983 21772800000 h 12 + · · · except he erroneously reports 189983/21772800000 as the coefficient of q 6 . This is an irreducible fraction and not equal to the correct coefficient.