Abstract
Our main aim is the accurate computation of a large number of specified eigenvalues and eigenvectors of Mathieu’s system as a multiparameter eigenvalue problem (MEP). The reduced wave equation, for small deflections, is solved directly without approximations introduced by the classical Mathieu functions. We show how for moderate values of the cutoff collocation parameter the QR algorithm and the Arnoldi method may be applied successfully, while for larger values the Jacobi–Davidson method is the method of choice with respect to convergence, accuracy and memory usage.
Authors
CălinIoan Gheorghiu
(Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy)
M. E. Hochstenbach
(Department of Mathematics and Computer Science, TU Eindhoven)
B. Plestenjak
(Department of Mathematics, University of Ljubljana)
J. Rommes
(NXP Semiconductors, The Netherlands)
Keywords
References
See the expanding block below.
Cite this papaer as
C.I. Gheorghiu, M.E. Hochstenbach, B. Plestenjak, J. Rommes, Spectral collocation solutions to multiparameter Mathieu’s system, Appl. Math. Comput., 218 (2012) 1199012000.
doi: 10.1016/j.amc.2012.05.068
http://www.win.tue.nl/~hochsten/pdf/mathieu.pdf
About this paper
Print ISSN
00963003
Online ISSN
Google Scholar Profile
google scholar link
Paper (preprint) in HTML form
Spectral collocation solutions to multiparameter Mathieu’s system
Abstract
Our main aim is the accurate computation of a large number of specified eigenvalues and eigenvectors of Mathieu’s system as a multiparameter eigenvalue problem (MEP). The reduced wave equation, for small deflections, is solved directly without approximations introduced by classical Mathieu functions. We show how for moderate values of the cutoff collocation parameter the QZ algorithm and the Arnoldi method may be applied successfully, while for larger values the JacobiDavidson method is the method of choice with respect to convergence, accuracy and memory usage.
Key words: Mathieu’s system; Chebyshev collocation; multiparameter eigenvalue problems; Jacobi–Davidson; tensor Rayleigh quotient iteration.
1 Introduction
Mathieu’s system is often used in the literature as a motivating example for the introduction of multiparameter eigenvalue problems, see, for example, [25]. This MEP approach of this well known system is the most natural one. The system is obtained if separation of variables is applied to solve the vibration of a fixed elliptic membrane. However, to the best of our knowledge, the accurate numerical solution of Mathieu’s system as a twoparameter eigenvalue problem was never studied in detail. Ruby in his paper [20] provides some examples from science and technology arguing that they deserve accurate solutions of Mathieu’s system. From the paper of Igbokoyi and Tiab [15], as well as from other papers quoted there, it is apparent that the case of an ellipse with the minor axis approaching zero is of tremendous importance in petroleum engineering.
The main aim of this paper is to find a large number (more than four hundred) of even and odd, $\pi $ and $2\pi $, eigenfrequencies and eigenmodes of Mathieu’s system as a MEP very accurately. Neither the radial Mathieu’s equation nor the Mathieu’s system, as a twoparameter eigenvalue problem, have been solved using Chebyshev collocation (pseudospectral) method. This approach, in conjunction with various methods to solve the (discretized) algebraic MEP, is considered in this paper. For these methods, as well as for the normality of the discretization matrices, is paid a particular attention. Thus, for small to moderate values of the cutoff collocation parameter, $N$, the QZ algorithm and shiftandinvert Arnoldi method work satisfactory. For larger $N$ they fail. The remedy is a JacobiDavidson based method which accurately and efficiently solves these cases.
In fact the literature concerning the numerics of the second problem is rather poor. The paper of Neves [16] contains some numerical results along with a Klein oscillation theorem for the multiparameter Mathieu’s system. These numerical results came from an ad hoc method. It involves a shooting scheme based on RungeKutta method used to solve a twopoint boundary value problem. Troesch and Troesch find in their paper [24] two lowest eigenvalues using for the representation of the Mathieu’s functions the Bessel functions. Gutiï¿$\frac{1}{2}$rrezVega et. al. in their paper [9] use the Fourier representation in order to find classical Mathieu’s functions. Without the need of special functions, Wilson and Scharstein [28] use a Fourier collocation method to find a ”wide range” of eigenfrequencies, i.e., first one hundred modes. They do not solve a MEP. Instead they solve a sequence of two generalized eigenvalues and this seems to affects the accuracy of their solutions.
In contrast, the angular Mathieu’s equation is solved by Trefethen in [23] and Weideman and Reddy in [26] by Fourier collocation and thoroughly analyzed by Boyd in his monograph [2]. More recently, Shen and Wang [21] provide approximation results (in Sobolev spaces) for the eigenmodes of the first Mathieu equation. They also solve the second Mathieu equation by a spectral Galerkin method and eventually by a Legendre spectralGalerkin method they solve a Helmholtz and a modified Helmholtz equation.
The paper is organized as follows. In Sect. 2 we introduce the Mathieu system as a MEP, i.e., the four possible differential MEPs. We comment on the twoparameter algebraic eigenvalue problems in Sect. 3 and provide an overview of the JacobiDavidson method to solve such problems in Sect. 4. The Chebyshev collocation discretization as well as a finite difference discretization of the Mathieu’s system as a MEP is considered in Sect. 5. Our numerical results are presented in Sect. 6. Sect. 7 concludes.
2 Mathieu’s system
The coupled system of Mathieu’s angular and radial equations in which $a$ and $q$ are independent parameters will be thought of now as a multiparameter (two parameter) eigenvalue problem. The following four MEPs can be formulated with respect to Mathieu’s system:

•
a $\pi $even problem
$$ (1) 
•
a $2\pi $even problem
$$ (2) 
•
a $\pi $odd problem
$$ (3) 
•
a $2\pi $odd problem
$$ (4) These coupled systems of twopoint boundary value problems come from the problem of a vibrating elliptic membrane $\mathrm{\Omega}$ with fixed boundaries $\partial \mathrm{\Omega}$,
$$\left(\mathrm{\Delta}+{\omega}^{2}\right)\psi (x,y)=0,(x,y)\in \mathrm{\Omega},\psi (x,y)=0,(x,y)\in \partial \mathrm{\Omega},$$ (5) when the separation of the variable, i.e., $\psi (x,y):=F\left(\xi \right)G\left(\eta \right)$ is used in the elliptical coordinates $\xi $ and $\eta ,$
$$ (6) Thus
$${\xi}_{0}:=\text{arccosh}\frac{\alpha}{h},$$ where $h:=\sqrt{{\alpha}^{2}{\beta}^{2}}$ is half the distance between the foci of the membrane. The parameter $q$ is related to the eigenfrequency $\omega $ by
$$q:=\frac{{h}^{2}{\omega}^{2}}{4},$$ and $a$ is the parameter arising in the separation of variables. Volkmer analyses in details in his monograph [25] the above four systems. For these right definite MEP, he provides results concerning the existence and countability of eigenvalues, numbers of zeros of eigenfunctions and the completeness of even and odd sets of eigenmodes. His analytical results are exhaustive. A Klein oscillation theorem is also available in [16] and some other useful comments on the formulations above can be found in [4]. The Mathieu system can also be ”embedded” in the most general setting of the multiparameter eigenvalue problem for ordinary differential equations formulated by Sleeman in [22].
3 Algebraic twoparameter eigenvalue problem
An algebraic twoparameter eigenvalue problem has the form
$$\{\begin{array}{c}{A}_{1}x=\lambda {B}_{1}x+\mu {C}_{1}x,\hfill \\ {A}_{2}y=\lambda {B}_{2}y+\mu {C}_{2}y,\hfill \end{array}$$  (7) 
where ${A}_{i},{B}_{i}$, and ${C}_{i}$ are given ${n}_{i}\times {n}_{i}$ complex matrices, $\lambda ,\mu \in $, $x\in {}^{{n}_{1}}$, and $y\in {}^{{n}_{2}}$. A pair $(\lambda ,\mu )$ is called an eigenvalue if it satisfies (7) for nonzero vectors $x$ and $y$. The tensor product $x\otimes y$ is then the corresponding eigenvector.
The twoparameter eigenvalue problem (7) can be expressed as two coupled generalized eigenvalue problems as follows. On the tensor product space $S:={}_{}{}^{{n}_{1}}\otimes {}^{{n}_{2}}$ of the dimension $N:={n}_{1}{n}_{2}$ we define so called operator determinants
${\mathrm{\Delta}}_{0}$  $=$  ${B}_{1}\otimes {C}_{2}{C}_{1}\otimes {B}_{2},$  
${\mathrm{\Delta}}_{1}$  $=$  ${A}_{1}\otimes {C}_{2}{C}_{1}\otimes {A}_{2},$  
${\mathrm{\Delta}}_{2}$  $=$  ${B}_{1}\otimes {A}_{2}{A}_{1}\otimes {B}_{2}$ 
(for details on the tensor product and relation to the multiparameter eigenvalue problem, see, for example, [1]). The twoparameter eigenvalue problem $(\text{7})$ is nonsingular when ${\mathrm{\Delta}}_{0}$ is nonsingular. In this case the matrices ${\mathrm{\Delta}}_{0}^{1}{\mathrm{\Delta}}_{1}$ and ${\mathrm{\Delta}}_{0}^{1}{\mathrm{\Delta}}_{2}$ commute and (7) is equivalent to a coupled pair of generalized eigenvalue problems
$$\{\begin{array}{c}{\mathrm{\Delta}}_{1}z=\lambda {\mathrm{\Delta}}_{0}z,\hfill \\ {\mathrm{\Delta}}_{2}z=\mu {\mathrm{\Delta}}_{0}z\hfill \end{array}$$  (8) 
for decomposable tensors $z\in S$, $z=x\otimes y$ (see [1]).
There exist some numerical methods for twoparameter eigenvalue problems. If $N$ is small, we can apply the existing numerical methods for the generalized eigenvalue problem to solve the coupled pair (8). An algorithm of this kind, which is based on the QZ algorithm, is presented in [12].
When $N$ is large, it is not feasible to compute all eigenpairs. There are some iterative methods that can be applied to compute some solutions. Most of them require good initial approximations in order to avoid misconvergence. One such method, that we apply in our numerical experiments, is the Tensor Rayleigh Quotient Iteration (TRQI) from [17], which is a generalization of the standard Rayleigh quotient iteration, (see, e.g., [6]). This method computes one eigenpair at a time.
In case when we are interested in more than just one eigenpair and we do not have any initial approximations, a method of choice is a Jacobi–Davidson type method [12]. The stateoftheart, which uses harmonic Ritz values [13], can be used to compute a small number of eigenvalues of (7), which are closest to a given target $({\lambda}_{\mathrm{T}},{\mu}_{\mathrm{T}})\in {}^{2}$. The method is presented in the next section.
4 Overview of Jacobi–Davidson method
For the numerical solution we exploit a Jacobi–Davidson method as developed in [12, 13, 19]. We will now present a brief overview.
In this method the eigenvectors $x$ and $y$ are sought in search spaces and , respectively. There are two main phases: expansion of the subspaces, and extraction of an approximate eigenpair from the search space. First consider the subspace expansion. Suppose that we have approximate eigenvectors $u\approx x$ and $v\approx y$ with corresponding approximate eigenvalue $(\sigma ,\tau )\approx (\lambda ,\mu )$; for instance, the tensor Rayleigh quotient. We are interested in orthogonal improvements $s\u27c2u$ and $t\u27c2v$ such that
${A}_{1}(u+s)$  $=$  $\lambda {B}_{1}(u+s)+\mu {C}_{1}(u+s),$  (9)  
${A}_{2}(v+t)$  $=$  $\lambda {B}_{2}(v+t)+\mu {C}_{2}(v+t).$  (10) 
Let
${r}_{1}$  $=$  $({A}_{1}\sigma {B}_{1}\tau {C}_{1})u,$  
${r}_{2}$  $=$  $({A}_{2}\sigma {B}_{2}\tau {C}_{2})v$ 
be the residuals of vector $u\otimes v$ and value $(\sigma ,\tau )$. We can rewrite (LABEL:problem2aa) and (9) as
$({A}_{1}\sigma {B}_{1}\tau {C}_{1})s$  $=$  ${r}_{1}+(\lambda \sigma ){B}_{1}u+(\mu \tau ){C}_{1}u$  
$+(\lambda \sigma ){B}_{1}s+(\mu \tau ){C}_{1}s,$  
$({A}_{2}\sigma {B}_{2}\tau {C}_{2})t$  $=$  ${r}_{2}+(\lambda \sigma ){B}_{2}v+(\mu \tau ){C}_{2}v$  
$+(\lambda \sigma ){B}_{2}t+(\mu \tau ){C}_{2}t.$ 
We neglect the secondorder correction terms $(\lambda \sigma ){B}_{1}s+(\mu \tau ){C}_{1}s$ and $(\lambda \sigma ){B}_{2}t+(\mu \tau ){C}_{2}t$. Let $V\in {}^{({n}_{1}+{n}_{2})\times 2}$ be a matrix with columns (for reasons of stability, preferably orthonormal) such that
$$\mathrm{s}pan(V)=\mathrm{s}pan(\left[\begin{array}{cc}{B}_{1}u& \\ {B}_{2}v& \end{array}\right],\left[\begin{array}{cc}{C}_{1}u& \\ {C}_{2}v& \end{array}\right]),$$ 
and let $W\in {}^{({n}_{1}+{n}_{2})\times 2}$ be
$$W=\left[\begin{array}{cc}u& 0\\ 0& v\end{array}\right].$$ 
With the oblique projection
$$P=IV{({W}^{T}V)}^{1}{W}^{T}$$ 
onto $\mathrm{s}pan{(V)}^{\perp}$ along $\mathrm{s}pan(W)$, it follows that
$$Pr=r\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}P\left[\begin{array}{cc}{B}_{1}u& \\ {B}_{2}v& \end{array}\right]=P\left[\begin{array}{cc}{C}_{1}u& \\ {C}_{2}v& \end{array}\right]=0.$$  (13) 
Therefore, we can project out the firstorder terms $(\lambda \sigma ){B}_{1}u+(\mu \tau ){C}_{1}u$ and $(\lambda \sigma ){B}_{2}v+(\mu \tau ){C}_{2}v$ using this oblique projection, reformulating (4) and (4) (without the neglected secondorder correction terms) as
$$P\left[\begin{array}{cc}{A}_{1}\sigma {B}_{1}\tau {C}_{1}& 0\\ 0& {A}_{2}\sigma {B}_{2}\tau {C}_{2}\end{array}\right]\left[\begin{array}{cc}s& \\ t& \end{array}\right]=\left[\begin{array}{cc}{r}_{1}& \\ {r}_{2}& \end{array}\right]$$  (14) 
for $s\perp u$ and $t\perp v$. We use (possibly) inexact solutions $s$ and $t$ to this linear system to expand the search spaces and .
Now we focus on the subspace extraction. As introduced in [13], the harmonic Rayleigh–Ritz extraction for the MEP extracts approximate vectors $u$, $v$ and corresponding values $\sigma $ and $\tau $ by imposing the Galerkin conditions
$$\begin{array}{ccc}{A}_{1}u\sigma {B}_{1}u\tau {C}_{1}u\hfill & \u27c2\hfill & ({A}_{1}\sigma {B}_{1}\tau {C}_{1}),\hfill \\ {A}_{2}v\sigma {B}_{2}v\tau {C}_{2}v\hfill & \u27c2\hfill & ({A}_{2}\sigma {B}_{2}\tau {C}_{2}).\hfill \end{array}$$  (15) 
This generally turns out to be a method of choice for interior eigenvalues near a target $(\sigma ,\tau )$. A very basic pseudocode for the method is given in the following algorithm, where RGS stands for repeated Gram–Schmidt, or any other numerically robust method to expand an orthonormal basis.
A Jacobi–Davidson type method for the MEP
Input: Starting vectors ${u}_{1}$, ${v}_{1}$ and a tolerance $\epsilon $
Output: An approximate eigenpair $(\theta ,\eta ,u,v)$ for the MEP
$s=u$, $t=v$, ${U}_{0}=[]$, ${V}_{0}=[]$
for $k=1,2,\mathrm{\dots}$
RGS(${U}_{k1},s$) $\to {U}_{k}$, RGS(${V}_{k1},t$) $\to {V}_{k}$
Extract an approximation $(u,v,\theta ,\eta )$ using (15)
Solve (approximately) $s\u27c2u$, $t\u27c2v$ from (14)
end
5 Discretization of Mathieu’s system as a MEP
Our previous numerical experiments concerning nonstandard, high order and singularly perturbed eigenvalue problems reported in [3, 7, 8] proved that Chebyshev collocation method is fairly accurate, flexible and implementable. It turned out to be superior to the spectral Galerkin or tau method also based on Chebyshev polynomials. The well known monograph of Fornberg [5] provides a thorough analysis with how, when and why this pseudospectral approach works.
Thus, the Chebyshev collocation discretization of MEP (1) reads
$$\{\begin{array}{c}\left({\left(\frac{4}{\pi}\right)}^{2}{\cdot}^{e,\pi}{D}_{n}^{2}+\left(a2q\cdot \mathrm{d}iag\left(\mathrm{cos}\pi \left({X}_{C}+1\right)/2\right)\right)\right)U=0,\\ \left({\left(\frac{2}{{\xi}_{0}}\right)}^{2}{\cdot}^{e,\pi}{D}_{nd}^{2}\left(a2q\cdot \mathrm{d}iag\left(\mathrm{cosh}{\xi}_{0}\left({X}_{C}+1\right)\right)\right)\right)V=0,\end{array}$$  (16) 
where ${}_{}{}^{e,\pi}D_{n}^{2}$ and ${}_{}{}^{e,\pi}D_{nd}^{2}$ are second order differentiation matrices in the Chebyshev nodes of the second kind ${X}_{C}$ , i.e.,
$${X}_{C}:=\{\mathrm{cos}\left(\frac{\left(k1\right)\pi}{N1}\right),k=1,2,\mathrm{\dots},N\}.$$ 
In the symbol ${}_{}{}^{e,\pi}D_{n}^{2}$ the upper indices $e$ and $\pi $ stand for even and $\pi $ period, and the lower index $n$ for the Neumann boundary conditions
$${G}^{\prime}\left(0\right)={G}^{\prime}\left(\pi /2\right)=0,$$ 
which are enforced. Similarly, in ${}_{}{}^{e,\pi}D_{nd}^{2}$ the mixed boundary conditions
$${F}^{\prime}\left(0\right)=F\left({\xi}_{0}\right)=0,$$ 
are introduced, so $n$ comes from the first and $d$ from the second boundary condition. We used the seminal paper of Weideman and Reddy [26] to obtain the entries of these two matrices and the simple and general strategy of Hoepffner [14] to impose all boundary conditions. These matrices are also available in the book of Trefethen [23]. The vectors $U$ and $V$ contain the unknown values of $G$ and $F$ in the nodes ${X}_{C}$.
Thus, the problem (16) is an algebraic MEP of type (7) with $(a,q)$ standing for $(\lambda ,\mu )$. The discretizations for the last three problems (2), (3) and (4) are analogous.
Unfortunately, the matrices ${}_{}{}^{e,\pi}D_{n}^{2}$ and ${}_{}{}^{e,\pi}D_{nd}^{2}$ are dense, nonsymmetric and have high condition numbers (see for instance [8]).
The pseudospectra (see [11] for definition and numerical code) of even problems (1) and (2) are depicted in Fig. 1. This picture shows two mildly nonnormal MEP with decreasing nonnormality for large $a\left(q\right)$.
It is worth noting at this moment that the curves $a\left(q\right)$ represent the solutions of the first SturmLiouville problems in (1) and (2) for $q\in [0,10].$ They are the interlaced quasi ”vertical” curves in Fig. 1. The family of curves $A\left(q\right)$ depicts the solutions of the second SturmLiouville problem in (1) or (2) for the same range of $q$. They are represented by the quasi ”oblique” curves. Their intersections localize the eigenpairs $(a,q)$ of MEP (1). Our Fig. 1 refines Fig. 1 from the paper of Neves [16].
The Chebyshev collocation discretization of MEP (4) reads
$$\{\begin{array}{c}\left({\left(\frac{4}{\pi}\right)}^{2}{\cdot}^{o,2\pi}{D}_{dn}^{2}+\left(a2q\cdot \text{diag}\left(\mathrm{cos}\pi \left({X}_{C}+1\right)/2\right)\right)\right)U=0,\\ \left({\left(\frac{2}{{\xi}_{0}}\right)}^{2}{\cdot}^{o,2\pi}{D}_{d}^{2}\left(a2q\cdot \text{diag}\left(\mathrm{cosh}{\xi}_{0}\left({X}_{C}+1\right)\right)\right)\right)V=0.\end{array}$$  (17) 
In ${}_{}{}^{o,2\pi}D_{dn}^{2}$ the upper index $o$ stands for odd property, the index $2\pi $ for period and $dn$ for the mixed boundary conditions
$$G\left(0\right)={G}^{\prime}\left(\pi /2\right)=0$$ 
. The matrix ${}_{}{}^{o,2\pi}D_{d}^{2}$ with the lower index $d$ involves the symmetric Dirichlet boundary conditions
$$F\left(0\right)=F\left({\xi}_{0}\right)=0.$$ 
The pseudospectra of odd problems (3) and (4) are depicted in Fig. 2.
It shows, two even more normal problem than (1) and (2). The explanation consists in the fact that the Dirichlet boundary conditions
$$F\left(0\right)=F\left({\xi}_{0}\right)=0$$ 
in (3) and (4) induce a sort of symmetry in the differentiation matrices, i.e. the matrix ${}_{}{}^{o,2\pi}D_{d}^{2}$ is more normal than ${}_{}{}^{e,\pi}D_{nd}^{2}{=}^{e,2\pi}{D}_{nd}^{2}$.
In order to evaluate the performances of our strategy we carried out numerical experiments on the finite difference discretization of our differential eigenvalue problems. Thus the usual finite difference of (1) reads
$$\{\begin{array}{c}\left({\left(\frac{2}{\pi}\right)}^{2}{\cdot}^{e,\pi}{D}_{n}^{2,FD}+\left(a2q\cdot \mathrm{d}iag\left(\mathrm{cos}\pi \left({X}_{N+1}\right)\right)\right)\right)U=0,\\ \left({\left(\frac{1}{{\xi}_{0}}\right)}^{2}{\cdot}^{e,\pi}{D}_{nd}^{2,FD}\left(a2q\cdot \mathrm{d}iag\left(\mathrm{cosh}2{\xi}_{0}\left({X}_{N}\right)\right)\right)\right)V=0,\end{array}$$  (18) 
where ${}_{}{}^{e,\pi}D_{n}^{2,FD}$ and ${}_{}{}^{e,\pi}D_{nd}^{2,FD}$ stand for the second order centered finite difference approximation of the second derivative in the $N+1$ equispaced nodes ${X}_{N+1}$.
The matrices ${}_{}{}^{e,\pi}D_{n}^{2,FD}$ and ${}_{}{}^{e,\pi}D_{nd}^{2,FD}$ are now symmetric and tridiagonal of order $N+1$ and $N$ respectively and the Neumann boundary conditions were introduced by mirror imaging technique described in the monograph of Quarteroni et. al. [18], p. 549. Despite these simplifications in (18) the numerical results provided in the next section are obviously inferior to those obtained by Chebyshev collocation (see Tab. 3 and Tab. 4).
It is important to point out at this moment that in their paper [28] Wilson and Scharstein use a Fourier collocation method (and not Galerkin!) in order to discretize the Mathieu’s system. Their shape (trial) functions are some trigonometric functions which implicitly satisfy the boundary conditions but it is not clear from this paper what is the distribution of their nodes and how they are clustered to the boundary. As this paper is detailed in the monograph [29] it seems that they use a uniform grid. Our strategy takes the advantage of the Chebyshev clustering to the boundary.
6 Numerical examples
In our numerical experiments we compute solutions of Mathieu’s systems (1)–(4) using discretizations from Sect. 5. For each of the four systems we know from Volkmer [25] that for every pair of nonnegative indices $(i,j)$ there exists a pair $({a}_{ij},{q}_{ij})$ with nonzero functions ${F}_{ij}$ and ${G}_{ij}$ such that ${G}_{ij}$ has exactly $i$ zeros on $(0,\pi /2)$ and ${F}_{ij}$ has exactly $j$ zeros on $(0,{\xi}_{0})$. This is one way how we can index the solutions.
Another option of indexing comes from the fact that Mathieu’s systems (1)–(4) are related to the problem of a vibrating elliptic membrane with fixed boundaries (5). Each solution $(a,q)$ gives an eigemode of (5) with the eigenfrequency $\omega =2\sqrt{q}/h$. The solutions of (1) and (2) give all even eigenmodes of (5). We order the even eigenmodes so that ${\omega}_{1}^{e}\le {\omega}_{2}^{e}\le \mathrm{\cdots}.$ To each even eigenmode (see, for example, [16] or [28]) we can associate an index $(k,l)$, where $k$ is the number of zeros of $G$ on $(0,\pi )$, and $l$ is the number of zeros of $F$ on $(0,{\xi}_{0})$. The eigenmode is then ${\psi}_{e}^{k,l}(x,y)=F(\xi )G(\eta ).$ In a similar way the solutions of (3) and (4) give the odd eigenmodes ${\psi}_{o}^{k,l}$ of (5).
In particular, if ${F}_{ij}$ and ${G}_{ij}$ are solutions of one of Mathieu’s systems (1)–(4), then
$${F}_{ij}(\xi ){G}_{ij}(\eta )=\{\begin{array}{cc}{\psi}_{e}^{2i,j}(x,y)& \mathrm{for}(\text{1}),\\ {\psi}_{e}^{2i+1,j}(x,y)& \mathrm{for}(\text{2}),\\ {\psi}_{o}^{2i+2,j}(x,y)& \mathrm{for}(\text{3}),\\ {\psi}_{o}^{2i+1,j}(x,y)& \mathrm{for}(\text{4}).\end{array}$$ 
The choice of the method to solve the algebraic MEP’s depends on the requested eigenvalue and the required accuracy. It is clear that if we want to compute a higher eigenfrequency very accurately, we need a larger $N$. Depending on the size of $N$, we propose to use one of the following methods:

a)
EIG$\mathrm{\Gamma}$: When $N$ is small, we can apply the existing numerical methods (for instance eig in Matlab) to the eigenvalue problem
$${\mathrm{\Delta}}_{0}^{1}{\mathrm{\Delta}}_{2}z=\mu z,$$ (19) where the matrix ${\mathrm{\Gamma}}_{2}:={\mathrm{\Delta}}_{0}^{1}{\mathrm{\Delta}}_{2}$ is of size ${N}^{2}\times {N}^{2}$. The obtained eigenvector $z$ is decomposable, i.e., $z=x\otimes y$, and it is easy to compute $x$ and $y$ from $z$.

b)
EIGS$\mathrm{\Gamma}$: When matrix ${\mathrm{\Gamma}}_{2}$ is too large for a), we can apply the implicit shiftandinvert Arnoldi (available as function eigs in Matlab) on (19). One can see that the matrix ${\mathrm{\Gamma}}_{2}$ is quite sparse. It has $N$ full blocks of size $N\times N$ on its diagonal, whereas all nondiagonal $N\times N$ blocks are diagonal matrices. In many cases, when we need just a small number of eigenvalues, b) is more efficient than a) even for small $N$.
If $N$ is very large, this approach is not feasible anymore. The first problem is that the $L$ and $U$ factors of the LU decomposition of the matrix ${\mathrm{\Gamma}}_{2}\sigma I$ are virtually full triangular matrices, and we run out of memory.
Although we could try to use another solver instead of the default LU decomposition in eigs, there is another problem when $N$ is large. Namely, as the matrix ${\mathrm{\Gamma}}_{2}$ has size ${N}^{2}\times {N}^{2}$, the method builds its search space by vectors of size ${N}^{2}$, which is time and memory consuming.

c)
JD$W$: When $N$ is too large for b), we can apply the Jacobi–Davidson method. An advantage of the Jacobi–Davidson method is that it works with matrices and vectors of size $N$. Therefore, the method might be applied when b) is too expensive.
The results were obtained using Matlab 2010b running on Intel Core Duo P8700 2.53GHz processor using 4GB of memory. In this environment, the approach EIGS$\mathrm{\Gamma}$ works up to $N=80$, for larger $N$ we have to use JD$W$. The method EIGS$\mathrm{\Gamma}$ might be more efficient than EIG$\mathrm{\Gamma}$ if many eigenvalues are required.
Example 1
We compare EigElip, which is a Matlab implementation of our method EIGS$\mathrm{\Gamma}$, to the Matlab function runelip by Wilson [27], which was used to compute the eigenfrequencies in [28].
The results in Table 1 show the results, obtained for the computation of the lowest $n$ even eigenfrequencies for the ellipse with given $\alpha $ and $\beta $. Parameters ${N}_{1}$ and ${N}_{2}$ for EigElip specify the number of points used for the discretization of Mathieu’s systems, which might be different for each of the two equations. The number ${N}_{1}$ is used for the angular equation and ${N}_{2}$ is used for the radial equation. The values are chosen so that the computed eigenfrequencies are correct to at least 8 decimal places. The parameter nrts$=({k}_{m},{l}_{m})$ in runelip specifies that the method computes all eigenfrequencies of index $(k,l)$ where $k\le {k}_{m}$ and $l\le {l}_{m}$. The values are minimal possible so that all of the lowest $n$ eigenfrequencies are among the computed ones. The computational times, which are given in seconds, show that the new method is considerably faster than runelip. The values in the last column present the maximum absolute error of the computed eigenvalues. One can see that runelip is not accurate for higher eigenfrequencies.
EigElip  runelip  

$\alpha $  $\beta $  $n$  $({N}_{1},{N}_{2})$  time  error  nrts  time  error 
$2$  $1$  $100$  $(50,25)$  2.2  5e10  (27,6)  14.8  5e10 
$2$  $1$  $200$  $(60,30)$  6.5  5e9  (39,9)  34.9  5e9 
$2$  $1$  $300$  $(70,35)$  17.2  1e8  (48,11)  52.7  1e8 
$4$  $1$  $100$  $(65,25)$  3.4  3e10  (27,6)  18.2  3e10 
$4$  $1$  $200$  $(80,30)$  11.3  5e9  (39,9)  31.1  5e9 
$4$  $1$  $300$  $(90,30)$  19.7  3e3  (39,9)  51.9  3e3 
$8$  $1$  $100$  $(80,20)$  3.2  3e3  (48,3)  14.9  1e9 
$\mathrm{cosh}(2)$  $\mathrm{sinh}(2)$  $100$  $(35,40)$  2.6  3e3  (24,9)  13.4  1e9 
Example 2
In this example we use Jacobi–Davidson with harmonic Ritz values, presented in Section 4, to compute eigenvalues close to a given target. Depending on the region of interest we do this for several targets. The results of this phase is a set of eigenpairs $(({\lambda}_{k},{\mu}_{k}),{x}_{k}\otimes {y}_{k})$ for $k=1,\mathrm{\dots},m$. For each obtained eigenvalue $({\lambda}_{k},{\mu}_{k})$ we compute its index $({i}_{k},{j}_{k})$, where ${i}_{k}$ and ${j}_{k}$ are the number of zeros of ${x}_{k}$ and ${y}_{k}$, respectively. Here we assume that vectors ${x}_{k}$ and ${y}_{k}$ are discrete approximations of continuous curves.
In the second phase we extend the obtained set by the TRQI. We exploit the following property of eigenvectors of Mathieu’s system. Let ${x}_{1}\otimes {y}_{1}$ and ${x}_{2}\otimes {y}_{2}$ be approximate eigenvectors belonging to eigenvalues with indices $({i}_{1},{j}_{1})$ and $({i}_{2},{j}_{2})$, respectively. If ${j}_{1}={j}_{2}$ and ${i}_{1}$ is close to ${i}_{2}$, then ${x}_{1}$ and ${x}_{2}$ do not differ much. The same applies to ${y}_{1}$ and ${y}_{2}$ when ${i}_{1}={i}_{2}$ and ${j}_{1}$ is close to ${j}_{2}$. This is displayed on Fig. 3, where the $x$ part of the eigenvector corresponding to the eigenvalue with index $(k,5)$ for $k=0,\mathrm{\dots},5$ are presented. So, for each pair of eigenvectors, such that ${i}_{1}$ is close to ${i}_{2}$ and ${j}_{1}$ is close to ${j}_{2}$, we can apply TRQI with an initial approximation ${x}_{1}\otimes {y}_{2}$ in order to compute the eigenpair with the index $({i}_{1},{j}_{2})$. This simple approach usually converges in a couple of steps.
We take Chebyshev discretization (16) with matrices of size $100\times 100$ and two targets: $(0,0)$ and $(100,0)$. For each target we do 200 outer iterations of the Jacobi–Davidson method with harmonic Ritz values. As preconditioners we take ${({A}_{i}{\sigma}_{T}{B}_{i}{\tau}_{T}{C}_{i})}^{1}$ for $i=1,2$, where $({\sigma}_{T},{\tau}_{T})$ is the current target. We apply 20 steps of GMRES to solve the correction equations. As a result, we get 36 eigenpairs for the target $(0,0)$ and 28 eigenpairs for the target $(100,0)$. From these eigenpairs we compute additional 20 eigenvalues with the TRQI, so that in the end we have approximations for all eigenpairs of (1) with indices $(i,j)$, such that $i\le 13$ and $j\le 5$.
The computed eigenvalues are presented on Fig. 4 in three different colors. The red and the green eigenvalues were computed using the Jacobi–Davidson method with targets $(0,0)$ and $(100,0)$, respectively, while the blue eigenvalues were computed using the TRQI. One can see that in a large majority the eigenvalues computed by the Jacobi–Davidson method are indeed the closest ones to the given target.
Example 3
Using the method EIGS$\mathrm{\Gamma}$ we can accurately compute higher eigenmodes than the previously reported in the literature.
For example we take $\alpha =4$ and $\beta =1$ and compute the lowest 300 even eigenmodes using EigElip with ${N}_{1}=120$ and ${N}_{2}=40$. The eigenmodes ${\psi}_{e}^{3,8}$ and ${\psi}_{e}^{52,3}$ for the eigenfrequencies ${\omega}_{298}^{e}=24.45490912$ and ${\omega}_{300}^{e}=24.53067377$ are presented in Figures 5 and 6, respectively.
It is important to underline that higher eigenmodes, which require larger matrices, could not be obtained by EIG and EIGS methods due to memory limitations, while JS$W$ was able to compute these eigenmodes up to the required accuracy.
Using ${N}_{1}={N}_{2}=500$ and the target $(1,7500)$ we computed the even eigenfrequency of the ellipse with $\alpha =2$ and $\beta =1$ that is closest to the target ${\omega}_{T}=100$. The result is $\omega =99.97702290$. Its eigenmode ${\psi}_{e}^{41,25}$ is presented on Fig. 7.
Example 4
There are certain applications that require the eigenmodes of an ellipse with a very large ratio $\alpha /\beta $. See, for instance, [15].
In this example we show that such problems can also be solved efficiently via the MEP approach. If we take the ellipse with $\alpha =1000$ and $\beta =1$, and set ${N}_{1}=200$ and ${N}_{2}=15$, then EigElip returns the 10 lowest even eigenfrequencies in Table 2. The 6th lowest even eigenmode is shown in Fig. 8.
$n$  eigenfrequency  $n$  eigenfrequency 

1  1.57129649  6  1.57630126 
2  1.57229680  7  1.57730317 
3  1.57329744  8  1.57830539 
4  1.57429840  9  1.57930793 
5  1.57529967  10  1.58031079 
Let us remark that runelip fails to compute the 10 lowest eigenfrequencies. It returns $0$, followed by 4 eigenfrequencies smaller than 1. Such results provide confidence in the validity of our approach.
Example 5
We compare the accuracy of the computed eigenfrequencies if we discretize the Mathieu’s system (1) by the Chebyshev collocation discretization (16) or by the standard finite differences (18). We take $\alpha =2$, $\beta =1$ and compute even eigenfrequencies ${\omega}_{1}^{e}$, ${\omega}_{50}^{e}$, and ${\omega}_{100}^{e}$. The errors for the Chebyshev collocation and for the finite differences are collected in Tables 3 and 4, respectively. It is obvious that in spite of better conditioned matrices involved (symmetric, tridiagonal, etc.) finite differences require much larger matrices to obtain accurate results. On the other hand, using the Chebyshev collocation we can compute eigenvalues quite accurately with relatively small matrices.
$N$  error for ${\omega}_{1}^{e}$  error for ${\omega}_{50}^{e}$  error for ${\omega}_{100}^{e}$ 

(20,10)  1.8e10  1.4e2  2.8e3 
(30,15)  8.1e14  5.9e6  1.4e4 
(40,20)  9.5e14  3.0e10  8.0e8 
(50,25)  2.1e13  1.2e14  1.2e11 
$N$  error for ${\omega}_{1}^{e}$  error for ${\omega}_{50}^{e}$  error for ${\omega}_{100}^{e}$ 

200  5.9e3  5.7e2  3.9e2 
400  3.0e3  2.6e2  1.2e2 
800  1.5e3  1.2e2  4.5e3 
1600  7.4e4  6.0e3  1.8e3 
7 Conclusions
The Mathieu’s system as a MEP is accurately solved for domains with geometrical aspect ranging from a circle to very deformed ellipse, i.e. the ratio of the major to minor axes of ellipses equals ${10}^{3}$. This means stability with respect to the geometry of the problem.
Accurate numerical computation of high frequencies is much harder than for low frequencies. We introduced a hierarchy of numerical methods that can deal with the corresponding algebraic eigenvalue problems for increasing $N$, and we are able to compute eigenfrequencies and the corresponding eigenmodes from the first ones to the order of about ${10}^{4}$. However, the accuracy varies from a quasi spectral one for the lowest mode to a moderate one for the highest mode.
With respect to the accuracy as well as to the time required, our algorithm is superior to those reported in literature. It is also stable with respect to the degree of the spectral approximation, i.e. $N$, as is apparent from our reported numerical experiments.
All in all, our novel algorithm can be used to solve the MEP associated to Mathieu’s system corresponding to a large variety of geometrical settings.
References
 [1] F.V. Atkinson, Multiparameter eigenvalue problems, Academic Press, New York, (1972).
 [2] J.P. Boyd, Chebyshev and Fourier Spectral Methods, SpringerVerlag, (1989).
 [3] F.I. Dragomirescu, C.I. Gheorghiu, Analytical and numerical solutions to an electrohydrodynamic stability problem, Appl. Math. Comput., 216, 37183727, (2010)
 [4] S.R. Finch. Mathieu eigenvalues. algo.inria.fr/csolve/mthu.pdf, (2008)
 [5] B. Fornberg. A Practical Guide to Pseudospectral Methods, Cambridge University Press, (1998).
 [6] G.H. Golub, C.F. Van Loan. Matrix Computations, 3rd ed., The Johns Hopkins University Press, Baltimore, MD, (1996).
 [7] C.I. Gheorghiu, Spectral Methods for Differential Problems, Casa Cartii de Stiinta, (2007).
 [8] C.I. Gheorghiu, F.I. Dragomirescu. Spectral methods in linear stability. Application to thermal convection with variable gravity field. Appl. Numer. Math. 59, 12901302, (2009)
 [9] J. Gutiï¿$\frac{1}{2}$rrezVega, S. Chï¿$\frac{1}{2}$vezCerda, R. RodriguezDagnino. Free oscillations in an elliptic mambrane. Rev. Mex. Fiz., 45,613622, (1999).
 [10] M.E. Hochstenbach, B. Plestenjak. A Jacobi–Davidson type method for a right definite twoparameter eigenvalue problem. SIAM J. Matrix Anal. Appl., 24, 392–410, (2002).
 [11] M.E. Hochstenbach, B. Plestenjak. Backward error, condition numbers, and pseudospectra for the multiparameter eigenvalue problem. Linear Algebra Appl., 375, 63–81, (2003).
 [12] M.E. Hochstenbach, T. Košir, and B. Plestenjak. A Jacobi–Davidson type method for the nonsingular twoparameter eigenvalue problem, SIAM J. Matrix Anal. Appl., 26, 477–497, (2005).
 [13] M.E. Hochstenbach, B. Plestenjak. Harmonic RayleighRitz for the multiparameter eigenvalue problem, Electron. Trans. Numer. Anal., 29 81–96, (2008).
 [14] J. Hoepffner. Implementation of boundary conditions. www.fukagata.mech.keio.ac.jp/ jerome/, (2007)
 [15] A. O. Igbokoyi, D. Tiab. New Method of Well Test Analysis in Naturally Fractured Reservoirs Based on Elliptical Flow. JCPT, 49, 1–15,(2010).
 [16] A.G.M. Neves. Eigenmodes and eigenfrequencies of vibrating elliptic membranes: A Klein oscillation theorem and numerical calculations. Commun. Pure Appl. Anal., 9, 611–624, (2004).
 [17] B. Plestenjak. A continuation method for a right definite twoparameter eigenvalue problem, SIAM J. Matrix Anal. Appl., 21, 1163–1184, (2000).
 [18] A. Quarteroni, R. Sacco, F. Saleri. Numerical Mathematics, 2nd ed., Springer Berlin Heidelberg, (2007).
 [19] J. Rommes. Arnoldi and JacobiDavidson methods for generalized eigenvalue problems $Ax=\lambda Bx$ with $B$ singular, Math. Comput., 77, 995–1015, (2008).
 [20] L. Ruby, Applications of the Mathieu equation. Am. J. Phys., 64, 3944, (1996).
 [21] J. Shen, LL. Wang, On spectral approximation in elliptical geometries using Mathieu functions. Math. Comput., 78, 815884, (2009).
 [22] B. D. Sleeman, Multiparameter spectral theory and separation of variables. J. Phys. A: Math. Theor. 41, (2008)015209
 [23] L.N. Trefethen, Spectral Methods in MATLAB, SIAM, (2000).
 [24] B. A. Troesch, H. R. Troesch, Eigenfrequencies of an Elliptic Membrane, Math. Comput., 27, 755765, (1973)
 [25] H. Volkmer, Multiparameter Problems and Expansion Theorems, Lecture Notes in Math. 1356, SpringerVerlag, NewYork, (1988)
 [26] J.A.C. Weideman, S.C. Reddy. A MATLAB differentiation matrix suite. ACM Trans. Math. Softw., 26, 465–519, (2000)
 [27] H. B. Wilson. Vibration modes of an elliptic membrane. Available from http://www.mathworks.com/matlabcentral/fileexchange. MATLAB File Exchange, The MathWorks, Natick, (2004).
 [28] H. B. Wilson, R. W. Scharstein, Computing elliptic membrane high frequencies by Mathieu and Galerkin methods, J. Eng. Math, 57, 4155, (2007)
 [29] H. B. Wilson, L. S. Turcotte, D. Halpern. Advanced Mathematics and Mechanics Applications Using MATLAB, 3rd ed., Boca Raton: Chapman and Hall/CRC, (2003).
[1] H. Volkmer, Multiparameter Problems and Expansion Theorems, Lecture Notes in Math. 1356, SpringerVerlag, New York, 1988.
[2] J. Meixner, F. W. Schafke, Mathieusche Funktionen und Spharoidfunktionen, SpringerVerlag, 1954.
[3] L. Ruby, Applications of the Mathieu equation, Am. J. Phys. 64 (1996) 39–44.
[4] A. O. Igbokoyi, D. Tiab, New method of well test analysis in naturally fractured reservoirs based on elliptical flow, J. Can. Pet. Technol. 49 (2010) 1–15.
[5] A. G. M. Neves, Eigenmodes and eigenfrequencies of vibrating elliptic membranes: A Klein oscillation theorem and numerical calculations, Commun. Pure Appl. Anal. 9 (2004) 611–624.
[6] B. A. Troesch, H. R. Troesch, Eigenfrequencies of an elliptic membrane, Math. Comput. 27 (1973) 755–765.
[7] J. GutierrezVega, S. ChavezCerda, R. RodrıguezDagnino, Free oscillations in an elliptic mambrane, Rev. Mex. Fiz. 45 (1999) 613–622.
[8] H. B. Wilson, R. W. Scharstein, Computing elliptic membrane high frequencies by Mathieu and Galerkin methods, J. Eng. Math. 57 (2007) 41– 55.
[9] L. N. Trefethen, Spectral Methods in MATLAB, SIAM, Philadelphia, 2000.
[10] J. A. C. Weideman, S. C. Reddy, A MATLAB differentiation matrix suite, ACM Trans. Math. Softw. 26 (2000) 465–519.
[11] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Lecture Notes in Engineering 49, SpringerVerlag, Berlin, 1989.
[12] J. Shen, LL. Wang, On spectral approximation in elliptical geometries using Mathieu functions, Math. Comput. 78 (2009) 815–884.
[13] S. R. Finch, Mathieu eigenvalues, algo.inria.fr/csolve/mthu.pdf (2008).
[14] B. D. Sleeman, Multiparameter spectral theory and separation of variables, J. Phys. A: Math. Theor. 41 (2008) 1–20.
[15] F. V. Atkinson, Multiparameter Eigenvalue Problems, Academic Press, New York, 1972.
[16] M. E. Hochstenbach, T. Kosir, B. Plestenjak, A Jacobi–Davidson type method for the nonsingular twoparameter eigenvalue problem, SIAM J. Matrix Anal. Appl. 26 (2005) 477–497.
[17] B. Plestenjak, A continuation method for a right definite twoparameter eigenvalue problem, SIAM J. Matrix Anal. Appl. 21 (2000) 1163–1184.
[18] G. H. Golub, C. F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, Baltimore, 1996.
[19] M. E. Hochstenbach, B. Plestenjak, Harmonic RayleighRitz for the multiparameter eigenvalue problem, Electron. Trans. Numer. Anal. 29 (2008) 81–96.
[20] J. Rommes, Arnoldi and JacobiDavidson methods for generalized eigenvalue problems Ax = λBx with B singular, Math. Comput. 77 (2008) 995–1015.
[21] F. I. Dragomirescu, C. I. Gheorghiu, Analytical and numerical solutions to an electrohydrodynamic stability problem, Appl. Math. Comput. 216 (2010) 3718–3727.
[22] C. I. Gheorghiu, Spectral Methods for Differential Problems, Casa Cartii de Stiinta, ClujNapoca, 2007.
[23] C. I. Gheorghiu, F. I. Dragomirescu, Spectral methods in linear stability. Application to thermal convection with variable gravity field, Appl. Numer. Math. 59 (2009) 1290–1302.
[24] B. Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge University Press, Cambridge, 1998.
[25] J. Hoepffner, Implementation of boundary conditions, www.fukagata. mech.keio.ac.jp/~jerome/ (2007).
[26] M. E. Hochstenbach, B. Plestenjak, Backward error, condition numbers, and pseudospectra for the multiparameter eigenvalue problem, Lin. Alg. Appl. 375 (2003) 63–81.
[27] A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, second ed., Texts in Applied Mathematics 47, Springer, Berlin Heidelberg, 2007.
[28] H. B. Wilson, L. S. Turcotte, D. Halpern, Advanced Mathematics and Mechanics Applications Using MATLAB, third ed., Chapman and Hall/CRC, Boca Raton, 2003.
[29] H. B. Wilson, Vibration modes of an elliptic membrane. Available from http://www.mathworks.com/matlabcentral/fileexchange. MATLAB File Exchange, The MathWorks, Natick, 2004.