Pure Tone Modes
for a $5\mathrel{\mathop{\ordinarycolon}}3$ Elliptic Drum

In 1868 Émile Mathieu published a memoir showing how to find the modes, often called standing modes, 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 that memoir in detail; indeed we asked my colleague Robert H. C. Moir to translate the memoir from French in order to make that job easier. We learned that Mathieu had solved the linear model PDE 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, including a 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 that survive 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 and location of nodal lines (hyperbolas and ellipses, actually) make more sense. In this note I show some of the pictures I computed for myself using the (somewhat quixotic) software that I wrote for the Mathieu functions and for solving the Mathieu equation.

We parameterize the elliptic geometry as Mathieu did (he credited the change of coordinates to Lamé who called the result “thermometric” coordinates) through what is now called the confocal elliptic coordinates $x=c\cosh\beta\cos\alpha$, $y=c\sinh\beta\sin\alpha$, or equivalently $x+iy=c\cos(\alpha-i\beta)$. The real parameter $c$ gives the distance $2c$ between the foci of the ellipse. These coordinates are singular in the case $c=0$ when the ellipse becomes a circle, which is a bit of a headache: the proper analysis of the reduction to the circular case requires a double limit process (see Appendix B of [1] for a sketch of such an analysis). Mathieu starts with a wave equation, containing a material parameter $m^{2}$ representing the ratio of stress to density in the membrane, and then by separation of variables $w=P(\alpha)Q(\beta)\sin 2\lambda mt$ with the frequency of oscillation parameter $\lambda m$ dealing with the physical parameter; here $P$ is a solution of the so-called Mathieu differential equation and $Q$ is a solution of the so-called modified Mathieu equation. These are also called the angular and radial Mathieu equations.

The Mathieu differential equation is

and the modified Mathieu differential equation is

where $q$ is a parameter that depends on the focal parameter $c$ and the frequency of oscillation $2\lambda m$; the parameter $a$ is an eigenvalue of the differential equation (1). The eigenvalue $a$ depends on the value of $q=\lambda^{2}c^{2}$, and must be chosen to make the solution $y(\alpha)$ periodic with period either $\pi$ or $2\pi$. Mathieu also divided solutions up into even and odd categories, analogous to cosine and to sine. As in [1] we follow Mathieu and denote the Mathieu functions by $\mathrm{ce}_{g}(q;\alpha)$ for the even solutions and $\mathrm{se}_{g}(q;\alpha)$ for the odd, and the modified Mathieu functions by $\mathrm{Ce}_{g}(q;\beta)$ and $\mathrm{Se}_{g}(q;\beta)$. Here $g$ is an integer (Mathieu used the letter $g$ in this way, and it might seem unusual to modern eyes that are now used to the Fortran I-N convention of single-letter integer variables being $i$, $j$, $k$, $\ell$, $m$, or $n$). We have $g\geq 0$ for the even solutions and $g\geq 1$ for the odd. We normalize the eigenfunctions as follows:

We normalize the modified Mathieu functions analogously, so that $\mathrm{Ce}_{g}(q;\beta)=\mathrm{ce}_{g}(q;i\beta)$ and $\mathrm{Se}_{g}(q;\beta)=-i\mathrm{se}_{g}(q;i\beta)$. This normalization is used in some places in the literature, but most people use a different one involving the square integral over the period. See [1] for a discussion of why we chose this one.

The frequency of vibration associated with a particular standing mode depends also on the material properties of the membrane as encapsulated in the parameter $m^{2}$, the ratio of tension to density in the (homogeneous) membrane. If the membrane is oscillating like $\sin 2\lambda mt$ then via $q=\lambda^{2}c^{2}$ we have

Some authors use $h^{2}$ in place of $q$.

It is the so-called modified Mathieu functions $\mathrm{Ce}_{g}(q;\beta)=\mathrm{ce}_{g}(q;i\beta)$ and $\mathrm{Se}_{g}(q;\beta)=-i\mathrm{se}_{g}(q;i\beta)$ which determine, through the boundary conditions (the membrane is fixed at the edges, so $y=0$ there), the value of $q$. Mathieu established that every natural mode of the membrane is of the form

or

The dependence on the frequency of oscillation can be made clearer by writing a mode in full as

That is, multiple values of $q$ are possible for the same PDE with fixed parameter $m$ because different modes of vibration can be excited at the same time. Therefore a general solution composed of a sum of different excited modes would be a sum over modes with different values of $\lambda$ and therefore of $q$.

If we have a fixed elliptical drum in mind, then it will have a given aspect ratio or equivalently eccentricity. Say, for the sake of argument, that the aspect ratio is $5\mathrel{\mathop{\ordinarycolon}}3$. Then we have (putting the major axis on the $x$-axis)

defining the outer boundary of the ellipse (this also defines the units we use): this is quickly solved¹¹1$\tanh\beta_{0}=3/5$ so $\beta_{0}=\ln 2$ and $5^{2}-3^{2}=c^{2}=16$. to get $c=4$ and $\beta_{0}=\ln 2$. Part of the attraction of choosing this aspect ratio for demonstration is 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 solve either

or

for the parameter $q$. These are (of course) nonlinear equations, and for real values of $\lambda$ and $c$ turn out always to have an infinite number of solutions, for each value of the index $g$. See figure 1 where we show a representative pair of modified Mathieu equations and see figure 2 where we fix $\beta=\beta_{0}=\ln 2$ and vary $q$, exhibiting an oscillatory function which apparently has an infinite number of zeros. In practice we are only interested in the first few; the very large values of $q$ giving distant zeros will give impractically high-frequency standing modes anyway.

We wish to find, for a given nonnegative integer $g$ and choice of “even” or “odd”, the values of $q$ for which $Q(q;\ln 2)=0$, where $Q(q;\beta)$ is either $\mathrm{Ce}_{g}(q;\beta)$ (even case) or $\mathrm{Se}_{g}(q;\beta)$ (odd case) for the appropriate $g$. Either way this is a nonlinear equation in the real parameter $q$ that has an infinite number of solutions, as evidenced by figures 2(a)–2(d). In those figures, it is important that we chose to normalize by making $\mathrm{Ce}_{g}(q;0)=1$ while its derivative with respect to $\beta$ was zero there; if instead we had chosen to use the integral normalization then these values seem to decay root-exponentially, i.e. like $\exp(-2\sqrt{q})$.

We also wish to find $q$ so that $\ln 2$ is the first zero of the function $Q(q;\beta)$, and then another value of $q$ so that $\ln 2$ 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_{g,k}^{e}$ for the for the value of $q$ that makes $\ln 2$ the $k$th zero of $\mathrm{Ce}_{g}(q;\beta)$ (the even case), for $k=1$, $2$, $3$, $\ldots$. Similarly we will write $q_{g,k}^{o}$ for the value of $q$ that makes $\ln 2$ the $k$th zero of $\mathrm{Se}_{g}(q;\beta)$ (note that $0=\mathrm{Se}_{g}(q;0)$ for all $q$ and we can say this is the $0$th zero, if we like: then here as with $\mathrm{Ce}$ the index $k$ starts at $1$).

By trying to plot $Q(q;\beta)$ as a function of $\beta$ for various $q$ (we need first to find the appropriate eigenvalue $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, monotonic. Nonetheless, by simply plotting the function for a few values of $q$ one can quickly get an initial estimate $q_{0}$ of the value of $q$ for which $Q(q;\ln 2)=0$.

Once one has a decent initial estimate $q_{0}$, one would like to use Newton’s method to find the value of $q$ more precisely. This requires derivatives with respect to $q$, which ordinarily would require a Fréchet derivative and the use of a Green’s function (which works perfectly well) but I can report with pleasure that the simple Squire-Trapp formula [9] for numerical differentiation works beautifully here (see [5, p. 479] for a discussion of the excellent numerical properties of this formula). The derivation of the formula is simple: if $f$ is analytic and $x$ and $h$ are real, then

Therefore, taking the imaginary part gives

So all we need to do to compute $f^{\prime}(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 be computed 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 square root of machine epsilon (or whatever tolerance we are computing $Q(q;\beta)$ to) and compute instead with $q+ih$; we use the real part of the answer for $Q(q;\ln 2)$ and the imaginary part (divided by $h$) as $dQ/dq$. Since $h^{2}$ will then be smaller than machine epsilon, we can ignore the effects of the perturbation on the real part. This is a kind of finite difference formula, but one that for analytic $f$ suffers no instability as $h\to 0$. Really, it’s ridiculously effective.

We may then use Newton’s method

where the prime now denotes differentiation with respect to $q$. All of the results in tables 1 and 2 were computed in this way. I only report the values of $q$ to a few figures; more can be computed on demand.

The eigenvalues $a_{g}(q_{n})$ and $b_{g}(q_{n})$ were computed with an even faster iteration, the Inverse Cubic Iteration described in [4], essentially just for fun. The results are not tabulated here because they can be recomputed easily, given $q$.

Once one has the values of $q$ which are needed, it is straightforward to compute the mode shapes. Somewhat surprisingly, it is awkward to plot them (either in Maple or in Matlab) because the surface is defined parametrically; most contour plotting software wants a regular $x$-$y$ grid to work from. In Maple, this can be worked around by using plots[surfdata] which will quickly do a contour plot of irregular data as a decoration of its 3D plots, which can be looked at from the top down (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.

Instead, I inverted the coordinate transformation, so $\alpha-i\beta=\arccos((x+iy)/4)$, and simply evaluated my numerical solutions of the Mathieu equation (represented as strings of blends [6]) at the resulting $\alpha$ and $\beta$. 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 $400$ by $400$ grid. In contrast, the plots[surfdata] approach only takes about three seconds in total for a similar high-resolution graph. But in the meantime, the code was good enough to allow the production of the contour plots in this paper.

In figure 3 we see the contours corresponding to an even solution with $g=3$, where the value of $q$ was chosen so that $\ln 2$ was the first zero of $\mathrm{Ce}_{3}(q;\beta)$. Figures 4–14 plot only the contours of the first quadrant; by symmetry the other three can be deduced. We see confocal elliptic nodal lines, and hyperbolic nodal lines, in several figures; by comparing all the figures we may understand what Mathieu was talking about when he discussed the natural modes of vibration of an elliptic drum.

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 the fourth zero of all modes plotted in figures 4–14. In a “bouncing ball” mode, most of the vibration is confined to a central channel.

None of the pictures in this paper would have surprised Mathieu. He predicted the hyperbolic and elliptic nodal lines, and the exact numbers of the hyperbolic lines by an ingenious sign variation argument 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, although I did not read his description of the necessary boundary conditions all that closely (these boundary conditions are covered briefly in Chapter 28 of the Digital Library of Mathematical Functions https://dlmf.nist.gov/28 and more thoroughly in [7]).

From the $q$ data in tables 1 and 2 we may deduce the frequencies $\lambda$ of each of these modes. Writing the general solution as a sum, we have

where the relation $\lambda_{g,k}^{e,o}=\sqrt{q_{g,k}^{e,o}}/c$ ties each frequency to its mode; we have only tabulated the lowest frequency modes in this paper. The unknown coefficients $A_{g,k}$ and $B_{g,k}$ must be determined by the initial displacement impulse.