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 two-parameter 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 two-parameter 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 cut-off collocation parameter, $N$, the QZ algorithm and shift-and-invert Arnoldi method work satisfactory. For larger $N$ they fail. The remedy is a Jacobi-Davidson 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 Runge-Kutta method used to solve a two-point 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}$rrez-Vega 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 spectral-Galerkin 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 two-parameter algebraic eigenvalue problems in Sect. 3 and provide an overview of the Jacobi-Davidson 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.

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 two-point 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].

An algebraic two-parameter eigenvalue problem has the form

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 two-parameter 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

(for details on the tensor product and relation to the multiparameter eigenvalue problem, see, for example, [1]). The two-parameter eigenvalue problem $(\text{<a href="\#S3.E7" title="In 3 Algebraic two-parameter eigenvalue problem ‣ Spectral collocation solutions to multiparameter Mathieu’s system" class="ltx\_ref ltx\_font\_italic"><span class="ltx\_text ltx\_ref\_tag">7</span></a>})$ 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

for decomposable tensors $z\in S$, $z=x\otimes y$ (see [1]).

There exist some numerical methods for two-parameter 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 state-of-the-art, 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.

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⟂u$ and $t⟂v$ such that

Let

be the residuals of vector $u\otimes v$ and value $(\sigma ,\tau )$. We can rewrite (LABEL:problem2aa) and (9) as

We neglect the second-order 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

and let $W\in {}^{(n_1+n_2)\times 2}$ be

With the oblique projection

onto $\mathrm{s}pan{(V)}^{\perp }$ along $\mathrm{s}pan(W)$, it follows that

Therefore, we can project out the first-order 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 second-order correction terms) as

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

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.

Our previous numerical experiments concerning non-standard, 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

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

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

which are enforced. Similarly, in ${}_{}{}^{e,\pi }D_{nd}^2$ the mixed boundary conditions

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, non-symmetric 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 non-normal MEP with decreasing non-normality 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 Sturm-Liouville 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 Sturm-Liouville 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

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

. The matrix ${}_{}{}^{o,2\pi }D_d^2$ with the lower index $d$ involves the symmetric Dirichlet boundary conditions

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

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

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.

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

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:

1. 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$.

2. b)

   EIGS-$\mathrm{\Gamma }$: When matrix ${\mathrm{\Gamma }}_2$ is too large for a), we can apply the implicit shift-and-invert 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 non-diagonal $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.

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

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.