PPure Tone Modesfor a : Elliptic Drum
Robert M. CorlessScientific Director, Ontario Research Centre for Computer Algebramember, The Rotman Institute of PhilosophyThe University of Western OntarioAugust 18, 2020
In 1868 Émile Mathieu published a memoir showing how to find the modes, often called standingmodes , that generate pure tones in the vibration of an elliptic drum of arbitrary aspect ratio [8].In preparing [1], together with my co-authors Chris Brimacombe and Mair Zamir, I read thatmemoir in detail; indeed we asked my colleague Robert H. C. Moir to translate the memoir fromFrench in order to make that job easier. We learned that Mathieu had solved the linear modelPDE problem completely, and had given a detailed discussion of the possible nodal lines—that is,places on the membrane that would not move if the drum was oscillating in that mode, includinga complete zero-counting discussion using Sturm theory, which he extended to the periodic case.But he published no pictures; perhaps he had sketched them for himself, but we know of none thatsurvive from that time. Later authors, of course, have supplied such pictures: see for instance [7]or [2] or [3]. Once one has seen such pictures, Mathieu’s verbal descriptions of the number andlocation of nodal lines (hyperbolas and ellipses, actually) make more sense. In this note I showsome of the pictures I computed for myself using the (somewhat quixotic) software that I wrotefor the Mathieu functions and for solving the Mathieu equation.We parameterize the elliptic geometry as Mathieu did (he credited the change of coordinates toLamé who called the result “thermometric” coordinates) through what is now called the confocalelliptic coordinates x = c cosh β cos α , y = c sinh β sin α , or equivalently x + iy = c cos ( α − iβ ) . Thereal parameter c gives the distance between the foci of the ellipse. These coordinates are singular in the case c = when the ellipse becomes a circle, which is a bit of a headache: the proper analysisof the reduction to the circular case requires a double limit process (see Appendix B of [1] for asketch of such an analysis). Mathieu starts with a wave equation, containing a material parameter m representing the ratio of stress to density in the membrane, and then by separation of variables w = P ( α ) Q ( β ) sin with the frequency of oscillation parameter λm dealing with the physicalparameter; here P is a solution of the so-called Mathieu differential equation and Q is a solutionof the so-called modified Mathieu equation. These are also called the angular and radial Mathieuequations.The Mathieu differential equation is y (cid:48)(cid:48) ( α ) + ( a − cos ) y ( α ) = (1)1 a r X i v : . [ m a t h . HO ] A ug nd the modified Mathieu differential equation is y (cid:48)(cid:48) ( β ) − ( a − cosh ) y ( β ) = (2)where q is a parameter that depends on the focal parameter c and the frequency of oscillation ;the parameter a is an eigenvalue of the differential equation (1). The eigenvalue a depends on thevalue of q = λ c , and must be chosen to make the solution y ( α ) periodic with period either π or . Mathieu also divided solutions up into even and odd categories, analogous to cosine andto sine. As in [1] we follow Mathieu and denote the Mathieu functions by ce g ( q ; α ) for the evensolutions and se g ( q ; α ) for the odd, and the modified Mathieu functions by Ce g ( q ; β ) and Se g ( q ; β ) .Here g is an integer (Mathieu used the letter g in this way, and it might seem unusual to moderneyes that are now used to the Fortran
I-N convention of single-letter integer variables being i , j , k , (cid:96) , m , or n ). We have g ≥ for the even solutions and g ≥ for the odd. We normalize theeigenfunctions as follows: ce g ( q ; ) = ce g ( q ; α ) (cid:12)(cid:12)(cid:12)(cid:12) α = = se g ( q ; ) = se g ( q ; α ) (cid:12)(cid:12)(cid:12)(cid:12) α = = (3)We normalize the modified Mathieu functions analogously, so that Ce g ( q ; β ) = ce g ( q ; iβ ) and Se g ( q ; β ) = − i se g ( q ; iβ ) . This normalization is used in some places in the literature, but mostpeople use a different one involving the square integral over the period. See [1] for a discussion ofwhy we chose this one.The frequency of vibration associated with a particular standing mode depends also on thematerial properties of the membrane as encapsulated in the parameter m , the ratio of tensionto density in the (homogeneous) membrane. If the membrane is oscillating like sin then via q = λ c we have λ = √ qc = hc . (4)Some authors use h in place of q .It is the so-called modified Mathieu functions Ce g ( q ; β ) = ce g ( q ; iβ ) and Se g ( q ; β ) = − i se g ( q ; iβ ) which determine, through the boundary conditions (the membrane is fixed at the edges, so y = there), the value of q . Mathieu established that every natural mode of the membrane is of the form Ce g ( q ; β ) ce g ( q ; α ) (5)or Se g ( q ; β ) se g ( q ; α ) . (6)The dependence on the frequency of oscillation can be made clearer by writing a mode in full as P ( q ; α ) Q ( q ; β ) sin = P ( λ c ; α ) Q ( λ c ; β ) sin (7)That is, multiple values of q are possible for the same PDE with fixed parameter m because differentmodes of vibration can be excited at the same time. Therefore a general solution composed of a2um of different excited modes would be a sum over modes with different values of λ and thereforeof q .If we have a fixed elliptical drum in mind, then it will have a given aspect ratio or equivalentlyeccentricity. Say, for the sake of argument, that the aspect ratio is : . Then we have (puttingthe major axis on the x -axis) = c cosh β = c sinh β defining the outer boundary of the ellipse (this also defines the units we use): this is quickly solved to get c = and β = ln . Part of the attraction of choosing this aspect ratio for demonstrationis that the numbers are so simple.To find the permissible values of q that leave the edge of the drum fixed, then, we have to solveeither Ce g ( q ; ln ) = (8)or Se g ( q ; ln ) = (9)for the parameter q . These are (of course) nonlinear equations, and for real values of λ and c turnout always to have an infinite number of solutions, for each value of the index g . See figure 1where we show a representative pair of modified Mathieu equations and see figure 2 where wefix β = β = ln and vary q , exhibiting an oscillatory function which apparently has an infinitenumber of zeros. In practice we are only interested in the first few; the very large values of q givingdistant zeros will give impractically high-frequency standing modes anyway. q numerically We wish to find, for a given nonnegative integer g and choice of “even” or “odd”, the valuesof q for which Q ( q ; ln ) = , where Q ( q ; β ) is either Ce g ( q ; β ) (even case) or Se g ( q ; β ) (odd case)for the appropriate g . Either way this is a nonlinear equation in the real parameter q that hasan infinite number of solutions, as evidenced by figures 2a–2d. In those figures, it is importantthat we chose to normalize by making Ce g ( q ; ) = while its derivative with respect to β was zerothere; if instead we had chosen to use the integral normalization then these values seem to decayroot-exponentially, i.e. like exp (− √ q ) .We also wish to find q so that ln is the first zero of the function Q ( q ; β ) , and then anothervalue of q so that ln is the second zero, and so on. Each choice of g and of “even” versus “odd’will give a countably infinite number of values of q . There is no standard notation for these zeros.For now we will write q eg,k for the for the value of q that makes ln the k th zero of Ce g ( q ; β ) (theeven case), for k = , , , . . . . Similarly we will write q og,k for the value of q that makes ln the k th zero of Se g ( q ; β ) (note that = Se g ( q ; ) for all q and we can say this is the th zero, if welike: then here as with Ce the index k starts at ).By trying to plot Q ( q ; β ) as a function of β for various q (we need first to find the appropriateeigenvalue a , and how to do this is described in [1] for example) we see that in general increasing q brings the first zero closer to the origin. However, the motion of the zeros is not, in general, tanh β = so β = ln and − = c = . a) Ce ( ; β ) = ce ( ; iβ ) (b) Se ( ; β ) = (cid:61) ( se ( ; iβ )) Figure 1: (Left) A graph of Ce ( q ; β ) when q = . 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 locationsof the foci at ± c ). By construction, the first such zero occurs at β = , giving an aspect ratioof : . (Right) A graph of Se ( q ; β ) = (cid:61) ( se ( q ; iβ )) when q = . Again the first zero allowsa pure tone for an ellipse of aspect ratio : , by construction. (a) Ce ( q ; ln ) (b) √ q Se ( q ; ln ) (c) Ce ( q ; ln ) (d) √ q Se ( q ; ln ) Figure 2: Modified Mathieu functions for β = ln as a function of q . From the graphs, we wouldfind it natural that there are an infinite number of values of q for which Ce g ( q ; ln ) = andsimilarly Se g ( q ; ln ) = . With the normalization used in this paper, it seems that Ce g ( q ; ln ) eventually oscillates between − and , while Se g ( q ; ln ) decays like ( √ q ) .4onotonic. Nonetheless, by simply plotting the function for a few values of q one can quickly getan initial estimate q of the value of q for which Q ( q ; ln ) = .Once one has a decent initial estimate q , one would like to use Newton’s method to find thevalue of q more precisely. This requires derivatives with respect to q , which ordinarily wouldrequire a Fréchet derivative and the use of a Green’s function (which works perfectly well) butI can report with pleasure that the simple Squire-Trapp formula [9] for numerical differentiationworks beautifully here (see [5, p. 479] for a discussion of the excellent numerical properties of thisformula). The derivation of the formula is simple: if f is analytic and x and h are real, then f ( x + ih ) = f ( x ) + ihf (cid:48) ( x ) − h f (cid:48)(cid:48) ( x ) − i h f (cid:48)(cid:48)(cid:48) ( x ) + O ( h ) . (10)Therefore, taking the imaginary part gives (cid:61) ( f ( z + ih )) h = f (cid:48) ( x ) − h f (cid:48)(cid:48)(cid:48) ( x ) + O ( h ) , = f (cid:48) ( x ) + O ( h ) . (11)So all we need to do to compute f (cid:48) ( x ) is to compute f near to x but with a small purely imaginary perturbation. The formula does require f to be analytic because it effectively uses the Cauchy-Riemann equations (here Q is in fact entire) and requires f to be real when x is real and to becomputed to high relative accuracy in both the real and imaginary part. This holds in our case.That is, instead of computing with our real value of q , we choose an h smaller than the squareroot of machine epsilon (or whatever tolerance we are computing Q ( q ; β ) to) and compute insteadwith q + ih ; we use the real part of the answer for Q ( q ; ln ) and the imaginary part (divided by h ) as dQ/dq . Since h will then be smaller than machine epsilon, we can ignore the effects of theperturbation on the real part. This is a kind of finite difference formula, but one that for analytic f suffers no instability as h → . Really, it’s ridiculously effective.We may then use Newton’s method q n + = q n − Q ( q n ; ln ) Q (cid:48) ( q n ; ln ) where the prime now denotes differentiation with respect to q . All of the results in tables 1 and 2were computed in this way. I only report the values of q to a few figures; more can be computedon demand.The eigenvalues a g ( q n ) and b g ( q n ) were computed with an even faster iteration, the InverseCubic Iteration described in [4], essentially just for fun. The results are not tabulated here becausethey can be recomputed easily, given q . Once one has the values of q which are needed, it is straightforward to compute the modeshapes. Somewhat surprisingly, it is awkward to plot them (either in Maple or in Matlab) becausethe surface is defined parametrically; most contour plotting software wants a regular x - y grid towork from. In Maple, this can be worked around by using plots[surfdata] which will quickly doa contour plot of irregular data as a decoration of its 3D plots, which can be looked at from the topdown (use the keyword orientation=[-90,0] in the plots[surfdata] command). Unfortunately,this is somewhat “too fancy” and does not produce a nice-looking simple 2D contour plot.5able 1: Values of the parameter q for which Ce g ( q, ln ) = . The i th table column gives thevalues of q for which the i th zero of Ce g ( q ; β ) is β = ln . The corresponding eigenvalues are a g ( q ) . g first second third fourth0 1.7353 11.356 29.795 57.0111 3.3522 14.627 34.844 63.8482 5.6530 18.486 40.457 71.2413 8.6576 22.968 46.658 79.2014 12.368 28.100 53.463 87.7545 16.779 33.913 60.891 96.903Table 2: Values of the parameter q for which Se g ( q, ln ) = . The i th table column gives thevalues of q for which the i th zero of Se g ( q ; β ) is β = ln . The corresponding eigenvalues are b g ( q ) . g first second third fourth1 5.4300 19.478 42.306 73.9022 7.8189 23.636 48.248 81.6263 10.803 28.363 54.775 89.9144 14.406 33.696 61.800 98.7625 18.637 39.642 69.499 108.17Instead, I inverted the coordinate transformation, so α − iβ = arccos (( x + iy ) /4 ) , and simplyevaluated my numerical solutions of the Mathieu equation (represented as strings of blends [6]) atthe resulting α and β . At low resolution the contour plots take only a few seconds for each one,but take up to about a minute for high resolutions such as a by grid. In contrast, the plots[surfdata] approach only takes about three seconds in total for a similar high-resolutiongraph. But in the meantime, the code was good enough to allow the production of the contourplots in this paper.In figure 3 we see the contours corresponding to an even solution with g = , where the valueof q was chosen so that ln was the first zero of Ce ( q ; β ) . Figures 4–14 plot only the contoursof the first quadrant; by symmetry the other three can be deduced. We see confocal elliptic nodallines, and hyperbolic nodal lines, in several figures; by comparing all the figures we may understandwhat Mathieu was talking about when he discussed the natural modes of vibration of an ellipticdrum.The paper [3] intriguingly describes “whispering gallery” modes, two of which I plot in figure 15),and “bouncing ball” modes, which seem to be exhibited in the (d) figures corresponding to thefourth zero of all modes plotted in figures 4–14. In a “bouncing ball” mode, most of the vibrationis confined to a central channel. None of the pictures in this paper would have surprised Mathieu. He predicted the hyperbolicand elliptic nodal lines, and the exact numbers of the hyperbolic lines by an ingenious sign variationargument with Sturm sequences. He also noted that his methods could be used to solve the“confocal elliptic lake” problem—that is, a drum with a confocal hole taken out from the middle,6lthough I did not read his description of the necessary boundary conditions all that closely (theseboundary conditions are covered briefly in Chapter 28 of the Digital Library of MathematicalFunctions https://dlmf.nist.gov/28 and more thoroughly in [7]).From the q data in tables 1 and 2 we may deduce the frequencies λ of each of these modes.Writing the general solution as a sum, we have u ( x, y, t ) = (cid:88) g ≥ (cid:88) k ≥ A g,k Ce g ( q eg,k ; β ) ce g ( q eg,k ; α ) sin eg,k mt + (cid:88) g ≥ (cid:88) k ≥ B g,k Se g ( q og,k ; β ) se g ( q og,k ; α ) sin og,k mt , (12)where the relation λ e,og,k = (cid:113) q e,og,k /c ties each frequency to its mode; we have only tabulated thelowest frequency modes in this paper. The unknown coefficients A g,k and B g,k must be determinedby the initial displacement impulse. References [1] Chris Brimacombe, Robert M Corless, and Mair Zamir. Computation and applications ofMathieu functions: A historical perspective. arXiv preprint arXiv:2008.01812 , 2020.[2] L Chaos-Cador and E Ley-Koo. Mathieu functions revisited: matrix evaluation and generatingfunctions.
Revista mexicana de física , 48(1):67–75, 2002.[3] Goong Chen, Philip J Morris, and Jianxin Zhou. Visualization of special eigenmode shapes ofa vibrating elliptical membrane.
SIAM review , 36(3):453–469, 1994.[4] Robert M Corless. Inverse cubic iteration. arXiv preprint arXiv:2007.06571 , 2020.[5] Robert M Corless and Nicolas Fillion.
A graduate introduction to numerical methods .Springer-Verlag, 2013.[6] Robert M Corless and Erik Postma. Blends in Maple. arXiv preprint arXiv:2007.05041 , 2020.[7] Julio C Gutiérrez-Vega, RM Rodrıguez-Dagnino, MA Meneses-Nava, and S Chávez-Cerda.Mathieu functions, a visual approach.
American Journal of Physics , 71(3):233–242, 2003.[8] Émile Mathieu. Mémoire sur le mouvement vibratoire d’une membrane de forme elliptique.
Journal de mathématiques pures et appliquées , 13:137–203, 1868.[9] William Squire and George Trapp. Using complex variables to estimate derivatives of realfunctions.
SIAM Review , 40(1):pp. 110–112, 1998.7 a) Contour Plot (b) plots[surfdata]
Figure 3: A typical pure tone mode of a : ellipse. Both graphs show Ce ( q ; β ) ce ( q ; α ) where q = , x = cosh β cos α , y = sinh β sin α . The parameters − π < α < π and ≤ β ≤ ln .At β = ln , cosh β = and sinh β = , delimiting the elliptic boundary. The left figure is a fullversion of figure 9a. The right shows a colour 3D plot of the same. (a) q = (b) q = (c) q = (d) q = Figure 4: Contours of Ce ( q ; β ) ce ( q ; α ) in the first quadrant. (a) q = (b) q = (c) q = (d) q = Figure 5: Contours of Ce ( q ; β ) ce ( q ; α ) in the first quadrant.8 a) q = (b) q = (c) q = (d) q = Figure 6: Contours of Se ( q ; β ) se ( q ; α ) in the first quadrant. (a) q = (b) q = (c) q = (d) q = Figure 7: Contours of Ce ( q ; β ) ce ( q ; α ) in the first quadrant. (a) q = (b) q = (c) q = (d) q = Figure 8: Contours of Se ( q ; β ) se ( q ; α ) in the first quadrant.9 a) q = (b) q = (c) q = (d) q = Figure 9: Contours of Ce ( q ; β ) ce ( q ; α ) in the first quadrant. (a) q = (b) q = (c) q = (d) q = Figure 10: Contours of Se ( q ; β ) se ( q ; α ) in the first quadrant. (a) q = (b) q = (c) q = (d) q = Figure 11: Contours of Ce ( q ; β ) ce ( q ; α ) in the first quadrant.10 a) q = (b) q = (c) q = (d) q = Figure 12: Contours of Se ( q ; β ) se ( q ; α ) in the first quadrant. (a) q = (b) q = (c) q = (d) q = Figure 13: Contours of Ce ( q ; β ) ce ( q ; α ) in the first quadrant. (a) q = (b) q = (c) q = (d) q = Figure 14: Contours of Se ( q ; β ) se ( q ; α ) in the first quadrant.11 a) g = (b) g = (c) g = (d) g = Figure 15: Two “whispering gallery” modes, after [3]. These are contours of Ce g ( q ; β ) ce g ( q ; α ) for g = and q ≈ (left pair) and for g = with q ≈ (right pair). The name is becausemost of the vibration occurs near the boundary: the idea is that if these were rooms with ellipticshape, a person at one vertex could whisper at the appropriate frequency and be heard by a personat the other vertex, while another person in the middle could not hear what was being said. Thecolour plots were produced using plots[surfdata]plots[surfdata]