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 the classical Mathieu functions. We show how for moderate values of the cut-off 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ălin-Ioan 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

Mathieu’s system; Chebyshev collocation; multiparameter eigenvalue problem; Jacobi–Davidson method; tensor Rayleigh quotient iteration.

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) 11990-12000.
doi: 10.1016/j.amc.2012.05.068

PDF

http://www.win.tue.nl/~hochsten/pdf/mathieu.pdf

About this paper

Print ISSN

0096-3003

Online ISSN
Google Scholar Profile

google scholar link

Paper (preprint) in HTML form

Spectral collocation solutions to multiparameter Mathieu’s system

Spectral collocation solutions to multiparameter Mathieu’s system

C.I. Gheorghiu
“T. Popoviciu” Institute of Numerical Analysis,
PO Box 68, 3400 Cluj-Napoca 1, Romania
E-mail: ghcalin@ictp.acad.ro
M. E. Hochstenbach
Department of Mathematics and Computing Science, TU Eindhoven
PO Box 513, 5600 MB Eindhoven, The Netherlands
URL: www.win.tue.nl/~hochsten/
B. Plestenjak
Department of Mathematics, University of Ljubljana
Jadranska 19, SI-1000 Ljubljana, Slovenia
E-mail: bor.plestenjak@fmf.uni-lj.si
J. Rommes
NXP Semiconductors, The Netherlands
E-mail: rommes@gmail.com
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 cut-off collocation parameter the QZ 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.

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 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, π and 2π, 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ï¿12rrez-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.

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 π-even problem

    {G′′(η)+(a2qcos2η)G(η)=0,0<η<π/2,G(0)=G(π/2)=0,F′′(ξ)(a2qcosh2ξ)F(ξ)=0,0<ξ<ξ0,F(0)=F(ξ0)=0, (1)
  • a 2π-even problem

    {G′′(η)+(a2qcos2η)G(η)=0,0<η<π/2,G(0)=G(π/2)=0,F′′(ξ)(a2qcosh2ξ)F(ξ)=0, 0<ξ<ξ0,F(0)=F(ξ0)=0, (2)
  • a π-odd problem

    {G′′(η)+(a2qcos2η)G(η)=0, 0<η<π/2,G(0)=G(π/2)=0,F′′(ξ)(a2qcosh2ξ)F(ξ)=0, 0<ξ<ξ0,F(0)=F(ξ0)=0, (3)
  • a 2π-odd problem

    {G′′(η)+(a2qcos2η)G(η)=0, 0<η<π/2,G(0)=G(π/2)=0,F′′(ξ)(a2qcosh2ξ)F(ξ)=0, 0<ξ<ξ0,F(0)=F(ξ0)=0. (4)

    These coupled systems of two-point boundary value problems come from the problem of a vibrating elliptic membrane Ω with fixed boundaries Ω,

    (Δ+ω2)ψ(x,y)=0,(x,y)Ω,ψ(x,y)=0,(x,y)Ω, (5)

    when the separation of the variable, i.e., ψ(x,y):=F(ξ)G(η) is used in the elliptical coordinates ξ and η,

    x:=hcoshξcosη,y:=hsinhξsinη,0ξ<+,0η<2π. (6)

    Thus

    ξ0:=arccoshαh,

    where h:=α2β2 is half the distance between the foci of the membrane. The parameter q is related to the eigenfrequency ω by

    q:=h2ω24,

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

An algebraic two-parameter eigenvalue problem has the form

{A1x=λB1x+μC1x,A2y=λB2y+μC2y, (7)

where Ai,Bi, and Ci are given ni×ni complex matrices, λ,μ, xn1, and yn2. A pair (λ,μ) is called an eigenvalue if it satisfies (7) for nonzero vectors x and y. The tensor product xy 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:=n1n2 of the dimension N:=n1n2 we define so called operator determinants

Δ0 = B1C2C1B2,
Δ1 = A1C2C1A2,
Δ2 = B1A2A1B2

(for details on the tensor product and relation to the multiparameter eigenvalue problem, see, for example, [1]). The two-parameter eigenvalue problem (7) is nonsingular when Δ0 is nonsingular. In this case the matrices Δ01Δ1 and Δ01Δ2 commute and (7) is equivalent to a coupled pair of generalized eigenvalue problems

{Δ1z=λΔ0z,Δ2z=μΔ0z (8)

for decomposable tensors zS, z=xy (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 (λT,μT)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 ux and vy with corresponding approximate eigenvalue (σ,τ)(λ,μ); for instance, the tensor Rayleigh quotient. We are interested in orthogonal improvements su and tv such that

A1(u+s) = λB1(u+s)+μC1(u+s), (9)
A2(v+t) = λB2(v+t)+μC2(v+t). (10)

Let

r1 = (A1σB1τC1)u,
r2 = (A2σB2τC2)v

be the residuals of vector uv and value (σ,τ). We can rewrite (LABEL:problem2aa) and (9) as

(A1σB1τC1)s = r1+(λσ)B1u+(μτ)C1u
+(λσ)B1s+(μτ)C1s,
(A2σB2τC2)t = r2+(λσ)B2v+(μτ)C2v
+(λσ)B2t+(μτ)C2t.

We neglect the second-order correction terms (λσ)B1s+(μτ)C1s and (λσ)B2t+(μτ)C2t. Let V(n1+n2)×2 be a matrix with columns (for reasons of stability, preferably orthonormal) such that

span(V)=span([B1uB2v],[C1uC2v]),

and let W(n1+n2)×2 be

W=[u00v].

With the oblique projection

P=IV(WTV)1WT

onto span(V) along span(W), it follows that

Pr=randP[B1uB2v]=P[C1uC2v]=0. (13)

Therefore, we can project out the first-order terms (λσ)B1u+(μτ)C1u and (λσ)B2v+(μτ)C2v using this oblique projection, reformulating (4) and (4) (without the neglected second-order correction terms) as

P[A1σB1τC100A2σB2τC2][st]=[r1r2] (14)

for su and tv. 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 σ and τ by imposing the Galerkin conditions

A1uσB1uτC1u(A1σB1τC1),A2vσB2vτC2v(A2σB2τC2). (15)

This generally turns out to be a method of choice for interior eigenvalues near a target (σ,τ). 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 u1, v1 and a tolerance ε
Output: An approximate eigenpair (θ,η,u,v) for the MEP
s=u,  t=v,  U0=[],  V0=[]
for k=1,2,
M RGS(Uk1,s) Uk,  RGS(Vk1,t) Vk
M Extract an approximation (u,v,θ,η) using (15)
M Solve (approximately) su, tv from (14)
end

5 Discretization of Mathieu’s system as a MEP

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

{((4π)2e,πDn2+(a2qdiag(cosπ(XC+1)/2)))U=0,((2ξ0)2e,πDnd2(a2qdiag(coshξ0(XC+1))))V=0, (16)

where Dn2e,π and Dnd2e,π are second order differentiation matrices in the Chebyshev nodes of the second kind XC , i.e.,

XC:={cos((k1)πN1),k=1,2,,N}.

In the symbol Dn2e,π the upper indices e and π stand for even and π period, and the lower index n for the Neumann boundary conditions

G(0)=G(π/2)=0,

which are enforced. Similarly, in Dnd2e,π the mixed boundary conditions

F(0)=F(ξ0)=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 XC.

Thus, the problem (16) is an algebraic MEP of type (7) with (a,q) standing for (λ,μ). The discretizations for the last three problems (2), (3) and (4) are analogous.

Unfortunately, the matrices Dn2e,π and Dnd2e,π 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(q).

Refer to caption
Figure 1: The overlapped pseudospectra of problems (1) and (2), N=24,α=cosh(2),β=sinh(2)

It is worth noting at this moment that the curves a(q) represent the solutions of the first Sturm-Liouville problems in (1) and (2) for q[0,10]. They are the interlaced quasi ”vertical” curves in Fig. 1. The family of curves A(q) 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

{((4π)2o,2πDdn2+(a2qdiag(cosπ(XC+1)/2)))U=0,((2ξ0)2o,2πDd2(a2qdiag(coshξ0(XC+1))))V=0. (17)

In Ddn2o,2π the upper index o stands for odd property, the index 2π for period and dn for the mixed boundary conditions

G(0)=G(π/2)=0

. The matrix Dd2o,2π with the lower index d involves the symmetric Dirichlet boundary conditions

F(0)=F(ξ0)=0.

The pseudospectra of odd problems (3) and (4) are depicted in Fig. 2.

Figure 2: The overlapped pseudospectra of problems (3) and (4), N=24,α=cosh(2),β=sinh(2)
Refer to caption

It shows, two even more normal problem than (1) and (2). The explanation consists in the fact that the Dirichlet boundary conditions

F(0)=F(ξ0)=0

in (3) and (4) induce a sort of symmetry in the differentiation matrices, i.e. the matrix Dd2o,2π is more normal than Dnd2e,π=e,2πDnd2.

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

{((2π)2e,πDn2,FD+(a2qdiag(cosπ(XN+1))))U=0,((1ξ0)2e,πDnd2,FD(a2qdiag(cosh2ξ0(XN))))V=0, (18)

where Dn2,FDe,π and Dnd2,FDe,π stand for the second order centered finite difference approximation of the second derivative in the N+1 equispaced nodes XN+1.

The matrices Dn2,FDe,π and Dnd2,FDe,π 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 (aij,qij) with nonzero functions Fij and Gij such that Gij has exactly i zeros on (0,π/2) and Fij has exactly j zeros on (0,ξ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 ω=2q/h. The solutions of (1) and (2) give all even eigenmodes of (5). We order the even eigenmodes so that ω1eω2e. 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,π), and l is the number of zeros of F on (0,ξ0). The eigenmode is then ψek,l(x,y)=F(ξ)G(η). In a similar way the solutions of (3) and (4) give the odd eigenmodes ψok,l of (5).

In particular, if Fij and Gij are solutions of one of Mathieu’s systems (1)–(4), then

Fij(ξ)Gij(η)={ψe2i,j(x,y)for(1),ψe2i+1,j(x,y)for(2),ψo2i+2,j(x,y)for(3),ψo2i+1,j(x,y)for(4).

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-Γ: When N is small, we can apply the existing numerical methods (for instance eig in Matlab) to the eigenvalue problem

    Δ01Δ2z=μz, (19)

    where the matrix Γ2:=Δ01Δ2 is of size N2×N2. The obtained eigenvector z is decomposable, i.e., z=xy, and it is easy to compute x and y from z.

  2. b)

    EIGS-Γ: When matrix Γ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 Γ2 is quite sparse. It has N full blocks of size N×N on its diagonal, whereas all non-diagonal N×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 Γ2σ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 Γ2 has size N2×N2, the method builds its search space by vectors of size N2, 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-Γ works up to N=80, for larger N we have to use JD-W. The method EIGS-Γ might be more efficient than EIG-Γ if many eigenvalues are required.

Example 1

We compare EigElip, which is a Matlab implementation of our method EIGS-Γ, 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 α and β. Parameters N1 and N2 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 N1 is used for the angular equation and N2 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=(km,lm) in runelip specifies that the method computes all eigenfrequencies of index (k,l) where kkm and llm. 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
α β n (N1,N2) time error nrts time error
2 1 100 (50,25) 2.2 5e-10 (27,6) 14.8 5e-10
2 1 200 (60,30) 6.5 5e-9 (39,9) 34.9 5e-9
2 1 300 (70,35) 17.2 1e-8 (48,11) 52.7 1e-8
4 1 100 (65,25) 3.4 3e-10 (27,6) 18.2 3e-10
4 1 200 (80,30) 11.3 5e-9 (39,9) 31.1 5e-9
4 1 300 (90,30) 19.7 3e-3 (39,9) 51.9 3e-3
8 1 100 (80,20) 3.2 3e-3 (48,3) 14.9 1e-9
cosh(2) sinh(2) 100 (35,40) 2.6 3e-3 (24,9) 13.4 1e-9
Table 1: Comparison of EigElip and runelip.
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 ((λk,μk),xkyk) for k=1,,m. For each obtained eigenvalue (λk,μk) we compute its index (ik,jk), where ik and jk are the number of zeros of xk and yk, respectively. Here we assume that vectors xk and yk 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 x1y1 and x2y2 be approximate eigenvectors belonging to eigenvalues with indices (i1,j1) and (i2,j2), respectively. If j1=j2 and i1 is close to i2, then x1 and x2 do not differ much. The same applies to y1 and y2 when i1=i2 and j1 is close to j2. This is displayed on Fig. 3, where the x part of the eigenvector corresponding to the eigenvalue with index (k,5) for k=0,,5 are presented. So, for each pair of eigenvectors, such that i1 is close to i2 and j1 is close to j2, we can apply TRQI with an initial approximation x1y2 in order to compute the eigenpair with the index (i1,j2). This simple approach usually converges in a couple of steps.

Refer to caption
Figure 3: x-part of eigenvector of (1) corresponding to the eigenvalue with index (k,5) for k=0,,5.

We take Chebyshev discretization (16) with matrices of size 100×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 (AiσTBiτTCi)1 for i=1,2, where (σT,τ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 i13 and j5.

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.

Refer to caption
Figure 4: All eigenvalues of (1) with indices (i,j) for i13 and j5, computed by the Jacobi–Davidson method using targets (0,0) and (100,0), and extended by the TRQI.
Example 3

Using the method EIGS-Γ we can accurately compute higher eigenmodes than the previously reported in the literature.

For example we take α=4 and β=1 and compute the lowest 300 even eigenmodes using EigElip with N1=120 and N2=40. The eigenmodes ψe3,8 and ψe52,3 for the eigenfrequencies ω298e=24.45490912 and ω300e=24.53067377 are presented in Figures 5 and 6, respectively.

Refer to caption
Figure 5: Eigenmode ψe3,8 for the ellipse with α=4 and β=1.
Refer to caption
Figure 6: Eigenmode ψe52,3 for the ellipse with α=4 and β=1.

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 N1=N2=500 and the target (1,7500) we computed the even eigenfrequency of the ellipse with α=2 and β=1 that is closest to the target ωT=100. The result is ω=99.97702290. Its eigenmode ψe41,25 is presented on Fig. 7.

Refer to caption
Figure 7: Eigenmode ψe41,25 for the ellipse with α=2 and β=1 that has its eigenfrequency closest to ω=100.
Example 4

There are certain applications that require the eigenmodes of an ellipse with a very large ratio α/β. 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 α=1000 and β=1, and set N1=200 and N2=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
Table 2: The lowest 10 even eigenfrequencies for the ellipse with α=1000 and β=1.
Refer to caption
Figure 8: Eigenmode ψe6,1 for the ellipse with α=1000 and β=1.

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 α=2, β=1 and compute even eigenfrequencies ω1e, ω50e, and ω100e. 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 ω1e error for ω50e error for ω100e
(20,10) 1.8e-10 1.4e-2 2.8e-3
(30,15) 8.1e-14 5.9e-6 1.4e-4
(40,20) 9.5e-14 3.0e-10 8.0e-8
(50,25) 2.1e-13 1.2e-14 1.2e-11
Table 3: Accuracy of the Chebyshev collocation.
N error for ω1e error for ω50e error for ω100e
200 5.9e-3 5.7e-2 3.9e-2
400 3.0e-3 2.6e-2 1.2e-2
800 1.5e-3 1.2e-2 4.5e-3
1600 7.4e-4 6.0e-3 1.8e-3
Table 4: Accuracy of finite differences.

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 103. 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 104. 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, Springer-Verlag, (1989).
  • [3] F.I. Dragomirescu, C.I. Gheorghiu, Analytical and numerical solutions to an electrohydrodynamic stability problem, Appl. Math. Comput., 216, 3718-3727, (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, 1290-1302, (2009)
  • [9] J. Gutiï¿12rrez-Vega, S. Chï¿12vez-Cerda, R. Rodriguez-Dagnino. Free oscillations in an elliptic mambrane. Rev. Mex. Fiz., 45,613-622, (1999).
  • [10] M.E. Hochstenbach, B. Plestenjak. A Jacobi–Davidson type method for a right definite two-parameter 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 two-parameter eigenvalue problem, SIAM J. Matrix Anal. Appl., 26, 477–497, (2005).
  • [13] M.E. Hochstenbach, B. Plestenjak. Harmonic Rayleigh-Ritz 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 two-parameter 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 Jacobi-Davidson methods for generalized eigenvalue problems Ax=λBx with B singular, Math. Comput., 77, 995–1015, (2008).
  • [20] L. Ruby, Applications of the Mathieu equation. Am. J. Phys., 64, 39-44, (1996).
  • [21] J. Shen, L-L. Wang, On spectral approximation in elliptical geometries using Mathieu functions. Math. Comput., 78, 815-884, (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, 755-765, (1973)
  • [25] H. Volkmer, Multiparameter Problems and Expansion Theorems, Lecture Notes in Math. 1356, Springer-Verlag, New-York, (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, 41-55, (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, Springer-Verlag, New York, 1988.
[2] J. Meixner, F. W. Schafke, Mathieusche Funktionen und Spharoidfunktionen, Springer-Verlag, 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. Gutierrez-Vega, S. Chavez-Cerda, R. Rodrıguez-Dagnino, 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, Springer-Verlag, Berlin, 1989.
[12] J. Shen, L-L. 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 two-parameter eigenvalue problem, SIAM J. Matrix Anal. Appl. 26 (2005) 477–497.
[17] B. Plestenjak, A continuation method for a right definite two-parameter 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 Rayleigh-Ritz for the multiparameter eigenvalue problem, Electron. Trans. Numer. Anal. 29 (2008) 81–96.
[20] J. Rommes, Arnoldi and Jacobi-Davidson 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, Cluj-Napoca, 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.

2012

Related Posts