Computation and applications of Mathieu functions: A historical perspective
CCOMPUTATION AND APPLICATIONS OF MATHIEU FUNCTIONS:A HISTORICAL PERSPECTIVE ∗ CHRIS BRIMACOMBE † , ROBERT M. CORLESS ‡ , AND
MAIR ZAMIR ‡ Abstract.
Mathieu functions of period π or 2 π , also called elliptic cylinder functions, wereintroduced in 1868 by ´Emile Mathieu together with so-called modified Mathieu functions, in orderto help understand the vibrations of an elastic membrane set in a fixed elliptical hoop. Thesefunctions still occur frequently in applications today: our interest, for instance, was stimulated bya problem of pulsatile blood flow in a blood vessel compressed into an elliptical cross-section. Thispaper surveys and recapitulates the historical development of the theory and methods of computationfor Mathieu functions and modified Mathieu functions and identifies some gaps in current softwarecapability, particularly to do with double eigenvalues of the Mathieu equation. We demonstratehow to compute Puiseux expansions of the Mathieu eigenvalues about such double eigenvalues, andgive methods to compute the generalized eigenfunctions that arise there. In examining Mathieu’soriginal contribution, we bring out that his use of anti-secularity predates that of Lindstedt. Forinterest, we also provide short biographies of some of the major mathematical researchers involvedin the history of the Mathieu functions: ´Emile Mathieu, Sir Edmund Whittaker, Edward Ince, andGertrude Blanch. Key words.
Mathieu functions; modified Mathieu functions; historical survey; computation ofMathieu functions; double eigenvalues; frames; Puiseux series
AMS subject classifications.
1. Introduction.
What is the sound of an elliptic drum? In a memoir presentedat the Sorbonne in 1868, ´Emile Mathieu showed the way to find the answer, when hedescribed the solution of a mechanical vibration problem characterized by an ellipticboundary. The memoir was groundbreaking in that it introduced a new differentialequation whose eigenvalues and corresponding periodic solutions led to the definitionof a new class of functions. In 1912 Whittaker named these new functions in honour oftheir discoverer: the differential equation is now the Mathieu equation and the π and2 π periodic solutions (and only these periodic solutions) are the Mathieu functions.While work on theoretical and analytical aspects of Mathieu functions has con-tinued since the introduction of these functions in 1868, the focus has shifted in recentyears to work on numerical and computational aspects. This development, driven bysteep advances in digital technology, has given rise to heavy reliance today on “pack-aged software” for the evaluation of Mathieu functions. This practice comes with therisk of concealing as yet unresolved (or at least not completely resolved) analyticaland computational issues involved in the use of Mathieu functions.We review methods of computation of periodic Mathieu functions and claim thatall current published codes fall short at some “exceptional” or “double” points, asnoted both by [15] and by [40]. Although those objections were made more than fiftyand forty years ago respectively, and theoretical work on these was completed in [49],we believe that there still is no fully satisfactory code available, as we will detail. Wewill discuss what to do in practice at these double points, where the eigenvalues merge ∗ Submitted to the editors DATE.
Funding:
This work was supported by NSERC and, while RMC was visiting the Isaac NewtonInstitute during the programme Complex Analysis: Tools, techniques, and applications, by EPSRCGrant † University of Toronto, Toronto,
Canada ([email protected] ). ‡ Department of Applied Mathematics, Western University, London,
Canada ([email protected],[email protected]). 1 a r X i v : . [ m a t h - ph ] A ug C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR and we lose one independent eigenfunction and thus lose completeness of our set ofeigenfunctions. In modern terms this suggests using a frame containing the periodicMathieu functions as a method of solving partial differential equations (PDEs) withelliptic boundaries. We briefly explore one such frame, but only as a pointer to futureresearch.While we cover a lot of ground, we also point out several important referencesthat contain many details (and indeed many important results) that we omit. Ourfirst and most important reference is to Chapter 20 of the classic [1], written byGertrude Blanch, and its successor, Chapter 28 of the Digital Library of MathematicalFunctions (henceforth DLMF, https://dlmf.nist.gov/28.1) written by G. Wolf. TheDLMF is intended to be an updated replacement for [1] in this computer age. Bothare (perfectly legally) free online. The DLMF gives pointers to some, but not all,of the alternative notations used elsewhere in the rather substantial literature. Theoriginal Chapter 20 of [1] has even more references.
In section 2 we briefly outline three applica-tions 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 func-tions. 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 correspondingeigenvalue, and Puiseux series for the eigenvalue about those points. In section 5we look at algorithms for computing solutions of the Mathieu equations, includingMathieu functions (eigenfunctions), together with existing specialized algorithms justfor the 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 6. In appendix A we provide more details of theapplication that motivated us to undertake this work. In appendix B we discuss con-focal ellipses, and give a singular perturbation argument relating Mathieu functionsto Bessel functions in the limit as ellipses become circles. Finally, in appendix C wecompare Mathieu’s perturbative solution to the Mathieu equation with that producedby a computer algebra system, and apart from minor errors and typos in his paperconfirm his results.
2. Applications.2.1. Columns and Strings Under Periodically Varying Forces.
One of themore surprising physics experiments is the stabilization of a vertical hinged columnby up-and-down vibration. There are many YouTube videos demonstrating this, andafter finishing this paper the reader may choose to search some out: a good set ofkeywords is “stability of the inverted pendulum”. This phenomenon is explained usingtopological terms in [46]. The mathematics of it have been known at least since [71],and are concerned directly with the Mathieu equation if the periodic force is a simplesine or cosine.There are similar problems that while even more complicated still need the Math-ieu equation. We will briefly describe the model studied in [47], namely the stabilityof a column under periodically-varying compression or a string under tension witha periodically-varying tension force. These problems have infinitely many degreesof freedom, as opposed to the simple inverted pendulum above. Take the case ofa column under compression. If the force F is greater than the Euler load, thenit is well known that the column can buckle; the question that is at issue here is if F = P + H cos ωt , consisting of a steady part plus a periodically-varying part, can one ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE Fig. 1 . A flexible string with a time-varying tension F ( t ) . The figure indicates that the defor-mation w ( t ) is confined to the plane. We model this figure on Figure 1 of [47]. choose H and ω so that, even if P is larger than the Euler load, the column remainsstable?The answer is a qualified yes; one qualification is that at least part of the time, theforce must be less than the Euler load, which makes sense. The other qualificationsrequire the study of the stability of the solutions to the Mathieu equation, and requiresa good knowledge of the periodic solutions of the Mathieu equation. That paper alsostudies the motion of a string with a periodically-forced tension with a similar model;see figure 1. The equations of motion they derive are EI ∂ w∂x − F ( t ) ∂ w∂x + m ∂ w∂t , where they take F ( t ) = P + H cos ωt and the familiar Young’s modulus E characterizesthe relationship between stress and strain, and I is the moment of inertia of the body’scross-section.“In all of these problems the Mathieu equation (more properly, a sequence ofMathieu equations in the continuous systems) plays a central rˆole, since the decisionas to stability depends on the character of the solutions to such equations” [From theintroduction in [47].]In fact the study of the stability of these systems requires a little bit more thanwe are going to cover in this paper: it needs the Floquet theory, which we only lightlytouch in this paper, in addition to the study of the periodic solutions, where we spendmore time. Our motivation forthis present paper originated from a problem in pulsatile blood flow. Under normalcircumstances, blood flow occurs in vessels of circular cross sections, but under anumber of important pathological conditions the vessels are deformed by externalforces 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 crosssection as a simple, mathematically tractable, departure from a circular cross section,with the main focus being on the hemodynamic consequences of that departure. Inthe present paper, in contrast, our focus is on the mathematical and computationalconsequences.Fluid flow in a tube (see figure 2) is in general governed by a simplified form ofthe Navier-Stokes equations [75]. 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 axisof the tube. If the flow is pulsatile , as in the cardiovascular system, the velocity U consists of a steady part u plus an oscillatory part u ( t ) such that(2.1) U = u + u ( t ) , where t is time. Only the oscillatory component of velocity, u ( t ), is relevant to thepresent discussion since the transition from Bessel to Mathieu functions occurs in C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR
Fig. 2 . Flow in a tube of elliptic cross-section. The elliptic cross-section is meant to modelsome pathologies where the blood vessels are deformed so that they are no longer circular. Thisfigure is modelled after Figure 3.1, page 45, of [75].
Fig. 3 . Contours and nodal lines of a possible pure vibration node of an elliptic drum. Theaspect ratio of the ellipse is . The rightmost loops have positive height, and adjacent cells haveopposite signs. Notice the hyperbolic nodal lines where the membrane does not move, separating thecells. We have suppressed here the details of which Mathieu functions were used to produce thisfigure (we give them in section 3.1.5), but we followed the method outlined by Mathieu in his 1868paper [48]. the governing equations of this component of the flow as the cross section of thetube changes from circular to elliptic. We give a brief sketch with more details ofour motivating application in appendix A; it will be more understandable once basicnotions are introduced. In that appendix the focus is on the transition from Besselto Mathieu equations in the case of pulsatile flow in a tube as the cross section of thetube transitions from circular to elliptic geometry.While pulsatile flow in tubes of circular cross sections is governed by Bessel equa-tions, with solutions in terms of Bessel functions, the corresponding situation in tubesof elliptic cross sections involves the Mathieu equation and Mathieu functions. Froma mathematical perspective, this transition is not only a matter of curiosity but also amatter of practical importance because of the comparative difficulties involved in thenumerical computation of these two functions. The difficulties and pitfalls involvedin the computation of Mathieu functions are the focus of much of this present paper.
ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE One of the simplest phys-ical problems whose solutions involve the Mathieu equation and the Mathieu functionsis the sound made by an elliptical drum. This in fact was the problem Mathieu himselfstudied, and in section 3.1 we will look in detail at how he solved it. The physicalproblem being modelled is, in fact, remarkably easy to define: imagine an ellipticalhoop, fixed immovably, and a thin homogeneous membrane stretched tight across thehoop, making a drum. What are the natural modes of vibration of this drum, andhow could we describe mathematically its motion once struck?One possible natural mode for one particular drum is pictured in figure 3. Therewe see contours of vibration, including contours where there is no motion (“nodallines”). This figure was drawn using our own software to compute the relevant Mathieufunctions, but many software packages exist which could do this.
3. Historical overview, introducing notions and notation.
We introducethe Mathieu equation and the Mathieu functions in historical order, by discussingthe contributions of several of the main researchers involved. The result is a tour ofseveral aspects of late nineteenth-century and early twentieth-century mathematics.We also give some biographical details of these main figures. For reference, here isthe Mathieu equation in one common modern notation.(3.1) d ydx + [ a − q cos(2 x )] y = 0 , The parameter q is given by the physics or the geometry; the eigenvalue a must becalculated 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):(3.2) d ydz − [ a − q cosh(2 z )] y = 0 . The even solutions are conventionally written as Ce g ( z, q ) = ce g ( ± iz, q ). The eigen-values for even solutions are conventionally written a g ( q ). The odd solutions arewritten similarly, e.g. Se g ( z, q ) and Se g ( z, q ), and the odd eigenvalues are conven-tionally written b g ( q ). Here g is a nonnegative integer, and the solutions split intofurther classes if g is itself even or odd, as we will see. The use of the letter g for aninteger contradicts the usual I − N convention in use nowadays, from Fortran ; butMathieu used the letter g in this way and we find it convenient when its parity is asyet unspecified.If we write the general solution of the Mathieu equation (with no initial or bound-ary conditions applied) as y ( x ) = αw I ( x ; a, q ) + βw II ( x ; a, q )(3.3)for the pair satisfying w I (0; a, q ) = 1 with w (cid:48) I (0; a, q ) = 0 and w II (0; a, q ) = 0 with w (cid:48) II (0; a, q ) = 0, (using the notation of the DLMF) then the solution to the ModifiedMathieu equation must be y ( z ) = αw I ( ± iz ; a, q ) + βw II ( ± iz ; a, q ) . (3.4)Some software packages—for instance, Maple—denote these functions C and S re-spectively. When a ( q ) or b ( q ) is an eigenvalue, the periodic Mathieu functions must C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR satisfy (now allowing x to be complex and renaming it z )ce m ( z, q ) = c m w I ( z ; a m ( q ) , q )(3.5) se n ( z, q ) = s m w II ( z ; b n ( q ) , q ) , (3.6)for some normalization constants c m and s m . In this paper we take those nor-malization constants to be . In theory, the use of Mathieu functions and modified Mathieu functions in thesolution of the aforementioned physical problems is attractive because the functionsare analogous to harmonic functions, and expansions in terms of them can be efficientin comparison with direct numerical solution of the PDE model. In practice, thereare annoying difficulties: the available software might be restricted to real arguments,or use a different normalization than the one you are thinking of (the authors of [30]claim there are at least three normalizations in common use; they themselves use thesame normalization that we do here), or fail to be accurate for “difficult” values of theproblem parameters, say for large values. A more serious problem is the approximationproperties of the expansion itself, followed by the numerical stability of the expansion.The first question is explored in [64] for real q , with all the power of the Sturm-Liouvilletheory (which indeed Mathieu himself used in his 1868 paper).Now that we have sketched where we want to go with Mathieu functions, let usrecapitulate their development. Mathieu began his 1868 discus-sion [48] of the vibrations of an elastic membrane held fixed by a hoop in the shape ofan ellipse by first considering the simpler problem when the hoop is, in fact, circular .Mathieu’s discussion of the circular case starts with the PDE d wdt = m (cid:18) d wdx + d wdy (cid:19) (Mathieu uses d and not our modern Russian ∂ for partial derivatives) and thentransforms x = r cos α , y = r sin α to polar coordinates to get d wdt = m (cid:18) d wdr + 1 r dwdr + 1 r d wdα (cid:19) . Thereafter Mathieu uses what is now the standard method of separation of variablesfor a pure oscillation w = sin(2 λmt ) u ( r, α ) and u ( r, α ) = P ( α ) Q ( r ) to give a harmonicequation for P ( α ) and what we now call Bessel’s equation for Q ( r ): r d dr Q ( r ) + r ddr Q ( r ) − (cid:0) n − λ r (cid:1) Q ( r ) . One then writes the solution to the original vibration problem as a linear combina-tion of products of these eigenfunctions, and determines the unknown coefficients bymatching to the boundary conditions using orthogonality.Although Bessel (1784–1846) did his work first (and although it was actuallyDaniel Bernoulli who first identified this equation, even earlier) Mathieu does notcall this Bessel’s equation, or give it a name at all, but merely solves it in series inwhat must have been common practice at the time. He then demonstrates that one This memoir was translated from its nineteenth century French for us by Dr. Robert H. C. Moir,and the translation—which we believe may be of interest on its own—will be made available on arxiv.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE zeros of the various (series for the) Bessel functions, and cites Bourget’sMemoirs for a method for doing so, in order to identify the eigenvalues and eigenmodesof the vibration of the membrane in the circular case.Following an exactly similar strategy in the elliptic case, but using a confocalelliptic transformation x = c cosh( β ) cos( α ) and y = c sinh( β ) sin( α ) where 2 c isthe distance between the foci of the ellipse , Mathieu arrives at an equation that onseparation gives two equations equivalent to those now known as the Mathieu equation and the modified
Mathieu equation: − P d Pdα + 4 λ c cos ( α ) = + 1 Q d Qdβ + 4 λ c cosh ( β ) . “Comme le premier membre ne peut renfermer de α , et le second que β , ils sont´egaux a une mˆeme constante N .” Readers may be pleased to learn that separation ofvariables reads the same mutatis mutandis in the 21st century as it did in the 19th.There is one remaining difference to the modern notation, and that is the use ofcos α and cosh β . Using double-angle identities, Mathieu later in this same papertransformed these to something that we write in modern notation as d Pdα + ( a − q cos 2 α ) P = 0(3.7) d Qdβ − ( a + 2 q cosh 2 β ) Q = 0 . (3.8)Notice that these two equations can be transformed into each other by the change ofvariable β = iα . Here q = λ c is a relabeling of the parameter that contained thephysics, 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 earliervariable N by changing from cos α = (1 + cos 2 α ) / β = (1 + cosh 2 β ) / h where we have q here, and that notation is still occasionally used.The boundary conditions of the original problem reflect the elliptic geometry. Theangular coordinate α runs from 0 to 2 π , requiring periodicity (or, of course, from − π to π ). Often the problem is taken to have π -symmetry, reflecting a symmetry betweenthe top and the bottom of the ellipse. Mathieu noted thatthe theory of Jacques Charles Fran¸cois Sturm (1803–1855) implied that (for real q )what we now call the Sturm-Liouville form of the Mathieu equation could be writtenas ddx (cid:18) L dydx (cid:19) + G ( q ; x ) y = N y ( x ) , (3.9) Mathieu initially spelled the hyperbolic functions out explicitly, and later used “Pour simplifier”the notation E ( β ) = ( e β + e − β ) / β ) and E ( β ) = ( e β − e − β ) / β ). Johann HeinrichLambert had introduced the notation that we use today already in the 18th century; it is interestingthat it was not universally used in the 19th. Notice that this coordinate transformation is singular if c = 0. Therefore the connection ofMathieu functions to Bessel functions, while present, requires art to tease out. See Appendix B. Sturm was the successor of Poisson in the chair of mechanics in the Facult´e de Sciences, Paris.This fact has a certain poignancy when read together with the obituary of Mathieu [26], where welearn that Mathieu had desired that chair.
C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR (here, trivially, L = 1 and G ( q ; x ) = − q cos 2 x , while the eigenvalue N has beenmoved to the right hand side) and that several useful properties naturally followed.First, there are a countable number of eigenvalues N = − a k which can be arranged insequence a ( q ) < b ( q ) < a ( q ) < b ( q ) < a ( q ) < · · · , with the eigenvalues tendingto infinity.Strictly speaking, Mathieu had to extend the now-classical theory of Sturm to thecase of periodic boundary conditions (see e.g. [60] for a summary of this classical the-ory) in order to establish this much; more, he showed how to compute the eigenvaluesby an interesting perturbative argument which we will take up in section 3.1.3. Beforethat, he established an inequality around each neighbourhood of g , the square of aninteger; to us, it seems a convincing argument for the existence of each eigenvalue(and thus of the infinite collection).Next, to each (real) eigenvalue there corresponds a unique eigenfunction. Theseeigenfunctions are now called Mathieu functions, and they come in four classes. Ifthey are even and period π they are denoted by ce k ( x ), for k ≥
0. If they are evenand period 2 π they are denoted by ce k +1 ( x ), for k ≥
0. If they are odd and period π they are denoted by se k ( x ) for k ≥
1. If they are odd and period 2 π they are denotedby se k +1 ( x ) for k ≥
1. See figure 4 for a representative graph of a few low-frequencyMathieu functions, with q = 1 .
5. Mathieu established that these eigenfunctions areorthogonal with respect to the inner product defined by (cid:104) y k , y (cid:96) (cid:105) := (cid:90) p y k ( x ) y (cid:96) ( x ) dx = const · [ k = (cid:96) ] . (3.10)Eigenfunctions of one class are orthogonal to eigenfunctions of the other class withthe same period, using this inner product. Here p = π or p = 2 π as appropriate.We have used the Iverson convention [ k = (cid:96) ] to mean 1 if k = (cid:96) and 0 otherwise.The exact method used for normalization is a matter of convention. Finally, theseclasses of orthogonal eigenfunctions are complete : every reasonable even/odd period- p function f ( x ) can be expanded in a convergent series f ( x ) = (cid:88) k ≥ α k y k ( x ) . (3.11)Periodic functions which are neither even nor odd can of course be written f ( x ) = f ( x ) + f ( − x )2 + f ( x ) − f ( − x )2as a sum of an even function and an odd function and thus will use both classes of thegiven period in its expansion. Since the coefficients α k can be determined by orthogo-nality, this series is expected to be practical and to enable us to find computationallyuseful solutions to the original PDE by matching the boundary condition at the edgeof the ellipse.We emphasize that Mathieu only established this for real q . The case of complex q is a different matter, as noted by [49] and by [15]. The Sturmian theory fails therebecause there may be (and in fact are) double eigenvalues, and in that case Blanchpoints out that (cid:104) y k , y k (cid:105) = 0 irrespective of normalization. We will return to this later. As with Fourier series, one can prove convergence for a quite wide class of functions; however,for utility and rapid convergence, one needs at least continuity, and the more smoothness the better.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE (a) Some ce k ( α ) for q = 1 . k ( α ) for q = 1 . k +1 ( α ) for q = 1 . k +1 ( α ) for q = 1 . Fig. 4 . The first few Mathieu functions when q = 1 . . The normalization shown here is that y (0) = 1 in the case of the even functions and y (cid:48) (0) = 1 in the case of the odd functions. Period π functions are shown in the top row, period π functions in the bottom row. Mathieu normalized the solutions of his equation in away that might seem curious to modern eyes. First, he noted that it is easy to seethat the solution y ( x ) of any linear second order ordinary differential equation (ODE)may be written as y = P + P , where P is either a maximum or a minimum at x = a while P is zero at x = a ( a is arbitrary). A moment’s reflection shows thatMathieu was correct, and that this is true for any a in the domain of definition of y ( x ).0 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR
What is curious is that this separates a second order equation into two parts, eachof which imposes only one condition: P (cid:48) ( a ) = 0, or P ( a ) = 0. Clearly any multipleof P or of P will satisfy the same condition, but only one combination of suchfunctions will equal y ( x ). This left Mathieu free to normalize his functions in any wayconvenient to him, and he took great advantage of it, in particular in his perturbativesolutions. To normalize his functions, he chose to make the coefficient of cos gα in itsseries expansion to be unity (and similarly the coefficient of sin gα in the case of oddeigenfunctions); see appendix C for a comparison with one modern normalization. AsInce pointed out later in [44], this is not always possible because for some values of q the coefficient of cos gα is actually zero; but it is at least almost always possible inthe 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 normaliza-tion discussed at the end of section 20.5 in [1], namely, y ( x ) = a ce g ( x ) + b se g ( x ) andwe specify that ce g ( x ) satisfy not only ce (cid:48) g (0) = 0 but also ce g (0) = 1, and similarlyse g (0) = 0 and se (cid:48) g (0) = 1. As Blanch states in the aforementioned reference, conver-sion between normalizations is “rather easy,” but we wanted one that would alwayswork.However, in many published papers and codes the norm using the inner prod-uct (3.10) is nominally enforced to be π (or 2 π ), which can (again) almost always bedone, but not universally: at double eigenvalues, the “norm” must vanish. This isa 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 meansthat in modern terms the norm is an “indefinite” norm [5], and requires some carein handling. To emphasize, in this paper we do not normalize by the inner product,but instead choose to enforce initial conditions as above; this is done elsewhere in theliterature, but not commonly. In the 1868 paperunder 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 or-der in h ). It is interesting to note that to do so he essentially used what we nowknow as anti-secularity : he chose series coefficients in the eigenvalue expansion inorder to eliminate secular terms in the expansion for the eigenfunction and therebyenforced periodicity of the solution. This notion is typically introduced nowadays asthe Lindstedt-Poincar´e method or instead using the method of multiple scales; in aneven more modern way, this can be done by renormalization . See for instance [21]and the references therein for more details of those methods.When Mathieu published his Memoire in 1868, Anders Lindstedt was only 14years old and his work on perturbation was many years in the future. Mathieu wouldhave good grounds for a claim to priority, even though (perhaps) Lindstedt’s workwas somewhat more general. Mathieu’s use of anti-secularity is clear, however, onceone tries to recapitulate his steps; it seems very natural, although Mathieu does notcomment on it explicitly. Indeed, his section 11 which details the perturbation solutionreads more like an informal summary of notes of how to proceed, with many detailsleft out. Nonetheless, using anti-secularity to enforce periodicity is exactly what hedid. He also made several elegant uses of his freedom to normalize in the problem inorder to reduce the labour involved. We have implemented his solution in a computeralgebra system, to recapitulate his steps and fill in details; some of our computationsare compared to his in Appendix C. ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE π solution of equation (3.7)when q = h was small and the eigenvalue a approached g , the square of an unspec-ified integer. The solution in his notation and with his normalization and to fewerterms than he calculated to is :ce g ( α ) = cos gα + (cid:18) cos ( g − α g − − cos ( g + 2) α g + 1) (cid:19) h + (cid:18) cos ( g − α g −
2) ( g −
1) + cos ( g + 4) α g + 2) ( g + 1) (cid:19) h + (cid:32) cos ( g − α g −
4) ( g −
2) ( g −
1) + (cid:0) g − g + 7 (cid:1) cos ( g − α g −
2) ( g + 1) ( g − − (cid:0) g + 4 g + 7 (cid:1) cos ( g + 2) α g + 2) ( g + 1) ( g − − cos ( g + 6) α g + 2) ( g + 3) ( g + 1) (cid:33) h + O ( h )(3.12)As he noted, this series is valid only for large enough integers g . He also correctlycomputed the corresponding eigenvalue (he called it R in this part of his paper) as(3.13) a = g + h g −
1) ( g + 1) + (cid:0) g + 7 (cid:1) h (32 g −
64) ( g + 2) ( g − ( g + 1) + · · · . Mathieu then goes on to show how to compute perturbation solutions for specific,smaller, frequencies g . See appendix C for details.The idea of a series expression for the Mathieu functions was, of course, naturalfor the time. Whether the idea of enforcing periodicity by expanding the eigenvaluein series was original to Mathieu, we do not know; but its presence in this papercertainly predates Lindstedt’s work.For Mathieu, q = h was real, and small (if the interfocal distance 2 c was small).In many modern applications, q might be complex, or large, or both. It took manyyears of further research by others to go beyond these series. After finding these pertur-bation expansions, Mathieu took a different tack: he changed variables, first with ν = cos α , whereupon the Mathieu equation becomes (equation 28.2.3 in the DLMF)(3.14) (1 − ν ) d Pdν − ν dPdν + ( a + 2 q (1 − ν )) P = 0and alternatively by ν (cid:48) = sin α , whereupon the Mathieu equation becomes(3.15) (1 − ν (cid:48) ) d Pdν (cid:48) − ν (cid:48) dPdν (cid:48) + ( a − q (1 − ν (cid:48) )) P = 0 . The DLMF gives yet another algebraic form, using the change of variables ζ = sin α .These are interesting for several reasons, and we will mention in section 3.2 some ofthe further properties that can be deduced from these equations. What Mathieu usedthem 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, Ince pointed out in [45] that for certain values of q the coefficient of cos gα in this expansioncould, in fact, be zero, and that therefore this normalization is impracticable. For generic values of a and q , however, it works. C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR of course, but an important difference is seen between the two forms: in the originalformulation, 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 inequation (5.3) in section 5, while the recurrence relations for the algebraic formulationsare already given in [48]. For instance, for equation (3.14) if P ( ν ) = (cid:88) k ≥ ρ k ( ν − ν ) k then (after the first few terms which have to be separately investigated), ρ n +4 = 1( ν − n + 3)( n + 4) (cid:32) qρ ( n ) + 8 ν qρ ( n + 1)+ (cid:0) ν q − n + a − n − q − (cid:1) ρ ( n + 2) − (2 n + 5) ( n + 3) ν ρ ( n + 3) (cid:33) . We computed this recurrence relation automatically from equation (3.14) by usingthe gfun package in Maple, specifically its diffeqtorec command [62]. Mathieugave the simpler form at ν = 0, and separated out the even and odd series so thateach recurrence relation involved only two prior terms. This means that the Mathieufunctions are what is now called D -finite or holonomic , and can therefore be computedto high precision with an asymptotically fast algorithm. See [9, 50, 51, 70].Mathieu then considered properties of the functions that could be deduced fromthese power series, which could also be interpreted as series in powers of cos α or ofsin α . In particular, he used them to count real roots. In order to solve the vibrating drumproblem, Mathieu had also to solve equation (3.8). He chose to do this in a wayanalogous 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 matchingthe fixed boundary condition at the elliptical rim of the membrane. In our terms, themodified Mathieu functions are simply the Mathieu functions with purely imaginaryargument: Ce g ( x ; q ) = ce g ( ix ; q ) and Se g ( x ; q ) = − i se g ( ix ; q ). We show two suchfunctions in figure 5.We may now discuss the details of figure 3. We chose a drum shape with aspectratio 5 : 3. We chose to look at an even mode corresponding to a ( q ), so this meansour pure tone will be described by Ce ( β ; q )ce ( α ; q ). We decided that it should be, forthe introduction, a simple picture with no elliptical nodal lines, only hyperbolic; thismeant that we were looking for the first zero of Ce ( β ; 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 (ln 2; q ) = 0. By zerofinding on q , we find that q ≈ . q . The physics of the membrane would then give the frequency of oscillation sin 2 λmt via q = λ c , with the membrane parameter m . This value of q gives the eigenvalue a ≈ . α = ± . ± π/ ± . ± [0 , , , , /
40, rememberingthat our normalization is so that ce (0; q ) = 1 = Ce (0; q ). ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE (a) Ce ( β ; 2) = ce ( iβ ; 2) (b) Se ( β ; 2) = (cid:61) (se ( iβ ; 2)) Fig. 5 . (Left) A graph of Ce ( q, β ) when q = 2 . As the argument β increases, the functionbecomes increasingly oscillatory. This value of q could be used for an elliptical drum whose verticaldimension was such as to coincide with a zero of this function (units depending on the locations ofthe foci at ± c .) (Right) A graph of Se ( β ; q ) = (cid:61) (se ( iβ ; q )) when q = 2 . ´Emile L´eonard Mathieu (15 May1835– 19 October 1890) attended the ´Ecole Polytechnique du Paris, taking the en-trance examination in 1854. He defended his doctoral thesis in pure mathematicsin 1859, before a committee consisting of Lam´e, Liouville, and Serret. He was well-regarded by the community at the time for some of his works, and indeed is stillknown today for what are called “Mathieu groups”. Apparently, though, he did notreceive the positions that he truly wanted; he then turned to applied mathematics(specifically, Mathematical Physics) to see if that would “more engage the interest ofscientific men” [26]. In spite of this change, and in spite of winning a Gold Medal in1867, he was repeatedly passed over, and in 1869 he left for a position at Besan¸con,becoming Chair of Pure Mathematics there in 1871. He remained there until 1873when he took up the Chair in Pure Mathematics at Nancy, where he remained untilhis death in 1890. His obituary is a very interesting read. In it Duhem [26] praisesMathieu’s achievements, calling him the natural successor of Poisson, and blaming“fashion” and “politics” for passing him over (thus suggesting that the fashion andpolitics 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 Duhemcontends, 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 theSorbonne, was awarded to Picard instead, who was Hermite’s son-in-law and Mathieuthought this was a scandal; to be fair, it could be argued that Mathieu’s record wassuperior at the time. But in other cases it is clearer now. Mathieu complained thathe came second to another favourites of Hermite for another Chair, in this case toHermite’s student Henri Poincar´e. The modern view must be different to Mathieu’s:It would be very difficult today to imagine choosing Mathieu over Poincar´e for anyChair.We find in [17] (which contains an interesting view of the tension between Parisand the provinces, and passages from Mathieu’s correspondence) still other reasons4 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR why Mathieu was perpetually not chosen, and it seems that the judgement of theEstablishment 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 powerhappened in 1851, and the Franco-Prussian war in which Napoleon III was capturedended when the Prussians took Paris in 1870, just as one example of how the largerworld may have intruded on academic life. It is quite believable that in these turbulentsocial circumstances many deserving people did not receive all the recognition thatthey had earned.In spite of all the difficulties of the times, however, Mathieu left a very significantbody of mathematical advances for posterity. We have surveyed only a small cornerof his work in this present paper, and even from just this it is clear today that he wasone of the best mathematicians of his age, which included some of the greatest everknown.
It was Whittaker, oneof the giants of nineteenth- and early twentieth-century mathematics, who bestowedthe name the Mathieu equation on equation (3.1) and the name the Mathieu functions on the even and odd periodic solutions of the Mathieu equation, “and these only” [74].According to Whittaker’s obituary [66], the name was given in the paper in the fifthICM [72]. Attributing a person’s name to an equation or a function is a significantevent in mathematics because people (even mathematicians) are social animals, andwe simply pay more attention when a person’s name is involved. Such namings oftenget it wrong, of course: “Stigler’s Law of Eponymy” states that no scientific discoveryis named after its original discoverer. For instance, Puiseux series are named afterVictor 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 yearsprior to that by Ricatti, but in any case was also discovered by Newton. In this casehowever we think Whittaker got it right, and Mathieu deserves the credit. The moredescriptive “elliptic cylinder functions” is also used on occasion, and was also usedby Whittaker. Nowadays Mathieu functions and Modified Mathieu functions are alsocalled Angular Mathieu functions and Radial Mathieu functions, in accordance withtheir roles in the elliptic coordinate system: β is more like a radius, and α more likean angle.A full chapter of “A Course of Modern Analysis” by Whittaker and Watson [74]is concerned entirely with Mathieu functions. Any serious study of these functionsshould begin with this book. For some reason, the authors rescaled the equationand have 16 q there instead of the 2 q (or 2 h ) that Mathieu had. This is of no realconsequence. More interestingly, they use the same normalization convention thatMathieu did in his perturbation series computations: they take the coefficient ofcos gz in the expansion of ce g ( z ) to be unity, and similarly the coefficient of sin gz inthe expansion of se g ( z ) to be unity. This differs from modern practice.Several important theorems are established in that chapter: the orthogonalityunder the inner product of equation (3.10) is proved (by appeal to results from anearlier chapter), several integral equations are established, and the Floquet theoryof the solutions to periodically-forced ODEs is applied to the non -periodic solutionsof the Mathieu equation, and more generally to Hill’s equation. The Floquet theoryshows that solutions must exist of the form exp( µz ) φ ( z ) where φ ( z ) is periodic with Achille Marie Gaston Floquet (15 December 1847, ´Epinal 7 October 1920, Nancy)ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE π (because the periodic forcing of the Mathieu equation has that period) and µ (actually, in modern works starting in [1], ν where µ = iπν ) is called the characteristicexponent . Regions where (cid:60) ( µ ) > even and therefore exp( − µz ) φ ( − z ) is also a solution;this implies that for both solutions to be stable, we must have (cid:60) ( µ ) = 0; moreoverthe regions in ( a, q ) parameter space where the solutions are stable are bounded by characteristic lines containing periodic solutions.Below is one of the integral equations established in [74]: if G ( η ) is an evenMathieu function, then there is a characteristic number λ for which (translating fromthe 16 q convention of Whittaker to the 2 q convention used in this paper)(3.16) G ( η ) = λ (cid:90) π − π e √ q cos η cos θ G ( θ ) dθ . We have not tried using this as a method of computing even Mathieu functions, al-though Whittaker and Watson claim that it “affords a simple manner of constructing”them.Whittaker and Watson attribute the change of variable ζ = cos z (the DLMFhas sin z but this is the same) which turns the Mathieu equation into an algebraicdifferential equation to Lindemann, and make further references to works of Abel,Stieltjes, Sylvester and others in the study that results. They make an asymptoticconnection of Mathieu functions to Bessel functions using this form of the equation:in this formulation it is more natural, but still a bit involved. Along the way, we see a continued fraction show up; but not the same continued fraction that we shall shortlyencounter.Whittaker and Watson also show that the Fourier series for Mathieu functionsare well-defined, at least for small enough q , by exhibiting convergent power series forthe coefficients. This marks an important step. We take the following material fromthe remarkably well-written and well-informed Wikipedia article on Whittaker, sup-plemented by the mathematical obituary written by G. Temple [66]. Sir EdmundTaylor Whittaker (24 October 1873–24 March 1956) studied maths and physics atTrinity College, Cambridge from 1892. He was elected a Fellow of the College in 1896and continued there until 1906; he was elected a Fellow of the Royal Society of Lon-don in 1905. he then became Royal Astronomer of Ireland and Andrews Professor ofAstronomy at Trinity College Dublin. In 1911 he joined the University of Edinburghand in 1912 was elected a Fellow of the Royal Society of Edinburgh, later becomingPresident. He was knighted by George VI in 1945 “for service to mathematics”.His mathematical obituary previously cited runs 22 pages, not counting the facingportrait or the four page list of Whittaker’s works at the end. It is well worth readingand conveys a substantial amount of mathematics in its own right. The topic head-ings are 1.—Algebra 2.—Interpolation (exhibiting significant work in statistics andin numerical analysis, which the anonymous Wikipedia author also takes particularnote of) 3.—Automorphic functions 4.—Astronomy 5.—Potential theory and specialfunctions [It is here that we learn that Whittaker’s 1912 paper established that theMathieu equation is “the simplest linear differential equation which is not reducibleto a particular or degenerate case of the hypergeometric differential equation”; and itis also here that we learn that E. L. Ince was a research student working with Whit-taker in Edinburgh and it was from this period that Ince gained his interest in Mathieu6
C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR functions.] 6.—Dynamics 7.—Relativity and Electromagnetic Theory 8.—QuantumTheory 9.—Scientific Books and Monographs and 10.—Historical and philosophicalwritings. The very final section, “11.—Influence” ends on a slightly ironic note: Tem-ple claims that Whittaker’s influence was felt not just by his works, but also by hiscoinage of mathematical terms, some of which are listed in the final paragraph. Sadly,of all those terms listed, namely ‘isometric circle’, ‘adelphic integral’, ‘cotabular func-tions’, ‘cardinal function’, ‘congruent hypergeometric function’, ‘Mathieu function’,and ‘calamoids’, few apart from the Mathieu functions are widely known today. Weconfess to curiosity as to what he meant by a “calamoid”: the term seems to havesurvived only in botany, and has to do with palm leaves.But the mountain of scientific achievement that Whittaker built still stands onits own.
The perturbation solutions to theMathieu equation and the series in powers of cos α or sin α suggest that it is naturalto think of using Fourier series to represent these periodic functions. But by the timeof A Course of Modern Analysis [74] it would be unthinkable not to use Fourier series.It seems, however, that it was Ince who first put them to use for Mathieu functions.The basic idea is simple: we expand our proposed periodic solution in Fourierseries. y ( x ) = (cid:88) k ≥ A k cos kx + (cid:88) k ≥ B k sin kx . (3.17)After looking at the perturbative solutions produced by Mathieu, this would havebeen a very natural idea. Insertion into the Mathieu equation (3.1) and using themultiplication identitiescos(2 x ) cos( kx ) = (cos( k + 2) x + cos( k − x )(3.18) cos(2 x ) sin( kx ) = (sin( k + 2) x + sin( k − x )(3.19)and then equating coefficients of cos kx and similarly sin kx gives us a collection of re-currence relations. By circumstance (which Mathieu made simpler by converting fromcos 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 slightlydifferent forms) produce a set of equations that can be written as follows: aA − qA = 0(3.20) − qA + ( a − A − qA = 0(3.21) − qA + ( a − A − qA = 0(3.22)and thereafter − qA k − + ( a − (2 k ) ) A k − qA k +2 = 0 . (3.23)The odd cosine coefficients must instead satisfy( a − − q ) A − qA = 0(3.24) − qA + ( a − A − qA = 0(3.25) ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE k = 1) − qA k − + ( a − (2 k + 1) ) A k − qA k +3 = 0 . (3.26)The even sine coefficients give( a − B − qB = 0(3.27) − qB + ( a − B − qB = 0(3.28)and thereafter (again that last equation already satisfies the following with k = 2) − qB k − + ( a − (2 k ) ) B k − qB k +2 = 0 . (3.29)Finally, the odd sine coefficients give( a − q ) B − qB = 0(3.30) − qB + ( a − B − qB = 0(3.31)and thereafter (again that last equation already satisfies the following with k = 1) − qB k − + ( a − (2 k + 1) ) B k +1 − qB k +3 = 0 . (3.32)All four of these sets of recurrence relations give rise to infinite tridiagonal matrices,all but the first one symmetric. By an artifice, namely replacing A by A √
2, we maysymmetrize even the first one.The infinite eigenvalue problem then becomes(3.33) √ q · · ·√ q q · · · q q · · · q q · · · q
64 . . .... ... ... ... . . . . . . √ A A A A A ... = a √ A A A A A ... The eigenvalues of this matrix are denoted a ( q ), a ( q ), a ( q ), . . . and indeed for real q these occur in increasing order: a ( q ) < a ( q ) < a ( q ) < · · · . For complex q it is morecomplicated, but the set of eigenvalues remains discrete and countable, essentiallybecause the diagonal of the matrix above contains entries that grow quickly enough.The other three tridiagonal matrices are constructed analogously and have anal-ogous properties. For completeness, we include them.(3.34) q q · · · q q · · · q q · · · q
49 . . .... ... ... . . . . . . A A A A ... = a A A A A ... C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR
The eigenvalues of equation (3.34) are denoted a k +1 ( q ).(3.35) q · · · q q · · · q q . . .0 0 q
64 . . .... ... ... . . . . . . B B B B ... = a B B B B ... The eigenvalues of equation (3.35) are denoted b k ( q ).(3.36) − q q · · · q q · · · q q · · · q
49 . . .... ... ... . . . . . . B B B B ... = a B B B B ... The eigenvalues of equation (3.36) are denoted b k +1 ( q ). Remark . If q is real, then these are real symmetric eigenvalue problems for which fastalgorithms are available. If q is complex, then the matrices are complex symmetric,not Hermitian. This has some important consequences that we will pursue in section 4.It may be simplest nowadays, given the modern familiarity with matrices, to thinkof these infinite complex symmetric tridiagonal matrix eigenproblems as determiningthe eigenvalues and corresponding eigenfunctions of the Mathieu equation, on trun-cation. This, or at least in terms of determinants, seems to have been the way thatInce thought of the problem [44], although he did not have direct numerical methodsto 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 andfor the coefficients A k and B k of the corresponding eigenfunctions. These series andtheir computation remain of interest, with (to quote Blanch [14]) “an algorithm forgenerating the successive [series coefficients], suitable for computers” being publishedin [61] . Subsequent papers give similar recurrence relations, and, if the order islarge enough and not too many terms are needed, explicit symbolic formulae. Such aformula is termed “generic” in [32]. This last paper is interesting: the author, a researcher at “TRG Incorporated in Melville, NY”appeared to be disappointed at their main result, which was a set of recurrence relations whichcould be used numerically to generate the numerical series coefficients of both eigenfunctions andeigenvalues: it seemed that they were really searching for a good test problem for early symboliccomputation programs! We were intrigued, but unable to find out much about TRG Incorporatedexcept that staff there carried out research on computational fluid simulations and lasers in the1960’s. As a further remark on this quotation, we are not sure if Blanch is referring to human computers or to machines, although we incline to the latter interpretation because she does use twokinds of IBM computers for that article.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE Fig. 6 . The first few eigenvalues a g ( q ) and b g ( q ) , that is, the eigenvalues of the infinite tridiag-onal matrices in equations (3.33) – (3.36) , as functions of the real variable q . This graph is modelledafter the similar one in Fig 28.1 of the DLMF. Ince also derived a continued fraction from the above recurrence relations. Wewill continue to use the eigenvalue idea for the moment, as is done in many numericalmethods today such as [76]. The paper [19] makes especially good use of this, andadvocates strongly for its conceptual and computational use.When q = 0 these eigenvalues are simply the squares of the whole numbers: 0,1, 4, 9, and so on. They are, technically, double eigenvalues (taking even and oddfunctions together), but in this case they retain two independent eigenfunctions cos kx and sin kx , except if k = 0. For real q the eigenvalues can be sorted in increasing order.See Figure 6.The evenness of the Mathieu equation, and its invariance under the transformation z → z ± π/ q → − q , mean that the eigenvalues have the following symmetries: a n ( − q ) = a n ( q ), a n +1 ( q ) = b n +1 ( q ), and b n +2 ( − q ) = b n +2 ( q ). These are equa-tions 28.2.26–28.2.28 in the DLMF. There is also the conjugate symmetry a ( q ) = a ( q )and similarly for b ( q ).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/eigenvectorformulation 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 The authors of [19] attribute the derivation of the continued fraction form to Heine in [39] butit seems likely that Ince’s derivation was independent. C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR(a) Convergence when q = 2 (b) Convergence when q = 15 + 4 i Fig. 7 . What dimension matrix is needed to get an accurate value of eigenvalue a k ? Here weplot the errors in a k ( q ) for ≤ k ≤ by using a matrix of dimension k + ∆ n where ∆ n is asindicated on the figure: either ≤ ∆ n ≤ as on the left, or twice that as on the right. That is,for a (2) (top symbol in each column of points in the left hand graph) we plot the errors in usingmatrices of dimension , , , , , and . This tells us that to get double precision weneed dimension (dimension almost works). In contrast, for a (15 + 4 i ) , (top symbol in eachcolumn of points in the right-hand graph) we consider matrices of dimension , , . . . , .In this case we need dimension to get double precision accuracy in this eigenvalue. Higher-ordereigenvalues 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 smallerones do. When k = 30 (lowest points in each column) the accuracy is significantly better; the matrixdimension for the low point on the right is . All computations done in Digits. software package, one has to pay attention to the choices that the authors have made.Mathieu’s original normalization consists in setting the coefficient corresponding tocos gx to 1 (or that of sin gx in the odd case), and this is convenient for perturbationcomputations here as well, when it is possible.The eigenvalues of truncations of these infinite matrices converge quite quicklyto the desired eigenvalues of the Mathieu equation although the rate of convergencedepends on the eigenvalue. See figure 7. There are convergent series containingeigenvalues that can be used as numerical checks on the accuracy, but for this graphwe simply computed high-precision values of the eigenvalues by a continued fractionmethod and compared the results to the truncated matrix eigenvalues. As stated, thematrices 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 [58].If q is complex, then the matrices are complex symmetric tridiagonal, which are harderto solve, but for which there are still interesting algorithms [5]. There are specializedalgorithms to find multiple eigenvalues in the infinite dimension case [52]. We have notfelt the need to resort to fast methods: an ordinary slow ( O ( n )) eigenvalue algorithmis perfectly fine, because the matrix dimensions are so small in modern terms.After all, the Fourier series converges spectrally to the Mathieu equation [64], andthus relatively few eigenvalues are needed. At least, for the real case. There is some asymptotic work for Mathieu functionsin [74], and indeed even a little in Mathieu’s original work; but the first serious studies
ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE h = √ q and s = 2 m + 1, then as h → ∞ with m = 0, 1, 2, . . . , a m (cid:0) h (cid:1) b m +1 (cid:0) h (cid:1)(cid:27) ∼ − h + 2 sh −
18 ( s + 1) − h ( s + 3 s ) − h (5 s + 34 s + 9) − · · · . 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 algebrasystems, for instance Maple. In [42], Inceproved that the Mathieu equation (and similarly the Hill equation) “can admit butone solution of period π or 2 π ”. This is generally taken to mean that the otherlinearly-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 [43]Ince established for small q that if, for instance, the periodic solution was ce p ( z ), thenthe second, necessarily non-periodic, solution would be y = se p ( z ) + Kz ce p ( z ) + φ ( z )where both K and φ ( z ) would be O ( q p ) in size and φ ( z ) would be periodic. As wecan see in figure 8(a) this is true even if q is not small (in fact we took q = 1 . i approximately, and its eigenvalue a ≈ . kπ , where k is an arbitrary integer greater than 2. For these eigenvalues, all solutionsare periodic, if one is.” She gives the reference [57] in Chapter 20 of [1]. This is alsocovered in [28]. Edward Lindsay Ince (30 November 189116 March 1941) was, as previously mentioned, a research student of Whittaker andhad spent the years 1909–1913 at the University of Edinburgh. He had a short butcolourful career, which included what would have then been an exotic episode atthe newly-founded Egyptian University of Cairo from 1926 to 1931. Returningto England for reasons of health and family, he took various positions, ultimatelybecoming Head of Technical Mathematics at the University of Edinburgh. He waselected a Fellow of the Royal Society of Edinburgh in 1923. He died of leukemiain Edinburgh, aged just 49. Just before his death he was awarded the Makdougall-Brisbane prize for his work on Lam´e functions, but he died before he learned of theaward.There are two obituaries of Ince that we know of: one by Whittaker [73], andanother by Aitken [3]. Both make interesting reading, although there is more math-ematics in the one by Whittaker. The one by Aitken is perhaps more personal,although it does contain the passage “Ince firmly believed that theoretical solutionsof problems, however abstractly elegant, were incomplete unless the mathematician Founded 1908, it is now Cairo University. C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR(a) A non-periodic solution (b) Unstable solutions
Fig. 8 . Some non-periodic solutions of the Mathieu equation. On the left, we have the non-periodic solution for the Mulholland-Goldstein double point q ≈ . i showing secular growth asInce proved. Black is the real part of the solution, red is the imaginary part. On the right weshow absolute values of solutions for a = 4 and q = 10 , well into an unstable region in ( a, q ) spaceby the Floquet theory, showing exponential growth. Indeed the solutions are of the predicted form exp( µz ) φ ( z ) where Re ( µ ) > and φ ( z ) is periodic; we can see that this particular periodic functionhas two zeros by the “spikes” in the graph. either tabulated the solving functions himself or rendered them tabulable. Perusalof his papers will show that in his chosen field of research he achieved both of theseobjects”. As an aside, Aitken’s biography is itself worth reading; in addition to hismathematical accomplishments he was elected to the Royal Society of Literature forhis 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 theFirst 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 outreading mathematics under Professor George Chrystal, who we principally know nowa-days as the person who wrote the text that inspired Ramanujan; after Chrystal diedin 1911, “Ince continued under his successor”, namely Whittaker himself. Whittakerthen carefully describes Ince’s achievements with the Mathieu functions, includingInce’s proof [42] in 1922 that there could not be two independent solutions with pe-riod π or 2 π to the same Mathieu equation. Whittaker echoes Aitken’s remark aboutInce’s sensibility with regard to computation, and amplifies it: “Ince held that animportant part of a pure mathematician’s duty is to provide tables for the use ofphysicists and astronomers, and he was well aware that the possibility of construct-ing such tables without a colossal expenditure of time and energy depends on theprogress of theoretical analysis.” Whittaker remarks that Ince’s 1932 tables of theMathieu functions their zeros and turning points was “A splendid piece of work, per-formed single-handed save for some help by an Egyptian assistant”, without givingthe name of the assistant. Ince himself was more gracious, and acknowledged MansyShehata, who was then an Assistant in Pure Mathematics at the Egyptian Universityin Cairo. Ince also acknowledged grant support in purchasing calculating machines,which seemed to be of significant use. ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE
Blanch’s first pub-lication on Mathieu functions was [11]. This was reprinted in 1950 by the NationalBureau of Standards “with the permission of the editors of the Journal of Mathematicsand Physics to meet a continuing demand”. This refined and improved the contin-ued fraction method that Ince used and in particular allowed the error in the Fouriercoefficients to be estimated. She wrote a paper in the Transactions of the AmericanMathematical Society on the asymptotics of the odd periodic Mathieu functions [12],extending and correcting the works of others. Blanch wrote an influential paper inSIAM Review on continued fractions, namely [13]. This paper contained a detailed lin-earized rounding error analysis for continued fractions, and used both Bessel functionsand eigenvalues of Mathieu functions as examples. She argued that continued frac-tions gave the preferred method for computing Mathieu function eigenvalues in [14].Then Blanch & Clemm [15], still later improved by [40], went on to systematicallycompute the double eigenvalues .We believe that Blanch’s choice of the method of continued fractions was at leastpartly because of the kind of computing she typically did, which in her own wordswas “experimental in nature” and which by habit were continually checked as thecomputation went on. [Our computations for this paper were carried out in a similarstyle.] In the first part of her career the computations were, like Ince’s, carried outby hand and by desk calculating machine. A significant difference is that while Incehad only one assistant when he was in Cairo, Blanch organized a group of as many as450 assistants, and later used digital computers. But her computations were in manyways what we might call “artesenal” with a great deal of human involvement.Today people generally prefer software packages that can run without care onthe part of the user: in Blanch’s words, in a “robot-like” manner [13]. Most peoplewant 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 writinggeneral-purpose software for evaluation of the Mathieu functions, therefore, one might not choose continued fractions (as we will see) because their occasional instabilities inthe complex plane are somewhat variable; on the real line Blanch had analyzed thisbehaviour 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 becontrolled simply by using extra precision.Blanch’s variant of the continued fraction algorithm can be described as follows.We suppose first that a and q are given (later we will perform iterations, looking forvalues of a and sometimes also of q that make the continued fraction equal to zero).We will need the auxiliary quantities(3.37) V m = a − m q for m = 0, 1, 2, . . . . These are not defined if q = 0, but if q = 0 then we already know Blanch 1966 quotes using an “I.B.M. 1620 and 7094”. We merely note the printed periods inthe abbreviation “I.B.M.” used there, as opposed to the simpler trademark IBM used nowadays, andrefrain from comparing the speed and memory of the early computers with what is available today(or what is coming). C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR that the eigenvalues are g for integers g . Equation (3.37) is equation (1.05) in [14].We use this notation in the recurrence relations (3.20)–(3.32), and further with theauxiliary quantities c m where c = 2 and all other c k = 1, in order to look after thefirst edge case. We rewrite the recurrence relations (apart from the base cases) as(3.38) A m +2 + c m − A m − − V m A m = 0 . Putting G m = A m /A m − and dividing equation (3.38) by A m (Blanch was very carefulabout what happened when any A m was zero, but here we rely on IEEE arithmeticwith signed zeros and infinities to get everything right) we get(3.39) G m +2 + c m − G m − V m = 0 . This recurrence can be written either as(3.40) G m +2 = V m − c m − G m or(3.41) G m = c m − V m − G m +2 Of course, these recurrences must be started correctly by using the appropriate equa-tions (3.20)–(3.32). Running equation (3.41) until the denominator Q is so large that1 / | Q | ≤ ε where ε is our tolerance, starting from some M > c k = 1 for all k ≥ M − G M = 1 V M − V M +2 − VM +4 − ...and because the V m s grow like m this continued fraction converges for all a and all q (cid:54) = 0. Call the G M computed in this way G M, tail .Now we use equation (3.40) with increasing m starting from our known edge cases(depending on which class of Mathieu eigenvalue we wish to compute), and if and onlyif a is an eigenvalue then the two values of G M in the middle will agree. Call the G M computed in this way G M, head . Let(3.43) T ( a, q ) = G M, head ( a, q ) − G M, tail ( a, q ) .T ( a, q ) must be zero for a to be an eigenvalue corresponding to q . The edge casesfor G M, head determine whether this is an a k , a k +1 , b k , or b k +1 eigenvalue. Just which integer k depends, as a rule, on whether there is an unambiguous continuouspath 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 numericalinstability for real q .There are many methods one could use to find zeros of equation (3.43), butsince differentiation of T ( a, q ) with respect to a is simply a matter of differentiatingequations (3.40) and (3.41) (convergence is uniform in compact neighbourhoods andso differentiation is permissible), we choose Newton’s method. For a given q andstarting with an initial estimate a (0) for the eigenvalue , the iteration is a ( k +1) = a ( k ) − T ( a ( k ) , q ) T a ( a ( k ) , q ) . In order to find a good initial estimate, people usually use continuation from q near 0: one solvesfor q = q n , and then uses that eigenvalue as an initial estimate for the eigenvalue at q = q n + ∆ q .ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE q as initial estimates for theeigenvalues she required, and this worked very well, unless, of course, the eigenvaluewas not simple ( i.e. multiplicity 1), in which case even more derivatives of T turn outto be useful.The eigenvalues of the Mathieu equation are as previously stated almost alwayssimple, but for isolated complex values of q may have multiplicity 2. In particular,if q = is where s is real and i = −
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 . s = 6 . s = 1 . q for which two (or more) pairs of eigenvalues merge.One might be tempted to dismiss double points because they occur “with proba-bility zero” if one chooses the parameter q “at random”. But in applications requiringcomplex q the parameter will not usually be chosen at random, and indeed the prob-lem being modelled may call for a continuum of values of the parameter. In theapplication 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 somesense this means that the “probability zero” events actually occur with “probabilityone”. This reversal of expectation is a common occurrence in bifurcation phenomena,for instance. Blanch and Clemm used a variation on Newton’s methodin their systematic search for double points [15]. Their method was a not elegant,compared to two-dimensional Newton iteration, but it worked and it was the firstmethod to do so. Instead of just using one derivative with respect to a , they usedtwo, and expanding T ( a, q ) = T ( a ( k ) , q ) + T a ( a ( k ) , q )( a − a ( k ) ) + 12 T aa ( a ( k ) , q )( a − a ( k ) ) + · · · they set this to zero and solved the resulting quadratic for the update to a ( k ) , choosingthe smaller magnitude square root. This gives a superlinearly-converging iteration fordouble roots (even faster for simple roots), and together with interpolation enabledthem to compute the smallest 72 double eigenvalues to what we would now call doubleprecision. We confirmed their work and recomputed all their roots. We used theirtabulated values as initial estimates, and have plotted the results in figure 9. Blanchand Clemm carried out their computations on digital computers (we believe modelsIBM 1620 and IBM 7024). The works of Gertrude Blanch werepossibly not as familiar to the reader as were those of Mathieu, Whittaker, or Ince. Ifthat were true, we hope that this paper has helped to improve matters. As previouslystated, Gertrude Blanch wrote the chapter on Mathieu functions in the classic hand-book [1]. She apparently wrote in 1943 what might be the first modern textbook onnumerical analysis, and an updated version in 1982, according to [36]. Unfortunately,we have not seen a copy of either edition. We have cited above many of her paperson the Mathieu functions and related matters.We draw material for this section from several sources, including a transcriptof an interview with her in 1973 [67], an extensive biography by Grier [36], and6
C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR the collection of her papers at the Charles Babbage Institute at the University ofMinnesota. Gertrude Blanch’s extremely interesting life was well-documented. Wehave in this short biography concentrated on mathematical aspects of her life andleft out important non-mathematical aspects. The above resources are well worthconsulting 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 in1914. She changed her name, after her father died, to an Anglicized version of hermother’s family name, Blanc. She became an American citizen in 1921. She workedfor fourteen years to get enough money to attend university; her employer paid hertuition for New York University, where she graduated summa cum laude in 1932.She then went to Cornell, receiving her PhD in algebraic geometry in 1935 underthe guidance of Virgil Snyder and Wallie A. Hurwitz (she is not included in theMathematics Genealogy Project at the time we are writing this review; this accountis from her interview with Tropp and from the Cornell web page). After a short stintteaching at Hunter College for someone on sabbatical leave, she took an office clericaladministration and accounting job. This administrative job was to prove importantfor her later work with the Mathematical Tables Project, as detailed in [36]. Inorder to keep her mathematical interests alive, she took a night course in relativity atWashington Square College given by Arnold Lowan. When Lowan was asked to createthe Mathematical Tables Project under the New Deal Works Progress Administration,he asked Blanch to join. She became Technical Director, eventually organizing a groupof 450 (human) computers. According to Grier, she deserves much of the credit forthe success of the project, and part of that credit is due to her prior business-orientedadministrative experience. She published several papers during this time, includingone with Hans Bethe [16].“During her time at the Mathematical Tables Project, she particularly enjoyedworking with the Mathieu functions, and these functions would be central to the restof her career.” [36, p.23].When asked how she first got interested in Mathieu functions, she responded:“Morse was interested, for example, in Mathieu functions. I got started on Mathieufunctions because Morse needed them and there were any number of small things andsome special integrals that he came across within his field.”She remained interested in Mathieu functions for the rest of her career and indeedafter her retirement. Her ultimate academic appointment was as Head of Mathemat-ical Research in the Aerospace Research Laboratory at Wright Paterson Air ForceBase in Dayton, Ohio, where she wrote many of her papers. She became a Fellowof the American Association for Advancement in Science in 1962. She received theFederal Women’s Award from President Lyndon Johnson in 1964. She retired in 1967,and died in 1997, aged 99. There is an excellentconcise treatment of Mathieu functions, including the Floquet theory, in ChapterXVI of Volume 3 of [28], “The Bateman Manuscript”. In particular, it is there thatwe learn two interesting facts about µ = iπν , the Floquet characteristic exponent forthe Mathieu equation: first, that Poincar´e found a way to compute it from the twobasic solutions w I ( z ) and w II ( z ) using the evenness of the Mathieu equation: at least Philip Morse, at MIT, the first author of [53]. This monumental work also describes Mathieufunctions.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE w I ( z ) = e µz φ ( z ) + e − µz φ ( − z )2 φ (0)or w II ( z ) = e µz φ ( z ) − e − µz φ ( − z )2( φ (cid:48) (0) + µφ (0))will have nonzero denominator; now differentiate w II and compute w I ( π ) and w (cid:48) II ( π ).Since φ (0) = φ ( ± π ) and φ (cid:48) (0) = φ (cid:48) ( ± π ), we have that(3.44) cosh πµ = w I ( π ) = w (cid:48) II ( π ) . In the DLMF this equation (using cosine and not hyperbolic cosine, because ν isused instead, where µ = iν ) is called the characteristic equation , number 28.2.16.Blanch uses another convention, namely µ = iπν , defining ν differently, in [1]. Thesecond 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 solutionis periodic with period π and if n is odd the solution is periodic with period 2 π . Thecase when µ = i r where r is a rational integer generates those mysterious solutionsalluded to earlier which are periodic with greater periods. See section 20.3 in [1] formore discussion of this fact.The earlier book by M.J.O. Strutt (1932) was critically reviewed in [35], whostated “It is a useful book, both for pure mathematicians interested in the theoryof special functions, and for applied mathematicians compelled to use the functionsin their researches. It is worthy of consultation by both classes, but it is rathersuperficial . . . ”.The most extensive and thorough treatment of Mathieu functions is that of [49],which is based on an earlier book by Meixner and Sch¨afke that we have been unableto get hold of. The later edition handles significant generalizations of the theory ofMathieu functions, and as previously stated deals completely with the theory of doublepoints. Nearly every reference antedating this work cites it; indeed the BatemanManuscript 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 workremains the most extensive, but we can also mention [10]. Tables of integrals andseries for Mathieu functions are in [59]. The value of the numerical tables nowadaysis of course lessened by the availability of software; the tables of integrals and seriesmay also have been supplanted by computer algebra systems, at least to some extent.There are a large number of relatively modern books, papers, and software pack-ages for the computation of the Mathieu functions. Analytical work on the Mathieufunction is still extensive; to cite a paper almost at random, consider [55]. Some ofthe works ignore most of the others, and hence rediscover things, but others are veryinsightful. We particularly recommend [19] and [38] for exceptional visualizationsand clarity of exposition. On the purely computational side, there are three sets ofACM Transactions on Mathematical Software papers discussing implementations, thelatest being [29]. There are Python (scipy) and third-party Matlab implementationsof (real) Mathieu functions. The Mathematica implementation of complex Mathieufunctions may be very good (unfortunately, we have only had limited access to Math-ematica so we cannot be sure, but it may even have facilities for computing doubleeigenvalues) and is described in [68]. The Maple implementation, with which we are8 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR most familiar, encodes a substantial amount of analytical work including both q -seriesand asymptotic series.
4. Double points and Frames.
As previously noted, at isolated points inthe complex q -plane, for instance at the Mulholland-Goldstein double point q ∗ ≈ . i (reporting 16 digits ), we have a double eigenvalue: a = a = a ∗ ≈ . ( q ∗ , α ) = ce ( q ∗ , α ), leaving a gap in the com-pleteness of the set of eigenfunctions. This means that we will have to supplementthe set of eigenfunctions in some way in order to expand arbitrary functions in seriescontaining Mathieu functions.Second, near to these double points, the ordering of the eigenvalues becomesambiguous. For (cid:61) ( q ) < (cid:61) ( q ∗ ) we have a ( q ) < a ( q ) < a ( q ) < · · · , but at q ∗ equalityoccurs. For (cid:61) ( q ) > (cid:61) ( q ∗ ), both a ( q ) and a ( q ) are complex, and ordering is a matterof convention. The DLMF adopts the convention in this case that a ( q ) continuesas (cid:61) ( q ) increases by choosing the branch with negative imaginary part, while a ( q )takes the conjugate. See the visualization in section 28.7 of the DLMF.This convention makes some sense provided one thinks of q varying along theimaginary axis, and increasing. Approaching q ∗ along some other path in the q -planemay cause puzzlement: the ordering is conventional, no more.Also as q approaches q ∗ , the coefficients in the Mathieu series expansion for agiven function must become singular. For example, consider(4.1) cos 2 α = c ( q )ce ( q, α ) + c ( q )ce ( q, α ) + · · · . By numerical experimentation at high precision we find that c + c = 1 . − . i (4.2) c − c = − . . i √ q − q ∗ (4.3)when we use the normalization convention that ce k ( q,
0) = 1, and moreover that a = a ∗ + d · ( q − q ∗ ) / + O ( q − q ∗ )(4.4) a = a ∗ − d · ( q − q ∗ ) / + O ( q − q ∗ )(4.5)where d ≈ . . i = 1 . i ).These values were found by using orthogonality, which holds if q (cid:54) = q ∗ , and by high-precision computation of (cid:90) πx =0 ce k ( q, x ) dx = O ( q − q ∗ ) / . Then if y ∗ = ce , ( q ∗ , α ) is the coalesced solution to Mathieu’s equation, we can exam-ine the most important contribution to the perturbation by considering the solution This double point was the first found: studied in [54] and later computed by [18] to 3 digits andthen to double precision in [15]. We are slightly embarrassed to admit to how many figures we took these computations to: weworked at 250 decimal Digits in Maple, and solved the Mathieu equation with a tolerance of 10 − .We then worked with 5, 8, 13, 21, 34, and 55 Digit truncations of q ∗ and calculated the corresponding a and a to 120 Digit accuracy; this enabled us to identify the constants in this section to 50 Digitsor more. We only report the double precision values. Later, we confirmed these by simply computingthe Puiseux expansion of this point.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE y (cid:48)(cid:48) + ( a ∗ ± d · ( q − q ∗ ) / − q cos(2 α )) y = 0 . If y = y ∗ + u ( q ∗ , α ) · d · ( q − q ∗ ) / + · · · , then a short calculation dropping terms of O ( q − q ∗ ) and higher gives(4.7) u (cid:48)(cid:48) + ( a ∗ − q ∗ cos 2 α ) u + y ∗ = 0 . The function u must be periodic with the same period as y ∗ . It can therefore be com-puted numerically alongside y ∗ (by solving a boundary-value problem) or alternativelycan be expressed as an integral of Mathieu functions against y ∗ . This argument isextended and formalized in a short section in [49] starting on p. 82. Here and with justthis simple example, we see that u is essentially ∂y/∂a (found by solving a variationalequation). By combining the equations above, we can see that(4.8) cos 2 α = ( c + c )ce ∗ , ( q ∗ , α ) + d ( c − c ) √ q − q ∗ · u ( q ∗ , α ) + O ( q − q ∗ ) / . This arrangement is continuous as q → q ∗ because ( c − c ) √ q − q ∗ is O (1) in thatlimit. This shows explicitly that by adding ∂ ce , ( q ∗ , α ) /∂a to the collection of Math-ieu functions we are expanding with we ensure completeness of the set and the possi-bility of the expansion. We know of no freely-available software package for Mathieufunctions that provides for the computation of these extra functions. We discussmethods to compute u ( α ) in section 5.4. Here is a simple method for computingthe double eigenvalues (we found out afterwards that this method was also usedby [40]). When computing T ( a, q ), compute also the derivatives T a ( a, q ), T q ( a, q ), T aa ( a, q ) and T aq ( a, q ). Then a two-dimensional Newton iteration looks like this:(4.9) (cid:20) T a T q T aa T aq (cid:21) (cid:20) ∆ a ∆ q (cid:21) = (cid:20) − T − T a (cid:21) where all function evaluations and derivative evaluations occur at the current es-timates ( a ( k ) , q ( k ) ). Given sufficiently good initial guesses, this iteration convergesquadratically to double points ( a ∗ , q ∗ ). Remark q ∗ to double precision, one cannot naively compute a ∗ to the same precision: after all, it is sensitive to errors on the order √ q − q ∗ . Atdouble precision, this means one has a ∗ to about single precision, and it will havesplit into two simple eigenvalues. Their mean , of course, will converge well to thetrue double eigenvalue. This numerical split, however, might be used to advantage incomputation of independent eigenfunctions: the norms of these erroneous approximateeigenfunctions will be O ( q − q ∗ ) / and therefore the associated Fourier coefficients inthe eigenfunctions will be large, but perhaps this is tolerable. The resulting spectralexpansion of the function may well be accurate enough for one’s purposes. Perhapsthis is the real reason no-one has developed the tools to deal precisely with theseextra eigenfunctions in practice; of course in theory this has been understood since atleast [49].A significant advantage of Blanch’s version of the continued fraction algorithm isthat it can be carried out in series . One simply uses series arithmetic when adding,multiplying, or dividing. This automatically allows the computation of all derivatives0 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR needed. More, this allows computation of both local Taylor series for the eigenvalues,that is a ( q ) = (cid:88) k ≥ α k ( q − q ) k , by carrying out the Newton iteration in series [34] with q = q + x where x is theseries variable, and also Puiseux series a ( q ) = a ∗ + (cid:88) k ≥ α k ( q − q ∗ ) k/ for the eigenvalue about double points, by carrying out Newton iteration in series with q = q + x where x is the series variable and with the initial estimate a ( q ) = a ∗ + α x where α is one of the two roots of(4.10) 0 = T ( a ∗ , q ∗ )+ T a ( a ∗ , q ∗ )( α x + · · · )+ T q ( a ∗ , q ∗ ) x + 12 T a,a ( a ∗ , q ∗ )( α x ) + · · · . That is,(4.11) α = ± (cid:18) − T q T a,a (cid:19) / . Iterating in this manner will give eventual quadratic convergence—nearly doublingthe number of terms correct with each iteration—in computation of the Taylor orPuiseux series.
For the doubleeigenvalue a ∗ = 2 . . . . corresponding to the Mulholland-Goldsteindouble point q = q ∗ = 1 . . . . i , we have(4.12) a = a ∗ + α √ q − q ∗ + α ( q − q ∗ ) + α ( q − q ∗ ) / + · · · . Computation according to the method of the previous section gives that α ≈ ± . . . . (1 + i ) α ≈ − . iα ≈ α · ( − . i ) α ≈ − . α ≈ α · (0 . α ≈ − . iα ≈ α · (0 . i ) α ≈ . α ≈ α · (0 . . (4.13)Puiseux series can be computed about every double point, but we do not believe thatthese have been reported in the literature. We see a certain amount of unexplainedregularity in this series. In particular note that this particular α is a multiple of(1 + i ), which is a consequence of the purely imaginary character of this first doublepoint because √ i = (1+ i ) / √
2. In table 1 we tabulate the first few coefficients of somerepresentative series tabulated, essentially as an homage to all the great tabulators ofMathieu functions. Rounding errors will make the leading terms not exactly zero. It is important to trim them offand iterate only with the terms that should not be zero.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE Table 1
Some Puiseux series coefficients about double points. Coefficients computed in digits ofprecision and verified at that precision; rounded to five places and cut-and-pasted into this table.The largest double point off i R computed by Blanch and Clemm is also included. q ∗ a ∗ α α α α . i . ± . i ) − . i − . i α − . . i . ± . i ) − . i − . i α − . . i . ± . i ) − . i − . i α − . . . i ) (324 . . i ) ± (3 . . i ) (1 . − . i ) (0 . − . i ) α − (0 . . i ) Consider figure 9 where the smallestseventy-two double points , are plotted: as we said, there are a lot of them. Mor-ever, as we saw above, for values of q close to double points, the vanishing of the normof the affected Mathieu functions means that the expansion coefficients must becomesingular, and thus numerically troublesome: at the very least, there will be cancel-lation error entailed by the subtraction of large nearly equal quantities. We knowof 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 aselementary linear algebra: consider the analogous problem of finding the eigenvectorsof the following matrix, and expanding another vector as a linear combination thereof: A = (cid:20) a t a (cid:21) which if t (cid:54) = 0 has eigenvalues a ± √ t , and linearly independent eigenvectors v =[1 , √ t ] T and v = [1 , −√ t ] T . Expanding (say) [1 , T = c v + c v requires c + c = 1and c √ t − c √ t = 1 or c − c = 1 / √ t . This is obviously analogous to the situationabove. It is even more analogous when one considers the generalized eigenvector thatarises at t = 0: Au = a u + [1 , T . Exactly as in the Mathieu function case above,the numerical difficulties in expanding as a linear combination of eigenvectors showup for small nonzero t , but these are alleviated on adding the generalized eigenvectorto the mix, and writing instead(4.14) (cid:20) (cid:21) = c (cid:20) √ t (cid:21) + c (cid:20) −√ t (cid:21) + c (cid:20) (cid:21) Now, of course, the set is not linearly independent, and one has to choose the coeffi-cients in a sensible way. This leads to the modern theory of frames . See [2].
5. Algorithms for Mathieu functions and Modified Mathieu functions.
Once one has the eigenvalue, one needs a way to construct the eigenfunction. Hereare several methods. We took the tabulated values in [15], used Optical Character Recognition to convert them tocomputer-readable form, and ran our algorithm (which is the same as that of [40] and different fromthat of [15]) from section 4 to confirm them. After correcting several amusing OCR errors includingthe near-inexplicable occurrence of Russian characters masquerading as numerals, we found (asdid [40]) that all printed decimals in the tables of [15] were correct. Blanch and Clemm did notplot their double eigenvalues but only tabulated them, but Hunter and Guerrieri did plot some oftheirs. More, they computed farther into the complex plane than did Blanch and Clemm, and foundasymptotic formulae. In the figure, we see arcs of apparent square-root like curves spreading topositive infinity; we also see oval arcs of similar beads coming from the imaginary axis to the realaxis, forming a kind of peacock’s tail of double points. C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR(a) Blanch and Clemm double points (b) double points in all quadrants
Fig. 9 . In [15] some forty double eigenvalues a k and thirty-two double eigenvalues b k weretabulated. We plot their results here in figure 9(a). In the first quadrant, eigenvalues a k mergingwith eigenvalues a k +2 are plotted as blue diamonds. Eigenvalues a k +1 merging with eigenvalues a k − are plotted as black circles. Eigenvalues b k merging with b k +2 are plotted as red diamonds.Eigenvalues b k +1 merging with b k − are plotted as fuchsia circles. We see arcs of apparent square-root like curves spreading to positive infinity; we also see oval arcs of similar beads coming from theimaginary axis to the real axis, forming a kind of peacock’s tail of double points. Extension to all fourquadrants in figure 9(b) follows by conjugate symmetry and by the symmetries a k ( − q ) = a k ( q ) , b k ( − q ) = b k ( q ) , and a k +1 ( − q ) = b k +1 ( q ) , which last implies that the fuchsia circles and blackcircles exchange meaning in the left half plane: merging a k +1 and their conjugates are fuchsiacircles in the left half plane, black circles in the right half plane; merging b k +1 and their conjugatesare 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 , and they do so while keeping their independenteigenfunctions, which become cos kα and sin kα at q = 0 . The most straightforward method for com-puting an entire function is to use its Taylor series at a convenient point, say theorigin. 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 canget an accurate answer. For instance, the basic Mathieu function w I ( a, q, z ), called MathieuC in Maple, has the Taylor series beginning(5.1) w I ( a, q, z ) = 1 + a − q z + ( a − q ) − q z + · · · and since the function is entire, this series converges for all z . However, the series isimpractical, as we demonstrate now.Taking (say) the modest values a = 29 /
20 and q = 3 /
5, we still use several cpuminutes (on a Surface Pro) to compute all the terms in the above series up to O ( z ).This expense might be considered as pre-computation cost, and ignored, and so wedo. This truncated series can then be used to compute, for instance(5.2) w I ( a, q, . i ) = − . , using 46 Digits of precision in the computation; the final nonzero term of the truncatedseries is about − . · − and because all subsequent terms are smaller, and alternate,we would hope that the computed value is accurate. But the final 12 Digits, coloured ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE in sizeand 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 educationalpapers; see e.g. [23], which discusses it as an effect of the ill-conditioning of thetruncated Taylor series.The situation gets very much worse as | z | increases. Already by z = 4 . i thebest accuracy we can achieve is double precision (by using 19 + 16 = 35 Digits); by z = 4 . i only single precision by using 22 + 8 = 30 Digits; and by 4 . i only halfprecision, by using 23 + 4 = 27 Digits. By z = 5 . i it is already true that 800 termsare not enough in the series, no matter what precision we use because the truncationerror is too big. Moreover, even if we had enough terms, because the largest terms areat least 10 in size, the precision would need about 40 extra decimal digits anyway.The final conclusion is that this method is unaffordable. It can be rescued byusing analytic continuation instead of a single series (which of course is the underlyingbasis for most numerical methods to integrate ODE), and we look at that method insection 5.2. Most researchers instead choose to jump straight to a spectral method forthe periodic Mathieu functions, expanding in Fourier series; and then by a trick thesecan be used also for the modified Mathieu functions. We will discuss this method insection 5.3. One straightforward thing to try,valid for both the Mathieu functions and for the modified Mathieu functions, is nu-merical solution of the Mathieu equation. We could in principle use any method thatallows integration along a complex path, such as any Runge-Kutta method or mul-tistep method. Because these solutions are to be later used as functions, requiringevaluation at any point in their domain, we will need good interpolants and moderncodes automatically supply these.However, because the equation is linear, and because its variable coefficient has aknown Taylor series algorithm, it is a straightforward exercise to develop a specializednumerical solution by marching Taylor series, a method known as analytic continua-tion if infinite Taylor series are used and as a
Taylor series method if a finite orderis used at each step. Taylor series methods are often introduced in numerical analy-sis textbooks and then dismissed practically in the same breath as being too costlyor insufficiently general, but these methods do not actually suffer from those faultswhen properly implemented: see for example [56] which shows how to use them tosolve DAE but also gives many references to papers using them to solve initial-valueproblems (IVP) for ODE. Taylor series methods are of cost polynomial in the numberof bits of accuracy required, on finite intervals, as the number of bits of accuracyrequired goes to infinity [41].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 measured afterwards to validate thesolution. They can be implemented as variable stepsize and variable order methods.In the case of stiff problems they can be implemented as implicit methods [8]. Infact, 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 morenumerically stable for large z . Implicit methods are more effective for stiff problems,as is well-known.The Mathieu equation is not usually stiff , per se, but rather is frequently oscil-latory (especially along the imaginary z -axis, i.e. for modified Mathieu functions).4 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR
This causes a pure Taylor series method to suffer some instability. In contrast, theHermite-Obreschkoff method, which is implicit (something like a generalization of theimplicit midpoint method), performs more satisfactorily.We implemented this, not because we thought it would be the best method toevaluate Mathieu functions, but because we thought it would be a useful project toaid our understanding of them, and provide an independent check on other methodsthat we explore.
This de-scription is quite short because the basic technology is widely understood and ourimplementation is not the main focus of this paper. We provide details only for re-producibility. The key piece is the recurrence relations for the Taylor coefficients ofthe expansion of y ( t n + ∆ t ) = (cid:80) k ≥ w k ∆ t k . Simultaneously, we need the Taylorcoefficients for cos 2 t , which necessitates the Taylor coefficients for sin 2 t . This gives C k +1 = − k + 1 S k S k +1 = 2 k + 1 C k w k +2 = − aw k − q ( k + 1)( k + 2) k (cid:88) j =0 C j w k − j . (5.3) C = cos(2 t n ) and S = sin(2 t n ), while w = y ( t n ) and w = y (cid:48) ( t n ). This recurrenceneeds 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 just Hermite interpolation using the Taylor coefficients at each end of thestep [22]. A string of such blends put together as a piecewise polynomial is calleda “string of blends” because it is tied together at “knots”. Such interpolants arequite remarkably stable numerically, even at very high order, and can be evalu-ated in cost linear in the degree of the interpolant. The high order means thatthey are potentially competitive with spectral methods in efficiency. They are simpleto integrate and differentiate, and reasonably simple to multiply together, and thisgives inexpensive ways to evaluate integrals containing these numerical solutions, e.g. (cid:82) π f ( z )ce ( q, z ) dz . There are also simple methods based on companion matrices tofind zeros of blends; alternatively, there is an iterative scheme of higher order thanNewton’s method (2 + √ ≈ .
732 instead of 2) for the same cost per iteration [20].The Hermite interpolant (“blend”) essentially doubles the order of the methodfor the same number of Taylor coefficients taken at each step, and greatly improvesstability: when the degrees of the Taylor coefficients are the same at each end (as weuse) this is nearly an A -stable method.To take a step, we predict the stepsize by the PI step control suggested in [37],which uses data from the previous two steps (at the beginning we use a predictor basedon 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 (cid:48) ( t n + ∆ t ) = 0 and another(say y s ( t )) with y ( t n + ∆ t ) = 0 and y (cid:48) ( t n + ∆ t ) = 1; we then use collocation at thepoints τ = t n + ∆ t/ τ = t n + 3∆ t/ 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 tomake the residual r ( t ) = z (cid:48)(cid:48) − ( a − q cos( t )) z equal to zero at τ and τ . This gives ustwo linear equations in the two unknowns α and β . Solution of the 2-by-2 system is ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE z (cid:48)(cid:48) ( t ) of the blend z ( t ) is computed by programdifferentiation (actually hand-coded, so not technically “automatic” differentiation).We then measure the residual of the resulting blend at s = t n + ∆ t/
2. This isasymptotically (as ∆ t →
0) the location of the maximal residual of this blend. If thisresidual is small enough, smaller than the specified tolerance, we accept the step anduse this estimate in the PI control to predict the size of the next step.This therefore is merely a special-purpose numerical solver for the Mathieu equa-tion. In our experience it performs well, especially because we can measure the resid-ual r ( t ) separately at many points, to provide reassurance that we have obtained theexact solution to a differential equation very near to the Mathieu equation. Trivially,we have solved y (cid:48)(cid:48) + ( a − q cos 2 t ) y = r ( t )denoting the residual by r ( t ), and we have chosen stepsizes to ensure that | r ( t ) | issmall. By the Alexeev-Gr¨obner nonlinear variation of constants formula this meansthat the global error will also be small: there exists a function G so that (with thesame initial conditions for the reference solution y ( t ) and the computed solution z ( t )solving the above) y ( t ) − z ( t ) = (cid:90) tτ =0 G ( t, τ ) r ( τ ) dτ . Indeed because the Mathieu equation is linear we may write G ( t, τ ) explicitly as aGreen’s function:(5.4) G ( t, τ ) = w I ( τ ) w II ( t ) − w II ( τ ) w I ( t ) . We will use this again later. The Wronskian of the Mathieu equation is 1. We didnot plot G ( t, τ ) but instead merely estimated its effects by integration at differenttolerances. Outside of regions of exponential growth, i.e. where the real part of thecharacteristic exponent is zero, the Mathieu equation is well-conditioned. Even when (cid:60) ( µ ) > relatively small over shortenough integration paths.In short, we wrote a special-purpose numerical ODE solver to evaluate the Math-ieu functions. By varying the initial conditions, we may compute either of the basicsolutions, w I ( t ) or w II ( t ); once one has an eigenvalue given q this already allowscomputation of any of the Mathieu functions. We may compute the other indepen-dent (nonperiodic) solution as we did for figure 8(a). By integrating on a complexpath, we may compute any of the modified Mathieu functions. The modified Math-ieu functions, which are also needed already for the solution of the vibrating ellipticmembrane, are of course not periodic. We can compute Floquet solutions as we didfor figure 8(b), and estimate the Floquet exponent numerically [31] and compare itwith Poincar´e’s formula for µ (equation (3.44)), see figure 10. It is also possible to usethis method with the Green’s function in equation (5.4) to compute the generalizedeigenfunctions needed for double eigenvalues (see the discussion in section 5.4).In spite of these capabilities, we are not recommending our implementation ofthis method as a general-purpose Mathieu solver, both because we have not madeit bulletproof and because we believe that for most purposes it would not be asfast as a spectral solver based on Fourier series—after all, all the Fourier coefficientsare essentially immediately available once the eigenvectors have been computed—butfor our purposes, which included an independent check on existing software, it wasadequate. By varying the precision at which we computed and by taking tolerances6 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR(a) (cid:60) ( µ ) in the real ( a, q ) plane (b) Contours where (cid:60) ( µ ) = 0(c) (cid:60) ( µ ) in the ( a, s ) plane where q = is (d) Contours where (cid:60) ( µ ) = 0 showing doublepoints. Fig. 10 . Stability regions and characteristic exponent µ of solutions of the Mathieu equations.Here we track (cid:60) ( µ ) 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 for integers g ; Regions near the a axis are stable as s isincreased where q = is , except as before only neutrally so at a = g . Near those difficult-to-contourregions we supplemented the graph with Taylor series expansions of the characteristic curves, in red.Double eigenvalues are indicated with blue dots. to be extremely small we were able to buy accurate solutions by paying for computingcycles. That said, it is a variable order method that can tolerate quite high order without numerical difficulty, and so is quite fast. We have not made a detailed speed We have run it routinely with order 20, and sometimes as high as 80. That means that theresidual error on a subinterval of width h is O ( h ), asymptotically as h →
0. But in practice thecode can take quite large steps at this order, so the notion of an asymptotic order is not as helpfulas one might think. What makes this work is that the error coefficients are also very small, and ofcourse the solutions to the problem are analytic.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE
One interesting complication is thatwhen using the Hermite-Obreschkoff or Taylor method (or, indeed, when using any method), accurate solution becomes very expensive for large | z | . This is because solu-tions to the modified Mathieu equation become highly oscillatory along the real axis(so the solutions to the Mathieu equation become highly oscillatory along the imagi-nary axis). Indeed the asymptotics of the solution to the modified Mathieu equationcontain a term exp( i √ q cosh( z )) (see Formula 28.25.1 in the DLMF; we ignore thedenominator, which grows only more slowly). This analysis is confirmed by inspectionof the formulae in Chapter 28.23 of the DLMF (https://dlmf.nist.gov/28.23) where wesee the arguments 2 √ q cosh( z ) and 2 √ q sinh( z ) appearing inside the Bessel functionsused in the expansions for the modified Mathieu functions.Resolving any solution for graphical purposes requires computing a fixed numberof points in each cycle, say 5 or 10. One therefore sees that the number of computedpoints 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. Wetherefore see that the cost to directly compute—in this sense—the solution accuratelymust grow exponentially with the argument. This in part was why the simple Taylorseries at the origin failed already by z = 5 . i , in section 5.1.This does not contradict the theorem in [41] because that only holds on a fixedfinite 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 solutionnot only oscillates but grows doubly exponentially as x grows, containing the termsexp( ± √ q sinh( x ) sin( y )). This fact does not seem to have been explicitly remarkedon elsewhere. That fact suggests that most applications of the Mathieu functions willhave small | z | or else either purely real z or purely imaginary z ; that is, either thereal Mathieu functions or the real modified Mathieu functions will be the relevantfunctions.This simple complexity analysis suggests that any numerical method that we usewill only be helpful and efficient for | z | bounded by a modest constant, say 2 π . Forlarger values of (cid:60) z one must use the asymptotic formula instead. Remark | z | in the Taylor series method, akinto the instability induced for stiff problems [21, 65]. If the stepsize is not aggressivelyreduced as | z | increases, then eventually this instability takes over and the numericalsolution becomes meaningless and rapidly overflows. One would hope that insteadusing an implicit Taylor series method would alleviate this problem, and to a certainextent it does, but even with an implicit method the stepsize is forced to be so smallfor accuracy that failure is assured in practice when | z | is at all large. Remark a is an eigenvalue of the problem and the solution is periodic.Accurate computation of the series coefficients, together with accurate computationof the appropriate Bessel functions, seem to offer an inexpensive way to accuratelyevaluate a modified Mathieu function. Inspection of the derivative, however, showsthat the large argument z = x + iy also needs to be known to exponential accuracy8 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR in order to find an accurate value of the Mathieu function in question, especially neara zero. The difficulty seems to be intrinsic. Since hardly anyone complains aboutthis in the literature, we conclude that most people only want the values of Mathieufunctions and modified Mathieu functions for small or at most moderately large valuesof z .Mathieu himself thought there might be issues for large z and introduced thechange of variable ν = cos( z ), which leads to the following algebraic differentialequation (equation 28.2.3 in the DLMF).(5.5) (cid:0) ν − (cid:1) d dν y ( ν ) + ν ddν y ( ν ) + (cid:0)(cid:0) ν − (cid:1) q − a (cid:1) y ( ν ) = 0 . This algebraic differential equation (and indeed also the similar one in equation 28.2.2in the DLMF that arises on ζ = sin z ) has some interesting computational properties:for one, they are what is called D -finite or holonomic , meaning the recurrence relationfor the Taylor series terms is of fixed order, and which means accurate computationcan be done asymptotically quickly [70]. See also [50, 51, 9].Indeed near the imaginary axis this DE is also numerically easier—one may takemore nearly equal integration steps in the new variable, instead of the ever-decreasingones necessary in the original variable z —but somehow this is just “sweeping theproblem under the rug” because the value of ν itself becomes very large; if the cost ofthe stepsize selection is low (and the control is effective) one should get very nearlythe same numerical performance in the original variable, except that the recurrencerelation for the Taylor series terms in the original variable depends on all previousterms, not just a fixed number as for D-finite functions and so each step is also moreexpensive. However, the presence of the singularities ν = ± However, the method of choice, at least by popularvote, for computation of Mathieu functions is the use of what is effectively a spectralmethod. This is prima facie valid only for the case when a is an eigenvalue and thesolution is periodic. One computes the eigenvalues by the matrix method as above,and then the resulting eigenvector gives the coefficients in the Fourier series expansionof the Mathieu function (if using the continued fraction instead, the recurrence rela-tions can be used, although one has to take numerical care). The Fourier coefficientsdecay extremely rapidly, as is usual for Fourier expansion of smooth functions. Togive 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 whichcoefficient in its Fourier expansion):(5.6) A n m ( q ) A n ( q ) = ( − m ( m !) (cid:16) q (cid:17) m π (cid:0) O (cid:0) m − (cid:1)(cid:1) w II ( π ; a n ( q ) , q ) . This holds as m → ∞ , for fixed n . This states that the m th Fourier coefficientultimately decays like ( q/ m / ( m !) , much faster than exponentially. The other threetypes of Fourier coefficients decay similarly rapidly. This very rapid decay means thatgood approximations can be made with only a few terms of the Fourier series.For the corresponding modified Mathieu functions these Fourier series also con-verge, but now slowly. As an alternative, one can use the same Fourier coefficients in
ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE
39a Bessel function expansion, for instance equation 28.23.2 of the DLMF:me ν (cid:0) , h (cid:1) M ( j ) ν ( z, h ) = ∞ (cid:88) n = −∞ ( − n c ν n ( h ) C ( j ) ν +2 n (2 h cosh z ) , where the modified Mathieu function on the left can be approximated by the serieson the right; the notation therein is different to the notation in our paper, but the C s are Bessel functions and the c n s are the Fourier coefficients for the ordinary evenperiod- π Mathieu functions. These are not the only series one might use. In [30] (andin [1] and the DLMF) we find the following expansions, which the authors claim arerapidly convergent:Ce n ( q, x ) = ( − n A (cid:114) π (cid:88) k ≥ ( − k A k J k ( s ) J k ( t )Ce n +1 ( q, x ) = ( − n A (cid:114) π (cid:88) k ≥ ( − k A k +1 ( J k +1 ( s ) J k ( t ) + J k ( s ) J k +1 ( t ))Se n ( q, x ) = ( − n B (cid:114) π (cid:88) k ≥ ( − k B k ( J k +1 ( s ) J k − ( t ) − J k − ( s ) J k +1 ( t ))Se n +1 ( q, x ) = ( − n B (cid:114) π (cid:88) k ≥ ( − k B k +1 ( J k +1 ( s ) J k ( t ) − J k ( s ) J k +1 ( t )) . (5.7)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 Mathieufunction. 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 backat least to [24]. Erd´elyi thought the ‘coincidence’ of these series coefficients to besignificant [27]. Here we have translated to the notation of this paper—although thenormalization used in the above does not agree with the normalization here, eventhough the authors of [30] claim to use the same normalization that we do; insteadthe formulas are the same as those in [28]. We tried this, and it worked well. Indeed,such a Fourier-Bessel method is similar to those advocated by almost everyone, from[4] to [76]. Against this Fourier-Bessel method it is not clear that a straightforwardnumerical solution such as the one we have implemented would be competitive, butfor instance the authors of [63] seem to think that something like it might be. Severalauthors indeed claim that numerical evaluation of Bessel functions for large argumentsand high order is difficult or inefficient, but we do not believe this statement: inour experience the standard recurrence method performs well. There do seem to benumerical instabilities in some Bessel series that can cause difficulty, although thesecan be mitigated by careful choices among the expansions [69]. We have not done adetailed comparison of methods.In general one has to be a bit careful with rounding error if the recurrence rela-tions are used to compute the Fourier coefficients: the relations can be unstable, asnoted by many authors. A practical solution is to use forward recurrence for the firstfew and backward recurrence for the rest. This seems to have been first advocated byBlanch [11]. Similar considerations apply if numerical eigenvectors are used (after all,the matrix is simply an arrangement of the recurrence relations). A rule of thumb is0 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR that if one is working to d digits with the forward recurrence, then the Fourier coeffi-cients will decay down to about 10 − d/ times the magnitude of the largest coefficientand then rounding error will start to impact the results thereafter. Four ways suggest themselvesto compute the generalized eigenfunction u = ∂y/∂a at a double point q ∗ with doubleeigenvalue a ∗ , which is a solution of equation (4.7) that satisfies periodic boundaryconditions. We duplicate that equation here for convenience:(5.8) u (cid:48)(cid:48) + ( a ∗ − q ∗ cos 2 α ) u + y ∗ = 0 . The first and simplest way is undoubtedly what people actually use: one pretendsthat the eigenvalue is not actually a double one—typically it would have split into a ∗ + d √ q − q ∗ + · · · and a ∗ − d √ q − q ∗ + · · · where q is a floating-point approximationto q ∗ anyway—and then use the computed eigenfunctions from the matrix method,each with norm O ( √ q − q ∗ ) and simply live with the errors. That does not soundlike professional practice, but if it is done knowingly then we suspect that it willusually give perfectly reasonable answers. If done unknowingly then we disapprove,but they’ll likely get away with it.The second way is to compute a generalized eigenvector of the infinite tridiagonalmatrix A for the problem. This is scarcely harder than the crude approach above:first, one averages the two computed approximate eigenvectors for the double eigen-value split pair, and averages the eigenvalues, to get a more accurate double eigenvalue a and eigenvector c . For convenience in the exposition below, imagine that we haveput the pair of eigenvalues that arise on splitting the double eigenvalue numerically inpositions 1 and 2, and similarly put their associated eigenvectors in columns 1 and 2of the matrix of eigenvectors. If by some miracle the eigenvalue routine doesn’t splitthe double and instead produces a single, accurate, double eigenvalue and only onecorresponding eigenvector, then we use that.Next, one solves ( a I − A ) u = − c . This system is singular, but c is in its range and this is not difficult; one could use theSVD, for instance . Since the matrix is of low dimension, the expense of the SVD isno obstacle. Then one uses the entries of u to construct the generalized eigenfunction u ( z ) in the same manner one constructs the eigenfunction v ( z ) from the eigenvector c . To expand a given function f ( z ) as a sum of these eigenfunctions, notice thatthe norm of v is zero, but the inner product of v ( z ) with u ( z ) is nonzero. Thegeneralized eigenfunction and v ( z ) are each orthogonal to all other eigenfunctions,however. Put(5.9) f ( z ) = αv ( z ) + βu ( z ) + N (cid:88) j =3 γ k v k ( z ) . Then the γ k are easily found by orthogonality as usual: γ k = (cid:82) pz =0 f ( z ) v k ( z ) dz (cid:82) pz =0 v k ( z ) dz k = 3 , , . . . , N . It is straightforward to set up the recurrence relations for this generalized eigenvector: they aremerely forced versions of the recurrences in equations 3.20–3.36; but as mentioned these recurrencesare known to be unstable sometimes, while in contrast the SVD has the virtue of answering thequestion in a manner that relieves all doubts.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE not use the conjugate. The norm of anonzero function can be negative, complex, or indeed zero. Such norms are called indefinite norms, in for instance [5].To find α and β we use the fact that while the norm of v ( z ) is zero, the norm ofthe generalized eigenfunction u ( z ) is not zero, and also (cid:90) p v ( z ) u ( z ) dz (cid:54) = 0 . Thus the two equations (cid:90) π f ( z ) v ( z ) dz = α · β (cid:90) π u ( z ) v ( z ) dz (cid:90) π f ( z ) u ( z ) dz = α (cid:90) π v ( z ) u ( z ) dz + β (cid:90) π u ( z ) dz (5.10)give us a triangular two-by-two system (indeed with constant diagonal) to solve forthe unknown coefficients.For example, consider q = 1 . i , the Mulholland-Goldstein doublepoint again, and its associated eigenvalue a = 2 . N by N matrix, where we took N = 25:its split eigenvalues were 2 . ± . · − i ), wherewe have displayed the distinct real digits in red. We averaged the correspondingeigenvectors, and put v ( z ) = c √ (cid:88) k =2 c k cos(2( k − z ) . The √ √ v (0) = 1 by scaling.We index from 1 in the above equation because that is usual for matrices. Thereal and imaginary parts of this are plotted in figure 11(a); v ( z ) is, of course, anapproximation for ce ( q, z ) = ce ( q, z ), the coalesced eigenfunction. By examining itsresidual v (cid:48)(cid:48) + ( a − q cos 2 z ) v we see that it is accurate to 10 − (plot not shown).We verified that v ( z ) has numerically zero norm, as well: (cid:90) π v ( z ) dz ≈ − . · − − . · − i . We then use the SVD to solve the singular system ( a I − A ) u = − c , and form(5.11) u ( z ) = u √ n (cid:88) k =2 u k cos(2( k − z ) . We then removed a multiple of v ( z ) so that u (0) = 0. The real and imaginary partsof this generalized eigenfunction are plotted in figure 11(b).We now use these eigenfunctions to expand a smooth function of period π . Wetook(5.12) f ( z ) = e cos 2 z cos 6 z ≈ αv ( z ) + βu ( z ) + (cid:88) k =3 γ k v k ( z ) . C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR
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 (cid:90) π f ( z ) v ( z ) dz = α · β (cid:90) π u ( z ) v ( z ) dz we may identify β ≈ . . i . Now because (cid:90) π f ( z ) u ( z ) dz = α (cid:90) π v ( z ) u ( z ) dz + β (cid:90) π u ( z ) dz where both of the integrals on the right are nonzero and we now know β , this gives α ≈ . − . i . As you can see from the graph of the magnitudes of thecomputed γ k for 3 ≤ k ≤
25 in figure 12(a), these are appreciable. Notice also theexponential decay of the coefficients; this is typical for a spectral expansion of a smoothfunction. We plot the error f ( z ) − αv ( z ) − βu ( z ) − (cid:80) k ≥ γ k v k ( z ) in figure 12(b) wherewe see that it is less than 5 · − .The third way that we were thinking that one could compute these eigenfunctionsseems a little more complicated: we solve the initial-value problem sweeping forwardfor v ( z ) using the Hermite-Obreschkoff method described earlier, which chooses themesh; we record the local Taylor series for the two local functions satisfying y ( z j ) = 1, y (cid:48) ( z j ) = 0 and y ( z j ) = 0, y (cid:48) ( z j ) = 1. We then solve the boundary value problem for u ( z ) on that interval by imposing periodic boundary conditions and using collocationat two points in each interval. If there were M subintervals, this gives an almost blockdiagonal matrix of 2 M equations in the 2 M unknowns, namely the coefficients α k and β k of the linear combination of the two solutions at each interior node, togetherwith α = α M and β = β M . This sounds involved, but in fact it is straightforward:the truly hard part of solving BVP by collocation (see [6, 7]) is choosing the mesh,and here we don’t have that problem because the mesh will have been chosen by thePI control on the forward sweep. We suspect that few people will actually implementthis method, however; we haven’t, yet, either.A fourth way, similar but perhaps even simpler, is to use the Green’s functionfrom equation (5.4); we already have methods for integrating strings of blends andfor 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 atthe same time. The general solution of equation (5.8) is (suppressing the dependenceof w I ( z ; a, q ) on a and q and similarly w II for brevity) u ( z ) = αw I ( z ) + βw II ( z ) + w II ( z ) (cid:90) z w I ( ζ ) y ∗ ( ζ ) dζ − w I ( z ) (cid:90) z w II ( ζ ) y ∗ ( ζ ) dζ . If we are solving for a generalized eigenfunction for ce g ( q ∗ , z ) then y ∗ = ce g ( 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 = π u ( π ) = βw II ( π ) − w I ( π ) (cid:90) π w II ( ζ ) y ∗ ( ζ ) dζ . Since w I ( π ) = 0 because it, being the eigenfunction in question in this example, isperiodic, we see that β = 0 as well. Thus all we need to identify the eigenfunction isthe integral against the Green’s function. We tried this, and it confirmed the resultsin the top row of figure 11. The bottom row in that figure was computed this way. Each collocation point will give an equation involving four unknowns, the α s and β s of theendpoints of the interval containing the collocation point.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE (a) coalesced eigenfunctions (b) generalized eigenfunction(c) coalesced eigenfunctions (d) generalized eigenfunction Fig. 11 . Top row left: Real and imaginary parts of the coalesced eigenfunctions v ( z ) =ce ( q, z ) = ce ( q, z ) corresponding to the Mulholland-Goldstein double point q ≈ . i (real partin black, imaginary part in red). On the right, we have the corresponding generalized eigenfunctionobtained by solving y (cid:48)(cid:48) + ( a − q cos 2 z ) y + v = 0 . Bottom row left: Real and imaginary partsof the coalesced eigenfunctions v ( z ) = se ( q, z ) = se ( q, z ) corresponding to the next-largest pureimaginary double point q = 6 . . . . i with eigenvalue approximately . . On the right, wehave the corresponding generalized eigenfunction obtained by solving y (cid:48)(cid:48) + ( a − q cos 2 z ) y + v = 0 .
6. Concluding remarks.
We have completed a historical survey of the compu-tation of the Mathieu functions, which are defined to be the period π and period 2 π solutions of the Mathieu differential equation (3.1). Our original motivation for thiswas to solve a problem of pulsatile blood flow in a vessel of elliptic cross-section; toactually do that we used the spectral method, but with the brute force of multipleprecision to power through the double-eigenvalue difficulty. That was inelegant, andso we investigated the alternatives.To carry out our investigation, we implemented our own procedures (in Maple)for the computation of Mathieu functions (and generalized eigenfunctions) in order4 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR(a) Decay of coefficients (b) Approximation error
Fig. 12 . (Left) size of coefficients in the expansion of f ( z ) = exp(cos 2 z ) cos 6 z in terms ofMathieu functions at the Mulholland-Goldstein double point q = 1 . . . . i ; the coefficients seemto decay like (blue dashed line) exp(1 . − . k ) . (Right) the difference f ( z ) − αv ( z ) − βu ( z ) − (cid:80) k =3 γ k v k ( z ) . Here α ≈ . − . i and β ≈ . − . i . (black for real part,red for imaginary part). to explore some of the difficulties. Our code is at present, like we imagine Blanch’sto have been, “artesanal” and meant only for use with careful human supervision.For experimentation, of course, this is a feature, not a bug. The task of constructingfully general, bulletproof code for the Mathieu functions is one that calls for dedicatedeffort and analysis. We hope that this present paper encourages a team to undertakethe task.In 1957 G. Temple called Mathieu functions “indispensable but intractable in-struments of mathematical physics”; after this investigation we no longer think thatthey are so intractable. But they are somewhat involved, and the major remainingquestion is how best to compute the sometimes-needed generalized eigenfunctions.This suggests that it might be useful to think of Mathieu functions together with thegeneralized eigenfunctions as a frame for spectral expansion [2]. We suspect that thegeneralized eigenvector approach will provide a practical frame. Of course, we mustnot forget about the advances in direct numerical solution of the underlying PDE: themethods of [33] may be preferable in many applications over expansion in Mathieufunctions.Perhaps the most interesting results of our readings of the literature were, first,the discovery that Mathieu had anticipated Lindstedt’s anti-secularity perturbationmethod by several decades ; and that Blanch’s algorithm for computing eigenvaluescould be implemented (easily!) in series, thereby allowing one to compute seriesfor eigenvalues which would allow greater surety in continuation methods (Newtoniteration starting with estimates from nearby q is usually used) or to compute Puiseuxseries about double points. We believe that the computation of these Puiseux series isnew to this paper. The notion of generalized eigenfunctions for the Mathieu equation In [71] we find that in the Astronomy literature some people called the Mathieu equation the“Lindstedt equation”, so the ideas of Mathieu and of Lindstedt may be more connected than weknow.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE
Acknowledgments.
We were grateful for advice from John May, which spedup our Maple implementation of Mathieu’s anti-secular perturbation method. Wealso thank Erik Postma and J¨urgen Gerhard for many comments on the computationof Mathieu functions in Maple, and Marcus Webb at Manchester for discussions oncompleteness. The Department of Applied Mathematics provided useful support forthis project, as did the Rotman Institute of Philosophy at Western. RMC thanks theIsaac Newton Institute for Mathematical Sciences and the staff of both the UniversityLibrary and the Betty and Gordon Moore Library at Cambridge for support and hos-pitality during the programme Complex Analysis: Tools, techniques, and applications,when some of the work on this project was undertaken.
REFERENCES[1]
M. Abramowitz and I. A. Stegun , Handbook of mathematical functions: with formulas,graphs, and mathematical tables , vol. 55, Courier Corporation, 1964.[2]
B. Adcock and D. Huybrechs , Frames and numerical approximation , SIAM Review, 61(2019), pp. 443–473.[3]
A. Aitken , Dr. EL Ince , Nature, 148 (1941), pp. 309–310.[4]
F. A. Alhargan , A complete method for the computations of Mathieu characteristic numbersof integer orders , SIAM review, 38 (1996), pp. 239–255.[5]
P. Arbenz and M. E. Hochstenbach , A Jacobi–Davidson method for solving complex sym-metric eigenvalue problems , SIAM Journal on Scientific Computing, 25 (2004), pp. 1655–1673.[6]
U. Ascher, J. Christiansen, and R. D. Russell , COLSYS–a collocation code for boundary-value problems , in Codes for Boundary-Value problems in ordinary differential equations,Springer, 1979, pp. 164–185.[7]
G. Bader and U. Ascher , A new basis implementation for a mixed order boundary value odesolver , SIAM journal on scientific and statistical computing, 8 (1987), pp. 483–500.[8]
D. Barton , On Taylor series and stiff equations , ACM Transactions on Mathematical Software(TOMS), 6 (1980), pp. 280–294.[9]
A. Benoit, M. Joldes¸, and M. Mezzarobba , Rigorous uniform approximation of D-finitefunctions using Chebyshev expansions , Mathematics of Computation, 86 (2017), pp. 1303–1341.[10]
W. Bickley , The tabulation of Mathieu functions , Mathematical Tables and Other Aids toComputation, 1 (1945), pp. 409–419.[11]
G. Blanch , On the computation of Mathieu functions , Journal of Mathematicsand Physics, 25 (1946), pp. 1–20, https://doi.org/10.1002/sapm19462511, https://onlinelibrary.wiley.com/doi/abs/10.1002/sapm19462511, https://arxiv.org/abs/https://onlinelibrary.wiley.com/doi/pdf/10.1002/sapm19462511.[12]
G. Blanch , The asymptotic expansions for the odd periodic Mathieu functions , Transactionsof the American Mathematical Society, 97 (1960), pp. 357–366.[13]
G. Blanch , Numerical evaluation of continued fractions , Siam Review, 6 (1964), pp. 383–421.[14]
G. Blanch , Numerical aspects of Mathieu eigenvalues , Rendiconti del Circolo Matematico diPalermo, 15 (1966), pp. 51–97.[15]
G. Blanch and D. Clemm , The double points of Mathieus differential equation , Mathematicsof Computation, 23 (1969), pp. 97–108.[16]
G. Blanch, A. Lowan, R. Marshak, and H. Bethe , The internal temperature-density dis-tribution of the sun. , The Astrophysical Journal, 94 (1941), p. 37.[17]
E. Bolmont, P. Nabonnand, and L. Rollet , Les ambitions Parisiennes con-trari´ees d’ ´Emile Mathieu (1835–1890) , 2015, https://images.math.cnrs.fr/Les-ambitions-parisiennes-contrariees-d-Emile-Mathieu-1835-1890.html.[18]
C. J. Bouwkamp , A note on Mathieu functions , Proc. Kon. Nederland. Akad. Wetensch. v51,(1948), pp. 891–893.[19]
L. Chaos-Cador and E. Ley-Koo , Mathieu functions revisited: matrix evaluation and gen-erating functions , Revista mexicana de f´ısica, 48 (2002), pp. 67–75. C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR[20]
R. M. Corless , Inverse cubic iteration , arXiv preprint arXiv:2007.06571, (2020).[21]
R. M. Corless and N. Fillion , Backward error analysis for perturbation methods , in Algo-rithms and complexity in mathematics, epistemology, and science, Springer, 2019, pp. 35–79.[22]
R. M. Corless and E. Postma , Blends in Maple , arXiv preprint arXiv:2007.05041, (2020).[23]
R. M. Corless and L. Rafiee Sevyeri , The Runge example for interpolation and Wilkinsonsexamples for rootfinding , vol. 62, SIAM, 2020, pp. 231–243.[24]
J. Dougall , The solution of Mathieu’s differential equation , Proceedings of the EdinburghMathematical Society, 34 (1915), pp. 176–196.[25]
B. Duan and M. Zamir , Approximate solution for pulsatile flow in tubes of slightly noncircularcross-sections , Utilitas Mathematica, 40 (1991), pp. 13–26.[26]
P. Duhem , ´Emile Mathieu, his life and works , Bulletin of the American Mathematical Society,1 (1892), pp. 156–168.[27] A. Erd´elyi , On certain expansions of the solutions of Mathieu’s differential equation , in Math-ematical Proceedings of the Cambridge Philosophical Society, vol. 38, Cambridge Univer-sity Press, 1942, pp. 28–33.[28]
A. e. a. Erd´elyi , Higher Transcendental Functions , McGraw-Hill, 1953.[29]
D. Erricolo , Acceleration of the convergence of series containing Mathieu functions us-ing Shanks transformation , IEEE Antennas and Wireless Propagation Letters, 2 (2003),pp. 58–61.[30]
D. Erricolo and G. Carluccio , Algorithm 934: Fortran 90 subroutines to compute Math-ieu functions for complex values of the parameter , ACM Transactions on MathematicalSoftware (TOMS), 40 (2013), p. 8.[31]
T. F. Fairgrieve and A. D. Jepson , OK Floquet multipliers , SIAM journal on numericalanalysis, 28 (1991), pp. 1446–1462.[32]
D. Frenkel and R. Portugal , Algebraic methods to compute Mathieu functions , Journal ofPhysics A: Mathematical and General, 34 (2001), p. 3541.[33]
M. J. Gander and H. Zhang , A class of iterative solvers for the Helmholtz equation: factor-izations, sweeping preconditioners, source transfer, single layer potentials, polarized traces,and optimized schwarz methods , SIAM Review, 61 (2019), pp. 3–76.[34]
K. O. Geddes, S. R. Czapor, and G. Labahn , Algorithms for computer algebra , SpringerScience & Business Media, 1992.[35]
S. Goldstein , Lam´esche, Mathieusche und verwandte Funktionen in Physik und Technik. byM.J.O. Strutt. , The Mathematical Gazette, 17 (1933), p. 5960, https://doi.org/10.2307/3607965.[36]
D. A. Grier , Gertrude Blanch of the mathematical tables project , IEEE Annals of the Historyof Computing, 19 (1997), pp. 18–27.[37]
K. Gustafsson, M. Lundh, and G. S¨oderlind , A PI stepsize control for the numerical solu-tion of ordinary differential equations , BIT Numerical Mathematics, 28 (1988), pp. 270–287.[38]
J. C. Guti´errez-Vega, R. Rodrıguez-Dagnino, M. Meneses-Nava, and S. Ch´avez-Cerda , Mathieu functions, a visual approach , American Journal of Physics, 71 (2003), pp. 233–242.[39]
E. Heine , Handbuch der Kugelfunktionen I, II , Berlin: Reimer, 1878.[40]
C. Hunter and B. Guerrieri , The eigenvalues of Mathieu’s equation and their branch points ,Studies in Applied Mathematics, 64 (1981), pp. 113–141.[41]
S. Ilie, G. S¨oderlind, and R. M. Corless , Adaptivity and computational complexity in thenumerical solution of odes , Journal of Complexity, 24 (2008), pp. 341–361.[42]
E. Ince , A proof of the impossibility of the coexistence of two Mathieu functions , in Proc.Camb. Phil. Soc, vol. 21, 1922, pp. 117–120.[43]
E. Ince , The second solution of the Mathieu equation , in Mathematical Proceedings of theCambridge Philosophical Society, vol. 23, Cambridge University Press, 1926, pp. 47–49.[44]
E. Ince , The Mathieu equation with numerically large parameters , Journal of the LondonMathematical Society, 1 (1927), pp. 46–50.[45]
E. Ince , Xxii.tables of the elliptic-cylinder functions , Proceedings of the Royal Society ofEdinburgh, 52 (1933), pp. 355–423.[46]
M. Levi , Stability of the inverted penduluma topological explanation , SIAM review, 30 (1988),pp. 639–644.[47]
S. Lubkin and J. Stoker , Stability of columns and strings under periodically varying forces ,Quarterly of Applied Mathematics, 1 (1943), pp. 215–236.[48] ´E. Mathieu , M´emoire sur le mouvement vibratoire d’une membrane de forme elliptique. , Jour-nal de math´ematiques pures et appliqu´ees, 13 (1868), pp. 137–203.[49]
J. Meixner, F. W. Sch¨afke, and G. Wolf , Mathieu functions , Springer, 1980.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE [50] M. Mezzarobba , NumGfun: a package for numerical and analytic computation with D-finitefunctions , in Proceedings of the 2010 International Symposium on Symbolic and AlgebraicComputation, 2010, pp. 139–145.[51]
M. Mezzarobba , A note on the space complexity of fast D-finite function evaluation , in Inter-national Workshop on Computer Algebra in Scientific Computing, Springer, 2012, pp. 212–223.[52]
Y. Miyazaki, N. Asai, Y. Kikuchi, D. Cai, and Y. Ikebe , Computation of multiple eigenval-ues of infinite tridiagonal matrices , Mathematics of computation, 73 (2004), pp. 719–730.[53]
P. M. Morse and H. Feshbach , Methods of theoretical physics , American Journal of Physics,22 (1954), pp. 410–413.[54]
H. Mulholland and S. Goldstein , The characteristic numbers of the Mathieu equation withpurely imaginary parameter , The London, Edinburgh, and Dublin Philosophical Magazineand Journal of Science, 8 (1929), pp. 834–840.[55]
D. Naylor , On simplified asymptotic formulas for a class of Mathieu functions , SIAM journalon mathematical analysis, 15 (1984), pp. 1205–1213.[56]
N. S. Nedialkov and J. D. Pryce , Solving differential-algebraic equations by Taylor series(i): Computing Taylor coefficients , BIT Numerical Mathematics, 45 (2005), pp. 561–591.[57]
L. Onsager , Solutions of the Mathieu Equation of Period 4 Pi and Certain Related Functions ,PhD thesis, Yale University, 1935.[58]
B. N. Parlett , The Rayleigh quotient iteration and some generalizations for nonnormal ma-trices , Mathematics of Computation, 28 (1974), pp. 679–693.[59]
A. P. Prudnikov, Yu. A. Brychkov, and O. I. Marichev , Integrals and Series: More SpecialFunctions, Vol. 3 , Gordon and Breach Science Publishers, New York, 1990. Translatedfrom the Russian by G. G. Gould. Table erratum Math. Comp. v. 65 (1996), no. 215, p.1384.[60]
J. D. Pryce , Numerical solution of Sturm-Liouville problems , Oxford University Press, 1993.[61]
H. Rubin , Anecdote on power series expansions of Mathieu functions , Journal of Mathematicsand Physics, 43 (1964), pp. 339–341.[62]
B. Salvy and P. Zimmermann , Gfun: a Maple package for the manipulation of generatingand holonomic functions in one variable , ACM Transactions on Mathematical Software(TOMS), 20 (1994), pp. 163–177.[63]
M. Schneider and J. Marquardt , Fast computation of modified Mathieu functions appliedto elliptical waveguide problems , IEEE Transactions on Microwave Theory and Techniques,47 (1999), pp. 513–516.[64]
J. Shen and L.-L. Wang , On spectral approximations in elliptical geometries using Mathieufunctions , Mathematics of Computation, 78 (2009), pp. 815–844.[65]
G. S¨oderlind, L. Jay, and M. Calvo , Stiffness 1952–2012: Sixty years in search of a defini-tion , BIT Numerical Mathematics, 55 (2015), pp. 531–558.[66]
G. F. J. Temple , Edmund Taylor Whittaker, 1873-1956 , 1956.[67]
H. Tropp , Interview with Gertrude Blanch , 1973, https://sova.si.edu//details/NMAH.AC.0196
M. Trott , The Mathematica guidebook for symbolics , Springer Science & Business Media,2007.[69]
A. Van Buren and J. Boisvert , Accurate calculation of the modified Mathieu functions ofinteger order , Quarterly of applied mathematics, 65 (2007), pp. 1–23.[70]
J. van der Hoeven , Fast evaluation of holonomic functions near and in regular singularities ,Journal of Symbolic Computation, 31 (2001), pp. 717–744.[71]
B. van der Pol and M. J. O. Strutt , Ii. on the stability of the solutions of Mathieu’s equation ,The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 5(1928), pp. 18–38, https://doi.org/10.1080/14786440108564441, https://doi.org/10.1080/14786440108564441, https://arxiv.org/abs/https://doi.org/10.1080/14786440108564441.[72]
E. Whittaker , On the functions associated with the elliptic cylinder in harmonic analysis , inProceedings of the Fifth International Congress of Mathematicians, vol. 1, 1912, pp. 366–371.[73]
E. Whittaker , Edward Lindsay Ince: 18911941 , Journal of the London Mathematical Society,1 (1941), pp. 139–144.[74]
E. T. Whittaker and G. N. Watson , A course of modern analysis , Cambridge Univ. Press,1927.[75]
M. Zamir , Hemo-dynamics , Springer, 2016.[76]
C. Ziener, M. R¨uckl, T. Kampf, W. Bauer, and H. Schlemmer , Mathieu functions forpurely imaginary parameters , Journal of Computational and Applied Mathematics, 236(2012), pp. 4513–4524. C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR
Appendix A. Pulsatile flow in tube of elliptic cross-section.
The notationin this appendix differs slightly from that of the main paper, and instead agrees withthat of our other study (currently in progress).The equation governing the oscillatory component of the flow takes the generalform [75](A.1) ∂u∂t + 1 ρ ∂p∂z = µρ (cid:18) ∂ u∂x + ∂ u∂y (cid:19) , where x, y are rectangular coordinates within the cross section of the tube, z is alongthe axis of the tube and p is pressure. It is at this point that the geometry of the crosssectional boundary of the tube dictates the choice of coordinate system in which thegoverning equation is to be solved. In the case of a tube of circular cross section , the governing equation (A.1)takes the form(A.2) ∂u c ∂t + 1 ρ ∂p∂z = µρ (cid:18) ∂ u c ∂r + 1 r ∂u c ∂r (cid:19) , where subscript ‘ c ’ is being used to associate the results with a tube of circular crosssection.For an oscillatory pressure gradient of the form(A.3) ∂p∂z = k e iωt and by separation of variables(A.4) u c ( r, t ) = U c ( r ) e iωt the equation becomes(A.5) d U c dr + 1 r dU c dr − i Λ c a U c = k µ , where(A.6) Λ c = ρωa µ is a nondimensional frequency parameter and a is the radius of the tube.Equation (A.5) is a form of a Bessel equation with general solution(A.7) U c ( r ) = ik a µ Λ c + AJ ( ζ ) + BY ( ζ )where A and B are arbitrary constants and J and Y are Bessel functions of orderzero and of the first kind, respectively, satisfying the standard Bessel equations d J dζ + 1 ζ dJ dζ + J = 0(A.8) d Y dζ + 1 ζ dY dζ + Y = 0 . (A.9) ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE ζ is related to the radial coordinates by(A.10) ζ ( r ) = Ω ra , where Ω is a frequency parameter related to the nondimensional frequency parameter(A.11) Ω = (cid:18) i − √ (cid:19) (cid:112) Λ c . In the case of a tube of elliptic cross section the boundary conditions dictatea transformation to elliptic coordinates(A.12) x = d cosh ξ cos η, y = d sinh ξ sin η where 2 d is the focal distance and ξ , η are the elliptic coordinates, and equation (A.1)becomes(A.13) ∂u e ∂t + 1 ρ ∂p∂z = µρ d (cosh 2 ξ − cos 2 η ) (cid:18) ∂ u e ∂ξ + ∂ u e ∂η (cid:19) , where the subscript ‘ e ’ is now used to associate the results with a tube of elliptic crosssection.For an oscillatory pressure gradient of the form(A.14) ∂p∂z = k e iωt we use separation of variables(A.15) u e ( ξ, η, t ) = w ( ξ, η ) e iωt so that equation (A.13) can be formulated as an inhomogeneous Helmoltz equation(A.16) 2 d (cosh 2 ξ − cos 2 η ) (cid:18) ∂ u e ∂ξ + ∂ u e ∂η (cid:19) − iρωµ w = k µ . Using the translation(A.17) w ( ξ, η ) = v ( ξ, η ) − k iρω , the inhomogeneous term of equation (A.16) becomes(A.18) (cid:18) ∂ v∂ξ + ∂ v∂η (cid:19) − i e (cosh 2 ξ − cos 2 η ) v = 0where(A.19) Λ e = ρωd µ is the elliptic equivalent of the nondimensional frequency parameter.Applying separation of variables using(A.20) v ( ξ, η ) = f ( ξ ) g ( η ) , C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR equation (A.18) yields two separate Mathieu equations d gdη + ( s + 2 q cos 2 η ) g = 0(A.21) d fdξ − ( s + 2 q cosh 2 ξ ) f = 0(A.22)where s is a separating constant and(A.23) q = i Λ e . The foregoing analysis demonstrates clearly that the boundary and boundaryconditions play a key role in the transition from Bessel equations to Mathieu equationsin the case of pulsatile flow in a tube. However, while this explains the mathematicalaspects and origin of the transition from Bessel to Mathieu equations in this particularapplication, it does not provide a hint as to the corresponding origin of this transitionin the context of the physics of the flow. One of the aims of our study was tounderstand the link between the mathematical and the physical aspects and origin ofthis transformation.While, at first, properties of the flow in a slightly elliptic tube appear as a smallperturbation of the flow in a circular tube (and can indeed be analyzed in this man-ner [25]), a closer look at the shear stress on the boundary of the tube shows a moresignificant change that occurs as the cross section of the tube changes from circularto elliptic as shown in Figure 13. The figure shows that no matter how small thechange from circular to elliptic cross section is, the two vertices of the ellipse leadto periodicity in the distribution of shear stress on the boundary. This periodicity isthe origin of the transition from Bessel to Mathieu equations in the description of theflow. Appendix B. Confocal Ellipses verses Fixed Aspect Ratio Ellipses andBessel functions as a limit of Mathieu functions .
To help people visualizethe difference between a family of confocal ellipses, all of which have the same setof foci, and a family of fixed aspect ratio ellipses, which is likely what most peopleimagine when they think of a family of ellipses, we made figure 14. The point is todemonstrate how quickly confocal ellipses become circular-looking.Algebraically, the aspect ratio of a confocal ellipse is a/b = tanh( β ) where the con-focal parameterization is x = c cosh( β ) cos( α ), y = c sinh( β ) sin( α ). Asymptotically, a/b ≈ − − β ) so we see that a confocal family becomes circular exponentiallyquickly.This helps in understanding the limiting case in which the solutions of the modifiedMathieu equation approach Bessel functions. We modify the treatment slightly in [74]in what follows. We start from the Mathieu equation but replace q = k / − k /
4) to get(B.1) d ydz + ( a − k cos(2 z ) y = 0 , so k = √ q = 2 h . We change variables with ξ = k sin( z ), so that with M = a − k / A vertex of an ellipse is an endpoint of the major axis; a co-vertex is an endpoint of the minoraxis.ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE Fig. 13 . Shear stress at the boundary of a tube in pulsatile flow. The curves represent a rangeof tubes with different cross sections ranging from circular to increasingly elliptic. The numbers onthe 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 ofthe shear stress in a tube of circular cross section. we get, using cos 2 z = 1 − z = 1 − ξ /k ,(B.2) ξ d ydξ + ξ dydξ − ( ξ + M ) y ( ξ ) = − k d ydξ . If k is small, we recognize this as a small (admittedly singular) perturbation of Bessel’sequation, so the outer solution, away from ξ = 0, will be y ( xi ) = C I M ( ξ ) + C K M ( ξ ) + O ( k ) , where I and K denote Bessel functions. The next term, which is O ( k ), in the pertur-bation expansion is somewhat complicated, though small in size when we computedit, but in fact we don’t need it. The result is plotted in figure 15 and compared withthe relevant Mathieu function. Now, as k →
0, for ξ to remain O (1) we must havesin( z ) → ∞ ; that is, this matching will only be valid for large imaginary z . This isthe 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 / z = 2; and indeed already by q = 1 / ( q, iη ) (that is, we use a = − . . . . ) we see a marked resemblance to (cid:60) ( I M ( ξ )) in figure 15. Indeed if we choose a point ( η = 2 .
5) and choose a normaliza-tion factor C = 0 . C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR(a) Confocal ellipses and a circle (b) Fixed aspect ellipses and a circle
Fig. 14 . A comparison of confocal ellipses (left) with fixed aspect ratio ellipses (right). Theconfocal ellipses have foci at ± / . Their parametric equations are x = c cosh( β ) cos( α ) , y = c sinh( β ) sin( α ) . The four ellipses shown have β = [0 , / , / , . The circle (red dashed line)was chosen to have radius r = c exp(2) / , which as the mean of the semi-major and semi-minoraxes of the largest ellipse, corresponds well. We see that the largest confocal ellipse shown here isappreciably circular. Larger β would generate even more circular-looking ellipses. In contrast, thefixed aspect-ratio ellipses all have the same aspect ratio as the smallest nonsingular confocal ellipseshown on the left, namely a/b = tanh(2 / ≈ . . then the error ce ( q, iη ) − C (cid:60) ( I M ( k sin( iη ))) gets very small very quickly as shownin figure 15(b). Notice that here k ≈ . Appendix C. Comparing Mathieu’s perturbation series to Maple’s.
Inorder to compare the q -series produced by Maple’s series command with the handcomputation of Mathieu, we have to ensure that the same normalization is enforced.Since Mathieu’s normalization was to make the coefficient of cos gα equal to unity,and Maple’s normalization is to make (cid:82) p ce g ( α ) dα = π except when g = 0 when it is2 π , we compute Maple’s series and then divide by the coefficient of cos gα . Similarlyfor the se g ( α ) functions.For ce g ( q, t ) with symbolic g we get, on dividing by the coefficient of cos gt , thefollowing terms in the series. The coefficient of q is, provided g > g − α g − − cos ( g + 2) α g + 1) . This agrees perfectly with what Mathieu had. The coefficient of q is, provided g > g − α )32 ( g −
2) ( g −
1) + cos (( g + 4) α )32 ( g + 2) ( g + 1)Again this agrees perfectly with what Mathieu had. The coefficient of q is, provided ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE (a) ce ( q, iη ) (black solid line) for q = 1 / (cid:60) ( I M ( k sin( iη )) (red dashed line) (b) Difference ce − (0 . (cid:60) ( I M ) Fig. 15 . In the limit as q → and η → ∞ keeping ξ = 2 √ q sinh η constant, a Mathieufunction ce ( q, z ) rapidly approaches a Bessel function (cid:60) I M (2 √ q sin( z )) . In this figure, q = 1 / and M ≈ . i . Here z = iη is purely imaginary. We have ignored the O ( k ) term of the outersolution. Even so, the match seems good. g >
3, cos ( g − α
384 (( g − g −
2) ( g −
3) + (cid:0) g − g + 7 (cid:1) cos ( g − α
128 ( g + 1) ( g −
2) ( g − − (cid:0) g + 4 g + 7 (cid:1) cos ( g + 2) α
128 ( g + 2) ( g −
1) ( g + 1) − cos ( g + 6) α
384 ( g + 3) ( g + 2) ( g + 1)(C.3)Again this is in agreement with Mathieu, although he wrote 128 as 2 and 384 as2 ·
3. It is at this term that the difference between Mathieu’s form of the series and amodern series expansion of ce g ( q, α ) is most noticeable, because his work eliminatedthe cos gα term at this order (and altered the other coefficients). The coefficient of q is, provided g > g − α g −
1) ( g −
2) ( g −
3) ( g −
4) + (cid:0) g − g + 10 (cid:1) cos ( g − α
768 ( g −
2) ( g −
3) ( g + 1) ( g − + (cid:0) g + 5 g + 10 (cid:1) cos ( g + 4) α
768 ( g −
1) ( g + 3) ( g + 2) ( g + 1) + cos ( g + 8) α g + 4) ( g + 3) ( g + 2) ( g + 1) . (C.4)Mathieu wrote 6144 as 2 · ·
3; he also wrote g + 7 g + 20 g + 20 as thenumerator for the cos( g + 4) α term (and similar for the cos( g − α and had an extrafactor g + 2 in the denominator. Since g + 7 g + 20 g + 20 = ( g + 2) (cid:0) g + 5 g + 10 (cid:1) ,we see that his result was again correct though not in simplest form. The coefficientof q is (and instead of editing it for elegance, this time we leave it as automatically4 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR generated by Maple), provided g > α g − α )(122880 g − g −
1) ( g −
2) ( g −
3) ( g − (cid:0) g − g + 13 (cid:1) cos ( α g − α )(8192 g − g −
3) ( g −
4) ( g + 1) ( g − + (cid:0) g − g + 8 g − g + 47 g + 13 g + 232 (cid:1) cos ( α g − α )(3072 g − g −
3) ( g + 2) ( g + 1) ( g − − (cid:0) g + 5 g + 8 g + 8 g + 47 g − g + 232 (cid:1) cos ( α g + 2 α )(3072 g − g + 3) ( g + 2) ( g − ( g + 1) − (cid:0) g + 6 g + 13 (cid:1) cos ( α g + 6 α )(8192 g − g + 4) ( g + 3) ( g + 2) ( g + 1) − cos ( α g + 10 α )(122880 g + 614400) ( g + 4) ( g + 3) ( g + 2) ( g + 1)(C.5)Mathieu wrote 122880 as 2 · · ·
5, and similarly other large numbers in factoredform. He had g + 11 g + 49 g + 101 g + 78 = ( g + 3) ( g + 2) (cid:0) g + 6 g + 13 (cid:1) in thenumerator of the cos( g + 6) α term and an extra ( g + 3)( g + 2) in the denominator.Similarly for the cos( g − α term. For the cos( g − α term he had g + 7 g + 18 g +24 g + 63 g + 81 g + 206 g + 464 which is g + 2 times the numerator printed above;of course he had an extra factor g + 2 in the denominator to cancel it.Finally, he has an extra factor ( g −
4) in the denominator of the q ( h ) termof the eigenvalue and an apparently incorrect numerator, 9 g + 22 g − g − g + 22 g − g −
116 = ( g − g + 58 g + 29). Actually, on theline above Mathieu’s final form for R , the power is more clearly a 6: this isn’t even atypo, 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 pertur-bation series to ce g ( h , α ) /F up to and including terms of order h , on the under-standing that the factor F was chosen to make all coefficients of cos gα equal to zeroapart from the first one. Special series . The generic series is good only for “large enough” g . For specificsmall g , and indeed for any fixed g if one wants to compute enough terms, specialcomputations have to be done. We here give the results that Mathieu attempted, andcomment on any errors in his paper. ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE g = 2. We get (again removing the coefficients of cos 2 α as he did)ce ( q, α ) F = cos (2 α ) + (cid:18) − cos (4 α )12 (cid:19) q + cos (6 α )384 q + (cid:18) − −
43 cos (4 α )13824 − cos (8 α )23040 (cid:19) q + (cid:18)
293 cos (6 α )2211840 + cos (10 α )2211840 (cid:19) q + (cid:18) α )79626240 −
167 cos (8 α )66355200 − cos (12 α )309657600 (cid:19) q + (cid:18) − α )12740198400 + 629 cos (10 α )22295347200 + cos (14 α )59454259200 (cid:19) q + O (cid:0) q (cid:1) (C.6)where the factor F (not computed by Mathieu, but needed by us to compare a modernseries to his results) is(C.7) F = 1 − q
288 + 51191 q − q . Notice first that Ince was correct: there are values of q for which that factor F is zero,and therefore this normalization is not universally possible. Here we also see apparentarithmetic errors: Mathieu has 287 and not 293 as we have in the O ( q ) term. Hehas 21059 where we have 21041 in the O ( q ) term, and 41 instead of 167. The othertwo terms at that order are correct. He did not report the O ( q ) term.For the eigenvalue at g = 2 Mathieu reports the correct expansion,(C.8) a = 4 + 512 q − q + 100240179626240 q + O (cid:0) q (cid:1) , except he has 1002419 instead of 1002401.For g = 4, Mathieu reports everything correctly up until the constant (cos 0 α )term of the O ( q ) term: he gets − / / O ( q ) term we get(C.9) −
53 cos (2 α )124416000 − α )2419200000 −
53 cos (10 α )1032192000 − cos (14 α )1857945600while Mathieu gets(C.10) − cos 14 α − cos 10 α − cos 6 α − cos 2 α which has evident discrepancies at the cos 2 α and cos 14 α terms but is otherwisecorrect.For the eigenvalue, we get(C.11) a = 16 + 130 q + 433864000 q − q + O (cid:0) q (cid:1) which is nearly the same as Mathieu’s,(C.12) R = 16 + h + h − h + · · · except he erroneously reports 189983 / q . This is anirreducible fraction and not equal to the correct coefficient.6 C. BRIMACOMBE, R. M. CORLESS, AND M. ZAMIR
For g = 1ce ( q, α ) F = cos ( α ) − cos (3 α )8 q + (cid:18) − cos (3 α )64 + cos (5 α )192 (cid:19) q + (cid:18) cos (5 α )1152 − cos (3 α )1536 − cos (7 α )9216 (cid:19) q + (cid:18) − cos (7 α )49152 + 11 cos (3 α )36864 + cos (5 α )24576 + cos (9 α )737280 (cid:19) q + (cid:18) cos (9 α )3686400 − α )393216 + 49 cos (3 α )589824 − cos (7 α )983040 − cos (11 α )88473600 (cid:19) q + (cid:18) − cos (11 α )424673280 + 17 cos (7 α )39321600 + 55 cos (3 α )9437184 −
719 cos (5 α )141557760+ cos (9 α )70778880 + cos (13 α )14863564800 (cid:19) q + O (cid:0) q (cid:1) . (C.13)The correction factor is(C.14) F = 1 − q − q − q + 1211769472 q + 8105339738624 q + · · · The eigenvalue is(C.15) a = 1 + q − q − q − q + 1136864 q + 49589824 q + O (cid:0) q (cid:1) . Mathieu getsce ( h , α ) F = cos α − h cos 3 α + h (cid:0) − cos 5 α − cos 3 α (cid:1) − h (cid:0) cos 7 α − cos 5 α + cos 3 α (cid:1) + h (cid:0) cos 9 α − cos 7 α + cos 5 α + cos 3 α (cid:1) + · · · (C.16)and for the eigenvalue gets(C.17) R = 1 + h − h − h − h + h + · · · . All terms are in complete agreement with our results.For g = 3 we havece ( q, α ) F = cos (3 α ) + (cid:18) cos ( α )8 − cos (5 α )16 (cid:19) q + (cid:18) cos ( α )64 + cos (7 α )640 (cid:19) q + (cid:18) cos ( α )1024 − α )20480 − cos (9 α )46080 (cid:19) q + (cid:18) − cos (5 α )16384 − cos ( α )4096 + 17 cos (7 α )1474560 + cos (11 α )5160960 (cid:19) q + O (cid:0) q (cid:1) (C.18)where(C.19) F = 1 − q − q − q q O ( q ) . ATHIEU FUNCTIONS: HISTORICAL PERSPECTIVE a = 9 + 116 q + 164 q + 1320480 q − q − q + O (cid:0) q (cid:1) . Mathieu gets P = cos 3 α + h (cid:0) − cos 5 α + cos α (cid:1) + h (cid:0) cos 7 α + cos α (cid:1) + h (cid:0) − cos 9 α − cos 5 α + cos α (cid:1) + h (cid:0) · · · cos 11 α − · · cos 7 α − cos 5 α − cos α (cid:1) + · · · ; R = 9 + h + h + h − h + · · · . Mathieu’s eigenfunction seems already wrong at q ( h ) in several coefficients: 1 / /
8, 1 /
768 instead of 1 / O ( q )coefficient.Similarly, Mathieu’s series for se ( q, α ) /F is correct in every term, whilst his seriesfor se ( q, α ) /F is wrong already at O ( q4