SIAM J. NUMER. ANAL. Vol. 47, No. 5, pp. 3910–3937

c 2009 Society for Industrial and Applied Mathematics 

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF THE SPECTRAL ELEMENT METHOD∗ MARK AINSWORTH† AND HAFIZ ABDUL WAJID† Abstract. If the nodes for the spectral element method are chosen to be the Gauss–Legendre– Lobatto points and a Lagrange basis is used, then the resulting mass matrix is diagonal and the method is sometimes then described as the Gauss-point mass lumped finite element scheme. We study the dispersive behavior of the scheme in detail and provide both a qualitative description of the nature of the dispersive and dissipative behavior of the scheme along with precise quantitative statements of the accuracy in terms of the mesh-size and the order of the scheme. We prove that (a) the Gauss-point mass lumped scheme (i.e., spectral element method) tends to exhibit phase lag whereas the (consistent) finite element scheme tends to exhibit phase lead; (b) the absolute accuracy of the spectral element scheme is 1/p times better than that of the finite element scheme despite the use of numerical integration; (c) when the order p, the mesh-size h, and the frequency of the wave ω satisfy 2p + 1 ≈ ωh the true wave is essentially fully resolved. As a consequence, one obtains a proof of the general rule of thumb sometimes quoted in the context of spectral element methods: π modes per wavelength are needed to resolve a wave. Key words. spectral element method, mass lumped scheme, numerical dispersion AMS subject classifications. 65N50, 65N15, 65N30, 35A40, 35J05 DOI. 10.1137/080724976

1. Introduction. The spectral element method [6] is a spectrally accurate algorithm for solving a wide variety of partial differential equations on unstructured grids. The computational domain is typically broken into quadrilateral elements, within each of which the variables are approximated by high degree polynomials. A set of discrete equations for the coefficients is derived using a weak form of the problem in which the integrals are approximated using a quadrature rule. In particular, if the Gauss– Legendre–Lobatto quadrature rule is chosen in conjunction with a Lagrange basis for the approximation based at the nodes of the Gauss–Legendre–Lobatto rule, then the resulting mass matrix is diagonal in certain cases. For this reason, the same approach is sometimes described as the Gauss-point mass lumped finite element scheme [8, 10] in the finite element literature. The fact that the mass matrix is diagonal means that the method is particularly attractive for the efficient numerical simulation of wave phenomena, if used in conjunction with an explicit time-stepping scheme. The approach has found widespread application in a variety of areas ranging from acoustical waves [23], hydrostatic fluid flow [18], tumor angiogenesis [30], reaction-diffusion problems [24], edge finite element approximation of Maxwell’s equations [7], and seismic wave propagation [20, 21, 22]. A priori error estimates and the asymptotic convergence properties of the method have been well-studied in the literature; see, e.g., [6] and the references therein. Zampieri and Pavarino [32, 33] consider the spectral element method applied to the ∗ Received by the editors May 23, 2008; accepted for publication (in revised form) September 17, 2009; published electronically December 16, 2009. http://www.siam.org/journals/sinum/47-5/72497.html † Department of Mathematics, University of Strathclyde, 26 Richmond Street, Glasgow G1 1XH, Scotland ([email protected], [email protected]). The first author’s research was supported by the Engineering and Physical Sciences Research Council under grant EP/E040993/1. The second author’s research was supported by COMSATS Institute of Information Technology, Pakistan through a research studentship.

3910

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3911

acoustic wave equation used in conjunction with the Newmark scheme for advancing the solution in time. The stability of the scheme is established along with the dependence of the maximal admissible time-step on the order of the elements, and the method is shown to converge as the time-step is reduced and the order of the spectral elements increased. Related results will also be found in [31]. One of the key features demanded of a numerical method for wave propagation is the ability to accurately reproduce the dispersive and dissipative characteristic of the underlying problem. Indeed, some numerical schemes have been derived purely from a desire to minimize the numerical dispersion; see, e.g., [5, 10, 17]. The dispersive properties of many finite difference and low-order finite element schemes have been studied extensively [8, 9, 10, 14, 15, 16, 25, 29]. However, it was only relatively recently [1] that a complete, sharp analysis of the dispersive behavior of conforming finite element methods was given for finite elements of any order p as the mesh-size h is reduced and, as the order of the method p is increased, on a fixed mesh. The analysis of [1] was carried out in the context of the Helmholtz equation, but similar results were subsequently obtained for discontinuous Galerkin finite element methods [2, 4] and for edge element approximation of Maxwell’s equations [3]. The notion of increasing the order of the method on a mesh of fixed size has been proposed by a number of sources. For instance, Stanescu, Kopriva, and Hussaini [28] study the amplitude and phase errors for discontinuous spectral element methods applied to wave propagation problems with attention on the situation where the element size is larger than the wavelength. Despite the widespread usage of the spectral element method for computational wave propagation loc. cit., there seems to be no detailed study of the dispersive and dissipative properties of the scheme in the case of (a) arbitrary order approximation p as the mesh-size h is reduced or (b) as the order of approximation p is increased on a fixed mesh of size h. The goal of the present work is to study this issue in detail and to provide both a qualitative description of the nature of the dispersive and dissipative behavior of the scheme, along with precise quantitative statements of the accuracy, in the case of both (a) and (b). The approach is broadly similar to one used by [1] in the analysis of the conforming finite element method, where the key step in the analysis was the derivation of an explicit expression for the dispersion relation for the numerical scheme. The heart of the present analysis lies in the derivation of the corresponding discrete dispersion relation for the spectral element method. However, the discrete dispersion relation for the finite element case obtained in [1] took the form of a rational function expressed in terms of Pad´e approximants for the tangent and cotangent functions. In contrast, in the case of the spectral element method considered here, the discrete dispersion relation again assumes the form of a rational function, but this is no longer related to Pad´e approximants. Whilst this fact complicates the issue, it is nevertheless possible to carry out a fairly complete analysis. The main conclusions and comparisons with the finite element case for approximation of a wave of frequency ω are as follows: 1. For fixed order of approximation p, as the mesh-size tends to zero the following hold: (a) both the finite element and spectral element schemes give the same O(ωh)2p+1 accuracy for the phase error; (b) the constant appearing in front of the leading order term for the error in the spectral element case is −1/p times that for the finite element case.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3912

MARK AINSWORTH AND HAFIZ ABDUL WAJID

2. For a fixed mesh of size h with ωh  1, as the order of the scheme p is increased, the following hold: (a) for p = O(1), the phase error of the finite element scheme is O(1), whereas the spectral element scheme is of order 12 (ωh)2 ; (b) for ωh − o(ωh)1/3 ≤ 2p + 1 ≤ ωh + o(ωh)1/3 , the phase accuracy of both schemes undergoes a sharp transition between an error of order O(1) and a situation in which both schemes provide an essentially fully resolved numerical wave; (c) for 2p + 1  ωh, the error obtained with the spectral element scheme performs (ωh/4p)2 times better than that of the finite element scheme: (p)

ESE ≈



ωh 4p

2

(p)

EF E .

The first conclusion agrees with the practical observation that the Gauss-point mass lumped scheme (i.e., spectral element method) tends to exhibit phase lag, whereas the (consistent) finite element scheme tends to exhibit phase lead. More interestingly, the conclusion shows that the absolute accuracy of the spectral element scheme is 1/p times better than that of the finite element scheme despite the use of numerical integration. The second conclusion means that in the unresolved regime, the spectral element method behaves rather erratically (in numerical examples one sees drastic overshooting and undershooting of the true wave) in contrast with the less erratic (but still very poor) behavior of the finite element scheme. Both schemes exhibit a sharp transition whereby the true wave is essentially fully resolved by increasing the order from p to roughly order p + 1 or p + 2. The transition corresponding to the stage where the spectral element scheme essentially provides full resolution occurs when the order p, the mesh-size h, and the frequency ω are related by 2p + 1 ≈ ωh. A mesh of size h corresponds to there being 2π/ωh elements per wavelength. Each element in a pth order scheme involves p + 1 Gauss–Legendre–Lobatto points in each direction, or on average p + 1/2 degrees of freedom per element in each direction. Consequently, at the point when full resolution occurs, the scheme requires roughly (p + 1/2) × 2π/ωh = π degrees of freedom per wavelength. This agrees with the general rule of thumb sometimes quoted in the context of spectral element methods: π modes per wavelength are needed to resolve a wave, and the arguments presented here can be regarded as a rigorous proof of that fact. Finally, it is worth recalling the fundamental fact that at least two degrees of freedom per wavelength are needed for any scheme to resolve a wave. Consequently, the ability of the spectral element (and finite element) method to resolve a wave with π degrees of freedom per wavelength is close to optimal and perhaps helps to explain the popularity of such methods for computational wave propagation. The remainder of the paper is organized as follows. We start by developing a unidimensional model problem in section 2. In sections 3 and 4, discrete dispersion relations are derived for linear and higher order approximating elements for the model problem. Moreover, numerical results obtained with both spectral element and finite element method are shown. The remaining sections contain proofs of the results.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3913

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

2. Model problem and its discretization. Consider the one-dimensional model problem of wave propagation ∂2u ∂2u − 2 = 0, ∂t2 ∂x

(2.1)

x ∈ (0, 1),

t > 0,

subject to the boundary conditions u(0, t) = 0

∂u ∂u (1, t) + (1, t) = 0 ∂x ∂t

and

for t > 0,

and initial conditions u(x, 0) = u0 (x)

and

∂u (x, 0) = v0 (x) ∂t

for x ∈ (0, 1).

1 The variational formulation of the above problem is as follows: find u(t) ∈ HE (0, 1), t > 0, such that along with the initial and boundary conditions

d2 ∂u (1, t)v(1) = 0 (u, v) + (u , v  ) + dt2 ∂t

1 ∀v ∈ HE (0, 1),

t > 0,

1 where (·, ·) denotes the L2 -inner product on (0, 1) and HE (0, 1) = {v ∈ H 1 (0, 1) : 1 v(0) = 0}, where H (0, 1) is the usual Sobolev space. We construct a semidiscretization in space by introducing a partition Gh = {kh, k = 0, 1, . . . , K} of the interval 1 (0, 1) be the cor(0, 1) into K subintervals of equal length h = 1/K. Let Vhp ⊂ HE responding space of continuous piecewise polynomials of degree p ∈ N defined on Gh . We seek an approximate solution uhp ∈ Vhp such that

 ∂uhp  d2  (1, t)vhp (1) = 0 + (uhp , vhp ) + uhp , vhp 2 dt ∂t

∀vhp ∈ Vhp .

We can construct a basis for Vhp in terms of basis functions {Θ }p=0 defined on a reference element (−1, 1) as follows. Let −1 = ξ0 < ξ1 < · · · < ξp = 1 be distinct nodes on [−1, 1] and define Θ ∈ Pp (−1, 1) by the conditions Θ (ξm ) = δm . The N corresponding global basis functions for the entire mesh are denoted by {θi }i=1 and satisfy θi (1) = 0 for i = 1, 2, . . . , N − 1; then uhp can be written as uhp (x, t) =

N 

αi (t)θi (x),

x ∈ (0, 1),

t ≥ 0,

i=1

where αi are smooth complex-valued functions satisfying (2.2)

 N     dαi d2 αi (t)δiN δjN |θN (1)|2 = 0 (θi , θj ) 2 (t) + αi (t) θi , θj + dt dt i=1

for all j = 1, 2, . . . , N. By letting α denote the vector whose components are αi , i = 1, 2, . . . , N, (2.2) can be written in matrix form as (2.3)

M

d2 α d α + K α = 0, +C 2 dt dt

where, for i, j = 1, 2, . . . , N ,   (2.4) Kij = θi , θj , Mij = (θi , θj ) ,

and Cij = δiN δjN |θN (1)|2 .

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3914

MARK AINSWORTH AND HAFIZ ABDUL WAJID

We may define a fully discrete scheme by discretizing the temporal derivative using centered differences with step-size Δt > 0 to give     Δt Δt 2 C α n+1 = [2M − (Δt) K] C α n−1 , n = 1, 2, 3, . . . , αn − M − (2.5) M + 2 2 where α n ≈ α (nΔt). The matrices appearing in (2.4) are the consistent stiffness and mass matrices, which can be assembled in the usual fashion from the corresponding matrices defined on the reference element as follows: 1 1 dΘ (s) dΘm (s) ˆ m = h ˆ m = 2 ds and M Θ (s)Θm (s)ds (2.6) K h −1 ds ds 2 −1 for all , m ∈ {0, 1, 2, . . . , p}. Observe that at each time-step it is necessary to invert the matrix M + Δt 2 C. In practice the mass matrix M is often replaced by a lumped

mass matrix M [19], which means that each time-step involves the inversion of the

+ Δt C. The lumped mass matrix is obtained by employing the diagonal matrix M 2 spectral element method [6]. To derive the spectral element scheme we replace the integrals appearing in (2.6) with a numerical quadrature rule. The (p + 1)-point Gauss–Lobatto quadrature rule is defined by 1 p−1  2 [f (−1) + f (1)], f (s)ds ≈ Q(p) (f ) = w f (ξ ) + (2.7) p(p + 1) −1 =1

where are taken to be the zeros of Lp , and Lp is the pth order Legendre polynomial [13] with the weights given by [27, eqs. (4.10–26)] {ξ }p−1 =1

(2.8)

w =

2 p(p − 1)[Lp−1 (ξ )]2

∀ ∈ {1, 2, 3, . . . , p − 1}.

The quadrature rule is exact for polynomials of degree at most 2p − 1. Now, using (2.7) to approximate the integrals, the elemental stiffness and mass matrices appearing in (2.6) are replaced by the following forms:  dΘ (ξ ) dΘm (ξ )  ˆ m ≈ h ˆ m ≈ 2 w and M Θ (ξ )Θm (ξ )w K h ds ds 2 p

p

=0

=0

 dΘm for all , m ∈ {0, 1, 2, . . . , p}. As the product dΘ ds ds ∈ P2p−2 , the element stiffness matrix is not affected by the use of the quadrature rule, whereas the mass matrix will be different. A key idea in the spectral element method is that by choosing the nodes used to construct the basis functions to coincide with the nodes used in the quadrature rule, the elemental mass matrix becomes diagonal:

ˆ m ≈ h Q(p) (Θ Θm ) = h δm w ∀ , m ∈ {0, 1, 2, . . . , p}, M 2 2 Δt

and consequently, the matrix M + 2 C will be diagonal. The corresponding fully discrete scheme takes the form    

− Δt C α

+ Δt C α

− (Δt)2 K] M n+1 = [2M n−1 , n = 1, 2, 3, . . . , αn − M 2 2 where, thanks to the use of the quadrature rule, it is now necessary only to invert a diagonal matrix at each time-step.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3915

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3. Dispersive and dissipative behavior in space. In order to focus on the spatial discretization, we consider a time harmonic solution of the form u(x, t) = U (x)e−iωt for ω ∈ R fixed where the spatial component U satisfies −U  (x) − ω 2 U (x) = 0,

(3.1)

x ∈ (0, 1),

with boundary conditions U (0) = 1, U  (1) − iωU (1) = 0. Here, we choose nonhomogeneous Dirichlet data.

2

2 Exact real wave Numerical real wave obtained using the finite element method Numerical real wave obtained using the spectral element method

1.5

1

0.5

0.5 U

1

U

Exact real wave Numerical real wave obtained using the finite element method Numerical real wave obtained using the spectral element method Numerical real wave obtained using the weighted averaging method

1.5

0

0

−0.5

−0.5

−1

−1

−1.5

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

1

−1.5

0

0.1

0.2

0.3

0.4

0.5 x

(a)

0.6

0.7

0.8

0.9

1

(b)

Fig. 3.1. Numerical approximations of the solution to (3.1) obtained for ω = 20 with both finite element and spectral element methods using (a) 20 and (b) 40 linear elements.

The analytical solution of the above problem is U (x) = eiωx , so that both the real and imaginary parts oscillate. In Figure 3.1, we present the numerical approximations obtained using piecewise linear elements with both finite element and spectral element methods. The phase lead and phase lag of the numerical approximations are clearly visible. Furthermore, it is evident that the phase lead occurs when the solution is approximated using the finite element method (full integration), whereas phase lag corresponds to the approximation obtained with the spectral element method (numerical quadrature). Whilst it is possible to eradicate numerical dispersion and dissipation due to temporal discretization completely, by using, for example, an exponential integrator, this is much more problematic in the case of spatial discretization. In the one-dimensional case it is possible to modify the scheme to obtain a nondispersive approximation in space, but this is not possible in higher numbers of spatial dimensions [5]. 3.1. Numerical dispersion and dissipation. Let us study the dispersive and dissipative behavior of the finite element and spectral element schemes in more detail. Let Vh1 denote the set of continuous piecewise linear functions with nodes located at the nodes of the grid Gh . We seek an approximation Uh1 ∈ Vh1 of problem (3.1) satisfying Uh1 (0) = 1 such that (3.2)

B(Uh1 , vh1 ) = 0

1 ∀vh1 ∈ Vh1 ∩ HE (0, 1),

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3916

MARK AINSWORTH AND HAFIZ ABDUL WAJID

  where B(Uh1 , vh1 ) = (Uh1 , vh1 ) − ω 2 (Uh1 , vh1 ) − iωUh1(1)vh1 (1). We can write Uh1 in the form

Uh1 (x) =

N 

αi θi (x),

i=0

where θi (jh) = δij for all i, j = 0, 1, . . . , N are the usual piecewise linear hat functions associated with the nodes in the grid. Substituting Uh1 into (3.2) and taking vh1 = θj for j ∈ {0, 1, . . . , N }, we obtain α0 = 1, and N 

αi B (θi , θj ) = 0,

j = 1, 2, . . . , N.

i=0

Furthermore, since the mesh is uniform, we may express B(θi , θj ) in terms of B(θ0 , θ1 ) and B(θ1 , θ1 ) to obtain that the above conditions are equivalent to

(3.3)

α0 = 1, αj−1 B (θ0 , θ1 ) + αj B (θ1 , θ1 ) + αj+1 B (θ0 , θ1 ) = 0,   1 B (θ1 , θ1 ) − iω = 0. αN −1 B (θ0 , θ1 ) + αN 2

j = 1, 2, . . . , N − 1,

It is straightforward to see that the solution of the system (3.3) is given by αj =

[B(θ0 , θ1 ) cos(N − j)μ(1) h sin μ(1) h + iω sin(N − j)μ(1) h] [B(θ0 , θ1 ) cos N μ(1) h sin μ(1) h + iω sin N μ(1) h]

∀j = 0, 1, 2, . . . , N,

where μ(1) ∈ C is given by (3.4)

cos μ(1) h = −

1 B (θ1 , θ1 ) . 2 B (θ0 , θ1 )

All of the foregoing arguments apply equally well to the spectral element scheme leading to the same expression for the coefficients αj with the bilinear form B(·, ·) ˜ ·), where B(·, ˜ ·) denotes the bilinear form for the spectral element replaced by B(·, method obtained using reduced integration. We expect that μ(1) h → ωh as h → 0, corresponding to the fact that the frequency of the discrete approximation approaches the frequency of the true solution as the grid Gh is refined. For finite h > 0, the difference μ(1) h − ωh gives a measure of the dispersive and dissipative behavior of the numerical scheme. We can calculate μ(1) h explicitly for the above schemes as follows. The first order basis functions on the reference element are given by Θ0 = (1 − s)/2 and Θ1 = (1 + s)/2 for all s ∈ [−1, 1]. Hence, by applying a change of variable, we obtain 2 1   2κ2 + 3 B(θ0 , θ1 ) = (Θ0 Θ1 − κ2 Θ0 Θ1 )ds = − h −1 3h and 4 B(θ1 , θ1 ) = h



1

−1

2 2 (Θ2 1 − κ Θ1 )ds = −

2(4κ2 − 3) , 3h

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3917

where κ = ωh/2. Similarly,  1  2   1  2 ˜ B(θ0 , θ1 ) = Θ0 (ξ )Θ1 (ξ ) − κ Θ0 (ξ )Θ1 (ξ ) w = − h h =0

and

 1   4 2(1 − 2κ2 ) 2 2 ˜ 1 , θ1 ) = B(θ Θ2 1 (ξ ) − κ Θ1 (ξ ) w = h h =0

with ξ0 = −1, ξ1 = 1, and w0 = w1 = 1. Consequently, we obtain   3 − 4κ2 (3.5) μ(1) h = cos−1 3 + 2κ2 for the finite element scheme and   μ ˜ (1) h = cos−1 1 − 2κ2

(3.6)

for the spectral element scheme. If κ  1, then the discrete dispersion relations (3.5) and (3.6) may be expressed as a series in ωh to give the well-known results [29, eq. (41)]: μ(1) h − ωh = −

(ωh)3 + ··· 24

and μ ˜(1) h − ωh =

(ωh)3 + ··· . 24

Note that the leading terms in the phase differences in these expressions are identical in magnitude but have opposite signs. These expressions confirm the behavior observed in Figure 3.1, where it was observed that the finite element scheme exhibits phase lead whilst the spectral element scheme exhibits phase lag of equal magnitude. It is interesting to carry out a similar analysis in the case of order p elements on a uniform mesh of size h in the limit as ωh → 0. The following expressions are obtained for the dispersion errors: μ

(p)

 2 2p+1 p! (ωh) 1 + ··· h − ωh = − 2 (2p)! 2p + 1

and μ ˜(p) h − ωh =

 2 2p+1 p! (ωh) 1 + ··· 2p (2p)! 2p + 1

(see Theorem 4.2 which we will prove later). Observe that, as before, the signs of the leading terms in the error are opposite, but with the spectral element scheme the magnitude of the error is 1/p times that of the finite element scheme in the limit ωh → 0. We shall also investigate the nature of the convergence behavior of the schemes as the order p → ∞ on a fixed mesh. In Figure 3.2, we present the numerical approximations obtained on a fixed mesh of size h = 1 with frequency ω = 80 for orders of approximation p = 39 to p = 44. It is observed that for orders p ≤ 41, both finite element and spectral element schemes fail to resolve the wave, with the spectral element scheme exhibiting wild overshoots and undershoots. However, once

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3918

MARK AINSWORTH AND HAFIZ ABDUL WAJID

2

2 Exact real wave Numerical real wave obtained using finite element method Numerical real wave obtained using spectral element method

1

1

0.5

0.5

0

0

−0.5

−0.5

−1

−1

−1.5

−1.5

−2

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

Exact real wave Numerical real wave obtained using finite element method Numerical real wave obtained using spectral element method

1.5

U

U

1.5

−2

1

0

0.1

0.2

0.3

(a) p = 39

0.6

0.7

0.8

0.9

1

2 Exact real wave Numerical real wave obtained using finite element method Numerical real wave obtained using spectral element method

1.5

Exact real wave Numerical real wave obtained using finite element method Numerical real wave obtained using spectral element method

1.5

1

1

0.5

0.5

U

U

0.5 x

(b) p = 41

2

0

−0.5

0

−0.5

−1

−1

−1.5

−1.5

−2

0.4

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

0.8

0.9

−2

1

0

0.1

0.2

0.3

(c) p = 42

0.4

0.5 x

0.6

0.7

0.8

0.9

1

(d) p = 43 Case: ω h = 80 ω h >> 1

4

10

2 Exact real wave Numerical real wave obtained using finite element method Numerical real wave obtained using spectral element method

1.5

2

10

0

10

1 |R(p)(ω h)−cos(ω h)|

−2

U

0.5

0

−0.5

10

p=39 p=44

−4

10

−6

10

−8

10 −1

−10

10

−1.5

10

−2

10

−12

Error in the discrete dispersion relation for finite elements Error in the discrete dispersion relation for spectral elements

−14

0

0.1

0.2

0.3

0.4

0.5 x

0.6

0.7

(e) p = 44

0.8

0.9

1

0

10

20

30 p

40

50

60

(f) Error comparison

Fig. 3.2. Numerical approximations of the solution to (3.1) obtained for p = 39 to p = 44 with ωh = 80. The error in the discrete dispersion relations of both waves is compared.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3919

the order reaches p = 42, there is a dramatic improvement in the resolution of both schemes, and for p ≥ 43, both schemes provide an essentially exact approximation of the wave. We see that the convergence behavior of the higher order schemes, whereby h is fixed and p is increased, is quite different from that of fixed order approximation. In particular, over a very narrow range of p (in this case from p = 39 to p = 44), the approximate solution changes from being little more than garbage to providing an essentially exact representation of the true wave. Figure 3.2(f) sheds some light on the nature of this dramatic change in the qualitative behavior of the schemes over a very narrow range. In particular, we see that μ ˜(p) h − ωh mirrors precisely the kind of sharp transition around p ≈ ωh/2 seen in the approximate waves. Moreover, for p  ωh, it is seen that μ ˜(p) h for the spectral element scheme is of completely the wrong magnitude (compared with ωh), thereby accounting for the erratic overshoots and undershoots in the unresolved regime. 4. Higher order discrete dispersion relation. We now derive and study the exact discrete dispersion relation for elements of arbitrary order. Our objectives are twofold. First, we wish to compare the phase accuracy of finite element and spectral element schemes of fixed, but arbitrary, order p as the mesh-size becomes small. In particular, we shall show that the superiority of the spectral element scheme observed earlier is maintained for all orders p ≥ 2. Second, we have the more ambitious goal of explaining the dramatic behavior of the convergence of the schemes as the order is increased on a fixed mesh. The key to both of these analyses will be an explicit expression for the discrete dispersion relation for the spectral element scheme. Such expressions have been obtained (also above) for relatively low orders p = 1, 2, whereas our new result will be valid for arbitrary order p ∈ N. Let Vhp denote the set of continuous piecewise polynomials of degree p ≥ 2 on the grid Gh . In order to obtain the discrete dispersion relation for higher order elements, following [1], we introduce (p) basis functions {θ˜i }N i=0 ∈ Vhp satisfying the conditions (4.1)

(p) ˜ θ˜(p) , vhp ) = 0 θ˜i (jh) = δij , i, j ∈ Gh and B( i

 ∀vhp ∈ Vhp ,

(p)

 where Vhp = {vhp ∈ Vhp : vhp (jh) = 0, j ∈ Gh } with {θi }N i=0 defined similarly using ˜ ·). We seek a solution Uhp ∈ Vhp of (3.1) of the form B(·, ·) in place of B(·,

Uhp (x) =

N 

(p)

αi θ˜i (x)

i=0

˜ hp , vhp ) = 0 for all vhp ∈ Vhp ∩ H 1 (0, 1). Now, satisfying Uhp (0) = 1 such that B(U E following the arguments used earlier for linear elements, we arrive at the expressions (4.2)

cos μ ˜ (p) h = −

˜ θ˜(p) , θ˜(p) ) 1 B( 1 1 2 B( ˜ θ˜(p) , θ˜(p) ) 0

1

(p)

(p)

0

1

and (4.3)

cos μ(p) h = −

1 B(θ1 , θ1 ) . 2 B(θ(p) , θ(p) )

The analysis can be extended to the multidimensional problems on tensor product meshes in a way similar to section 2.3 of [1].

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3920

MARK AINSWORTH AND HAFIZ ABDUL WAJID

∞ Theorem 4.1. Let κ > 0 be given. Define sequences {ap }∞ p=1 and {bp }p=1 recursively by the rule

(4.4)

ap+1 = −

2p + 1 2p + 1 bp + ap−1 and bp+1 = ap + bp−1 κ κ

for p ∈ N with a0 = 1, a1 = 1, b0 = 0, and b1 = κ1 . Then, the discrete dispersion relation for the spectral element method is given, for p ∈ N, by cos μ ˜(p) h = R(p) (2κ), where   ap (κbp−1 + pap ) + bp (κap−1 − pbp ) (4.5) R(p) (2κ) = (−1)p . ap (κbp−1 + pap ) − bp (κap−1 − pbp ) The above expression is a rational function of κ since both the series defined in Lemma 5.2 are polynomials in powers of κ−1 , whilst the numerators and denominators are polynomials of degrees 2p and 2p − 2, respectively, in κ. Hence, the degree of R(p) (2κ) is [2p/2p − 2] for all p ∈ N. Consider the first order approximation (p = 1); then using (4.4) and (4.5), we find that cos μ ˜ (1) h = R(1) (2κ) = 1 − 2κ2 in agreement with our expression (3.6). Table 4.1 gives closed form expressions for R(p) (ωh) for p = 1, 2, 3 along with the leading term in the error for ωh  1. Table 4.1 The discrete dispersion relation R(p) (ωh) = cos μ ˜(p) h for order p approximation given in Theorem 4.1. We also indicate the leading term in the series expansion for the error when ωh  1 (see Theorem 4.2). Order p

R(p) (ωh) 1−

1

(ωh)2 2

μ ˜(p) h − ωh (ωh)3 24

2

(ωh)4 −22(ωh)2 +48 2((ωh)2 +24)

(ωh)5 2880

3

−(ωh)6 +92(ωh)4 −1680(ωh)2 +3600 2((ωh)4 +60(ωh)2 +1800)

(ωh)7 604800

4.1. Accuracy at small wave numbers. The following theorem (proved in section 5) gives the leading term for the error in the discrete dispersion relation, when ωh  1, for arbitrary order p ∈ N. Theorem 4.2. Let p ∈ N. Then, the error in the discrete dispersion relation for the spectral element method is given by (4.6)

cos μ ˜(p) h − cos ωh = −

 2 1 (ωh)2p+2 p! 2p+4 + O (ωh) 2p (2p)! 2p + 1

or, if ωh is sufficiently small, (4.7)

μ ˜

(p)

 2 1 (ωh)2p+1 p! 2p+3 + O (ωh) h − ωh = . 2p (2p)! 2p + 1

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3921

The leading terms in the series expansion for the error obtained using the spectral element method all have positive signs, thereby explaining the phase lag observed in the previous section. This result in the case p = 1 agrees with the expression (41) given in Thompson and Pinsky [29]. The case of general order p approximation is not conducted by Thompson and Pinsky, who nevertheless conjectured that the leading term in the expressions should be of order (ωh)2p . The correctness of this conjecture is confirmed by Theorem 4.2. More interestingly, for higher orders the spectral element scheme provides p-times better accuracy as compared to the discrete dispersion relation obtained with finite element scheme [1], where one finds that the corresponding result for the finite element scheme is μ(p) h − ωh = −

(4.8)

 2 p! (ωh)2p+1 1 + O (ωh)2p+3 . 2 (2p)! 2p + 1

It is evident from expression (4.7) that the error for the spectral element scheme has an additional multiple of 1/p for elements of order p ∈ N when compared with (4.8), again with a sign change. 4.2. Accuracy at large wave numbers. We now consider the error estimates for high wave numbers, i.e., ωh is large even though h is small. The next theorem gives a full description of the behavior of the error for large ωh as the order of the approximation p is increased. ˜(p) h − cos ωh Theorem 4.3. Suppose that ωh  1. Then the error E (p) = cos μ in the discrete dispersion relation for the spectral element method passes through four distinct phases as the order p ∈ N is increased: 2 1. Constant magnitude phase: For p = O(1), E (p) ≈ (−1)p (ωh) 2 . 2. Oscillatory phase: For 1  2p + 1 < ωh − o(ωh)1/3 , E (p) oscillates and decays to O(1) as p is increased. 3. Transition zone: For ωh − o(ωh)1/3 < 2p + 1 < ωh + o(ωh)1/3 , the error E (p) oscillates without further decrease. 4. Superexponential decay: For 2p + 1 > ωh + o(ωh)1/3 , E (p) decreases at a superexponential rate: (4.9)

E

(p)

sin(ωh) ≈ 2



ωh 4p

2

f ( 1 − (ωh/(2p + 1))2 )p+1/2 ,

where f : w → (1 − w)/(1 + w) exp(2w) so that in the case where 2p + 1 > ωhe/2, (4.10)

E

(p)

sin(ωh) ≈ 2



ωh 4p

2 

ωhe 2(2p + 1)

2p+1 .

Figure 4.1 shows the behavior of the actual error obtained with the spectral element scheme for different values of ωh with increasing order p. It is seen that the behavior is consistent with the predictions of Theorem 4.3. It is interesting to compare the results in Theorem 4.3 with those obtained in [1] for the finite element scheme. Observe that if EFp E denotes the corresponding error for the finite element scheme, then using the above result and Theorem 3.3 of [1], we have (4.11)

E

(p)

 ≈

ωh 4p

2

(p)

EF E .

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3922

MARK AINSWORTH AND HAFIZ ABDUL WAJID

4

6

10

10

4

10

2

10

2

10

Exponential Decay Phase

0

10

|R(p)(ω h)−cos(ω h)|

|R(p)(ω h)−cos(ω h)|

0

Oscillating and Decaying Phase

−2

10

−4

10

10

−2

10

−4

10

−6

10 −6

10

−8

10 −8

10

−10

Constant Magnitude Phase

10

Transition Zone (Oscillating)

−10

10

−12

0

2

4

6

8

10 p

12

14

16

18

20

10

0

10

20

30

(a) ωh = 20

40

50 p

60

70

80

90

100

(b) ωh = 160

Fig. 4.1. Behavior of the error in the discrete dispersion relation for the spectral element scheme at high wave numbers ωh  1 as the order p is increased. The transition region between the oscillatory decaying phase and the superexponential decay of the error is indicated in each case (cf. Theorem 4.3) and occurs when 2p + 1 ≈ ωh.

Broadly speaking, this means the performance of the two methods is similar in terms of (a) the fact that there is a sharp transition to superexponential convergence, and (b) the transition occurs at the same threshold. The most significant difference in the behavior of the spectral element and finite element schemes occurs in the unresolved regime where 2p+1 < ωh−o(ωh)1/3 . For the finite element scheme, the error is of order 1, whereas for the spectral element scheme, the error is of order 12 (ωh)2 . This behavior was observed in Figure 3.2, where it was found that the spectral element scheme overshoots and undershoots erratically in this range. In summary, in agreement with [1], we recommend that the order p and the mesh-size h be chosen so that 2p + 1 > ωh + C(ωh)1/3 , where C is a fixed constant and can be taken as unity in practice. 4.3. Nonphysical propagating modes. In the foregoing analysis, we have shown that in certain regimes, e.g., ωh → 0 or p → ∞, the numerical scheme admits discrete counterparts to the propagating modes e±iωx of the continuous equation. For example, it was shown that the pth order scheme admits Bloch waves with frequency μ ˜(p) , where cos μ ˜(p) h = R(p) (ωh) with R(p) defined in Theorem 4.1. In the special case of second order elements, we find that μ ˜(2) should satisfy (4.12)

cos μ ˜(2) h =

(ωh)4 − 22(ωh)2 + 48 , 2((ωh)2 + 24)

and hence, for ωh  1, there exists a solution of the form (˜ μ(2) h)2 = (ωh)2 +

1 (ωh)6 + · · · . 1440

However, if one rearranges (4.12) and seeks solutions for ωh in terms of μ ˜(2) h, then it is found that there exist two families of solutions given by  (2) μ ˜ h μ ˜(2) h μ ˜(2) h ± 2 36 − 36 sin2 + sin4 . (ωh)2 = 12 − 2 sin2 2 2 2

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3923

Following the approach of Cohen [10] and expanding these expressions, we find that   μ ˜(2) h 1 μ ˜(2) h 1 (ωh)2 = 4 sin2 + sin4 + · · · ≈ (˜ μ(2) h)2 − (˜ μ(2) h)6 + · · · 2 3 2 1440 and (ωh)2 = 24 − 8 sin2

μ ˜(2) h + ··· . 2

The first set of solutions obviously corresponds to the discrete Bloch wave described above. The second family of propagating modes exists only for frequencies ω satisfying ω 2 ∼ 24/h2 . The fact that these frequencies are mesh dependent shows that this second set of solutions is an artifact of the discretization process and is sometimes referred to as spurious, parasitic, or nonphysical propagating modes [10]. The situation persists for elements of all orders p ≥ 2. The fact that the leading ˜(p) term in the series for (ωh)2 corresponding to these frequencies is independent of μ (p) means that we can find the leading terms by formally inserting μ ˜ = 0 into the relation cos μ ˜(p) h = R(p) (ωh) to see that the spurious modes are of the form   ˜(p) h 2 2 μ (ωh) = β + O sin , 2 where β satisfies R(p) (β) = 1. It is not difficult to see that β = 0 is one solution (corresponding to the physical modes) whilst the remaining solutions correspond to the nonphysical modes. These nonphysical modes play little role in the behavior of the numerical scheme in the resolved regime where there exists a Bloch mode with frequency μ ˜ (p) exponentially close to ω. For instance, the exact solution of the problem −u − ω 2 u = f on (−1, 1) with appropriate boundary conditions takes the form u = ups + α1 eiωx + α2 e−iωx , where ups is a particular solution of the equation that depends only on the data f . Suppose we approximate the above equation using pth order elements on a mesh of K elements. This means we have Kp − 1 internal degrees of freedom plus two additional degrees of freedom corresponding to the endpoints of the domain. Conceptually, one could construct a basis for this discrete space consisting of standard hierarchic piecewise polynomials associated with the Kp − 1 interior degrees of freedom, plus an additional two nonlocal modes given by the discrete Bloch waves corresponding to the complementary solutions e±iωx . The internal modes provide an approximation to the particular solution ups and the nonlocal modes approximate the complementary function α1 eiωx + α2 e−iωx , whilst the nonphysical modes play no role and will not be activated by the numerical scheme in the resolved regime. 5. Proofs of the results. This section provides the proofs of the results for the error in the discrete dispersion relation for the spectral element method. ˆ ·) denote the bilinear 5.1. Basic polynomials. For u, v ∈ H 1 (−1, 1), let B(·, form (5.1)

ˆ v) = Q(p) (u v  ) − κ2 Q(p) (uv), B(u,

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3924

MARK AINSWORTH AND HAFIZ ABDUL WAJID

where Q(p) is the quadrature rule (2.7) and κ > 0 is a fixed constant. Polynomials Φp , Ψp ∈ Pp are defined for p ∈ N by (5.2)

ˆ p , v) = 0 Φp (1) = 1, Φp (−1) = (−1)p+1 : B(Φ

∀v ∈ Pp ∩ H10 (−1, 1)

ˆ p , v) = 0 B(Ψ

∀v ∈ Pp ∩ H10 (−1, 1).

and (5.3)

Ψp (1) = 1, Ψp (−1) = (−1)p :

These polynomials will play a central role in the derivation of the discrete dispersion relation. The following theorem provides explicit closed forms for the expressions ˆ p , Ψp ) and B(Φ ˆ p , Φp ), which we shall need later. B(Ψ Theorem 5.1. Let p = 2, 3, 4, . . . ; then ˆ p , Φp ) = −2κ ap B(Φ bp

(5.4) and

ˆ p , Ψp ) = −2κ (p + 1)ap−1 + pap+1 , B(Ψ (p + 1)bp−1 + pbp+1

(5.5)

where {ap } and {bp } are defined in Theorem 4.1. Proof. We begin by considering Φp . If p is odd, then Φp (1) = 1 = Φp (−1), and so Φp is an even function which implies Φp ∈ Pp−1 . Similarly, if p is even, then Φp (−1) = −1 = −Φp (1), and so Φp is an odd function which implies Φp ∈ Pp−1 . Hence, Φp ∈ Pp−1 for all p ∈ N. Let v ∈ Pp ∩ H10 (−1, 1); then vΦp and v  Φp ∈ P2p−1 , so that (5.1) integrates the function exactly when u = Φp , and we obtain 1  p   ˆ p , v) = Φ v − κ2 Φp v dx = 0 ∀v ∈ Pp ∩ H10 (−1, 1). B(Φ −1

p

Hence, we see that Φ = Φpe or Φp = Φpo , where Φpe and Φpo are polynomials analyzed in [1, eqs. (4.1)–(4.2)], and from [1, eqs. (4.13)–(4.14)], again using the fact that the quadrature rule is exact in this case, we have 1  p p  ap p p ˆ Φ Φ − κ2 Φp Φp dx = −2κ , B(Φ , Φ ) = bp −1 which proves the assertion for Φp . We now consider Ψp , and for the remainder of the proof superscripts will be omitted since no confusion is likely to arise. Since Ψ and v ∈ Pp and Ψ and v  ∈ Pp−1 , we can use the fact that the quadrature rule is exact for P2p−1 and integration by parts to obtain 1 (5.6) Q(p) (Ψ v  ) = Ψ v  dx = [vΨ ]1−1 − Q(p) (Ψ v) = −Q(p) (Ψ v) −1

since v(±1) = 0. Combining this result with (5.3) gives (5.7)

Q(p) (F v) =

p−1 

w F (ξ )v(ξ ) = 0

∀v ∈ Pp ∩ H10 (−1, 1),

=1

where F = (Ψ + κ2 Ψ) ∈ Pp . Fix J ∈ {1, 2, 3, . . . , p − 1}, and let v ∈ Pp ∩ H10 (−1, 1) be chosen such that v(ξ ) = δJ ; then (5.7) implies that F (ξJ ) = 0. Hence, F = 0 at

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3925

the zeros of Lp , where Lp is the Legendre polynomial of degree p. Hence, F takes the form F (x) = Ψ (x) + κ2 Ψ(x) = σxLp (x)

(5.8)

for some σ ∈ R. Now using the identity satisfied by Legendre polynomials xLp (x) =

 1  (p + 1)Lp−1 (x) + pLp+1 (x) , 2p + 1

obtained using equations (13)–(14) given in [11, section 10.10], (5.8) becomes  σ   (5.9) F (x) = pLp+1 (x) + (p + 1)Lp−1 (x) . 2p + 1 Define a polynomial Υn ∈ Pp by the rule n/2 

(5.10)

Υn (x) =

 j=0

1 − 2 κ

j+1

(2j+1)

Ln+1 (x).

It is elementary to verify that Υn (x) + κ2 Υn (x) = −Ln+1 (x). We may write Ψ(x) = αΥp (x) + βΥp−2 (x) for suitable constants α and β, and in addition     Ψ (x) + κ2 Ψ(x) = α Υp (x) + κ2 Υp (x) + β Υp−2 (x) + κ2 Υp−2 (x)   = − αLp+1 (x) + βLp−1 (x) . Comparing the last equation with (5.9), we are led to the conclusion that α = −σp/(2p + 1) and β = −σ(p + 1)/(2p + 1), and with these values, we obtain (5.11)

Ψ(x) = −

σ (pΥp (x) + (p + 1)Υp−2 (x)) . 2p + 1

Applying the boundary condition Ψ(1) = 1, we obtain σ = −(2p + 1)/η(1), provided that η(1) is nonzero, with η(1) = (p + 1)Υp−2 (1) + pΥp (1). Consequently, Ψ(x) = η(x)/η(1). We wish to obtain a closed form expression for ˆ B(Ψ, Ψ) = Q(p) (Ψ2 ) − κ2 Q(p) (Ψ2 ). The function Ψ2 ∈ P2p−2 is integrated exactly by the quadrature rule, so that 1 1 Q(p) (Ψ2 ) = Ψ2 dx = [Ψ Ψ]1−1 − ΨΨ dx = 2Ψ (1) − Q(p) (Ψ Ψ), −1

−1

and therefore, ˆ B(Ψ, Ψ) = 2Ψ (1) − Q(p) (Ψ Ψ) − κ2 Q(p) (Ψ2 ) = 2Ψ (1) − Q(p) (Ψ(Ψ + κ2 Ψ)). Equation (5.8) implies that Q(p) (Ψ(Ψ + κ2 Ψ)) = Q(p) (σxΨLp (x)), and therefore, (5.12)

2p + 1 (p) ˆ Q (xΨLp (x)). B(Ψ, Ψ) = 2Ψ (1) + η(1)

Now, using the quadrature rule (2.7), we obtain Q(p) (xΨLp (x)) =

p−1 

w Ψ(ξ )Lp (ξ )ξ − w0 Ψ(−1)Lp (−1) + wp Ψ(1)Lp (1),

=1

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3926

MARK AINSWORTH AND HAFIZ ABDUL WAJID

and then, since ξ for all = 1, 2, 3, . . . , p − 1 are the zeros of Lp , we get Q(p) (xΨLp (x)) = w0 (−1)p+1 Lp (−1) + wp Lp (1). Moreover, since w0 = wp = Q

(p)

(xΨLp (x))

2 p(p+1) ,

substituting Lp (±1) = (±1)(p+1) p(p + 1)/2 gives

= 2, which on insertion into (5.12) gives ˆ B(Ψ, Ψ) =

(5.13)

2 (η  (1) + 2p + 1) η(1)

after straightforward manipulations. Now, as in [1], using (5.10), the values of Υ and its derivatives at the boundary x = 1 are given by Υp (1) = ap+1 − 1, Υp−2 (1) = ap−1 − 1, Υp (1) = −bp+1 /κ, and Υp−2 (1) = −bp−1 /κ, where ap+1 , ap−1 , bp+1 , and bp−1 are the expressions obtained from the recurrence relation (4.4). A proof of this will be provided in Lemma 5.2 below. Consequently, the values of η(1) and η  (1) may be written in the form 1 η(1) = − ((p + 1)bp−1 + pbp+1 ) κ and η  (1) + 2p + 1 = (p + 1)ap−1 + pap+1 . Finally, inserting the above values into (5.13) and simplifying gives (5.5), which completes the proof. We now give closed form expressions and present some elementary properties of the coefficients ap and bp defined in (4.4). ∞ n Lemma 5.2. Let {an }∞ n=0 and {bn }n=0 be defined as in (4.4). Then κ an is a polynomial of degree n in κ, and n/2

(5.14)

an =

 (−1)k (n + 2k)! 1 , (2k)! (n − 2k)! (2κ)2k

n = 0, 1, 2, . . . ,

k=0

while κn bn is a polynomial of degree n − 1 in κ, and (n−1)/2

(5.15)

bn =



k=0

1 (−1)k (n + 2k + 1)! , (2k + 1)! (n − 2k − 1)! (2κ)2k+1

n = 1, 2, 3, . . . ,

with b0 = 0 and . denoting the integer part. Moreover, these series satisfy       πκ sin(κ − πn/2) cos(κ − πn/2) an Jn+1/2 (κ) (5.16) , = cos(κ − πn/2) − sin(κ − πn/2) bn 2 −Yn+1/2 (κ) where J and Y are cylindrical Bessel functions of the first and second kind, respectively [13]. Proof. It is elementary to prove (5.14) and (5.15) using mathematical induction. The statements regarding κn an and κn bn are then simple consequences of (5.14) and (5.15). The identity (5.16) follows at once from equation (4.12) in [1]. Note that a typographical error in [1] for (5.16) has been rectified here. ˆ p , Ψp ) Equations (5.4) and (5.5) provide compact representations for the terms B(Ψ p p p p ˆ ˆ and B(Φ , Φ ), respectively. We first consider B(Ψ , Ψ ) for even values of p. Using

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3927

Lemma 5.2, we see that κp ap and κp bp are polynomials of degrees at most p and p − 1 in κ, respectively, i.e., κp ap ∈ Pp , and κp bp ∈ Pp−1 . Consequently, it is not difficult ˆ p , Ψp ) to verify that κp ap+1 , κp ap−1 ∈ Pp and κp+1 bp−1 , κp+1 bp+1 ∈ Pp . Hence, B(Ψ is a rational function of degree [p + 2/p] in κ for even values of p. Similarly, for odd ˆ p , Ψp ) is a rational function of degree [p + 1/p − 1] in κ. Now, as values of p, B(Ψ p Φ ∈ Pp−1 for all p ∈ N, i.e., Φp is polynomial of degree p − 1 in x, the quadrature rule in bilinear form (5.1) is exact for Φp . Hence, from Theorem 4.2 in [1], it is clear that ˆ p , Φp ) is also a rational function of κ for even and odd values of p. As p = 2N B(Φ ˆ p , Φp ) is a rational function of for even values of p, from Theorem 4.2(2) in [1], B(Φ ˆ p , Φp ) is a rational function degree [p/p − 2] in κ. Similarly, for odd values of p, B(Φ of degree [p + 1/p − 1] in κ when 2N + 1 is replaced by p in Theorem 4.2(1) from [1]. ˆ p , Φp ) could be represented in terms of Moreover, in [1] it was shown that B(Φ Bessel functions as follows: ˆ p , Φp ) = −2κ Jp+1/2 (κ) sin κ − Yp+1/2 (κ) cos κ B(Φ Jp+1/2 (κ) cos κ + Yp+1/2 (κ) sin κ

with κ = mπ

and ˆ p , Φp ) = 2κ Jp+1/2 (κ) cos κ + Yp+1/2 (κ) sin κ B(Φ Jp+1/2 (κ) sin κ − Yp+1/2 (κ) cos κ

with κ = (m + 1/2)π.

ˆ p , Ψp ). The following result extends the previous results to B(Ψ Corollary 5.3. Let p = 2, 3, . . . ; then we have the following: 1. If p is even and κ = (m + 1/2)π for all m ∈ Z, ˆ p , Ψp ) B(Ψ (p + 1)(Jp−1/2 (κ) cot κ + Yp−1/2 (κ)) − p(Jp+3/2 (κ) cot κ + Yp+3/2 (κ)) . = 2κ (p + 1)(Jp−1/2 (κ) − Yp−1/2 (κ) cot κ) − p(Jp+3/2 (κ) − Yp+3/2 (κ) cot κ)

(5.17)

2. If p is odd and κ = mπ for all m ∈ Z, ˆ p , Ψp ) B(Ψ (p + 1)(Yp−1/2 (κ) cot κ − Jp−1/2 (κ)) − p(Yp+3/2 (κ) cot κ − Jp+3/2 (κ)) . = 2κ (p + 1)(Jp−1/2 (κ) cot κ + Yp−1/2 (κ)) − p(Jp+3/2 (κ) cot κ + Yp+3/2 (κ))

(5.18)

Proof. This corollary is proved separately for even and odd order polynomials. Consider first the case when p is even. Since the series {ap } and {bp } for nonnegative integers p satisfy identity (5.16), using (5.16) the values of ap−1 , ap+1 and bp−1 , bp+1 are given as follows:   πκ  ap−1 = Jp−1/2 (κ) cos κ + Yp−1/2 (κ) sin κ , 2   πκ  −Jp+3/2 (κ) cos κ − Yp+3/2 (κ) sin κ , ap+1 = 2   πκ  −Jp−1/2 (κ) sin κ + Yp−1/2 (κ) cos κ , bp−1 = 2 and

 bp+1 =

 πκ  Jp+3/2 (κ) sin κ − Yp+3/2 (κ) cos κ . 2

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3928

MARK AINSWORTH AND HAFIZ ABDUL WAJID

Now, inserting these values into (5.5) and rearranging gives (5.17), which completes the proof in the even case. Similarly, when p is odd, identity (5.16) together with expression (5.5) gives (5.18), which completes the proof. ˆ p , Φp ) was related to certain types of Pad´e approxiIn [1], it was shown that B(Φ ˆ p , Ψp ) is a Pad´e approximant. Nevertheless, mants. Here, it is not the case that B(Ψ we have the following result. Theorem 5.4. Let p ∈ N satisfy p ≥ 2. Then we have the following: 1. If p is an even integer, then i. if κ = (m + 1/2)π, m ∈ Z, (5.19) (p) EΨ

ˆ p , Ψp ) + 2κ tan κ = − 1 = B(Ψ 2

  2 (2κ)2p+2 p! 1 + O(κ2p+4 ); 1+ p (2p)! 2p + 1

ii. if κ = mπ, m ∈ Z, (5.20)

(p) EΦ

 ˆ p , Φp ) − 2κ cot κ = 2 = B(Φ

p! (2p)!

2

(2κ)2p + O(κ2p+2 ). 2p + 1

2. If p is an odd integer, then i. if κ = mπ, m ∈ Z, (5.21) (p) EΨ

  2 1 (2κ)2p p! p p ˆ + O(κ2p+2 ); = B(Ψ , Ψ ) − 2κ cot κ = −2 1 + p (2p)! 2p + 1

ii. if κ = (m + 1/2)π, m ∈ Z, (5.22)

(p)

ˆ p , Φp ) + 2κ tan κ = EΦ = B(Φ

 2 p! (2κ)2p+2 1 + O(κ2p+4 ). 2 (2p)! 2p + 1

Proof. The quadrature rule is exact for Φp in the bilinear form (5.1). Moreover the estimates (5.20) and (5.22) are the same as (4.15) and (4.16) in [1] for p = 2N + 1 and p = 2N , respectively. Therefore, we have  2 p! (2κ)2p (p) + O(κ2p+2 ) Eo(p) = EΦ = 2 (2p)! 2p + 1 and Ee(p)

=

(p) EΦ

 2 p! (2κ)2p+2 1 + O(κ2p+4 ). = 2 (2p)! 2p + 1

Now, to prove estimates (5.19) and (5.21) consider the case when p is even. Straightforward manipulations beginning with (5.17) give  −1 ˜ p+3/2 (κ) 1 − Q ˜ p+3/2 (κ) tan κ ˆ p , Ψp ) + 2κ tan κ = − 2κ Q , (5.23) B(Ψ cos2 κ ˜ p+3/2 (κ) is studied in ˜ p+3/2 (κ) = (p+1)Jp−1/2 (κ)−pJp+3/2 (κ) . The behavior of Q where Q (p+1)Yp−1/2 (κ)−pYp+3/2 (κ) the appendix, where the following estimate is proved in Lemma A.1:   2 1 (2κ)2p+1 1 p! p+3/2 ˜ (5.24) Q + ··· (κ) = 1+ 2 p (2p)! 2p + 1

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3929

for κ  1. With the aid of this estimate, we obtain (5.19) as required. The assertions concerning polynomials of odd order are proved starting with (5.18) in a similar fashion, and we have  −1 ˜ p+3/2 (κ) 1 + Q ˜ p+3/2 (κ) cot κ ˆ p , Ψp ) − 2κ cot κ = − 2κ Q . (5.25) B(Ψ sin2 κ Finally, inserting (5.24) into the above equation and performing ordinary manipulations gives (5.21) as required. 5.2. Proof of Theorem 4.1. The proof of Theorem 4.1 now follows using a virtually identical argument to the proof of Theorem 3.1 in [1]. For the sake of completeness, we give an outline here. (p) Proof. For x ∈ (0, h), the functions {θ˜i }N i=0 can be written in terms of the basic polynomials Φ and Ψ using s = 2x/h − 1: (−1)p p 1 (p) [Ψ (s) − Φp (s)] and θ˜1p (x) = [Ψp (s) + Φp (s)], θ˜0 (x) = 2 2 where the expressions for θ˜0p (x) and θ˜1p (x) take the correct values at the boundary points x = 0 and x = h. Moreover, we define V ∈ Pp ∩ H01 (−1, 1) by V (s) = vhp (x),  is supported on (0, h). Now, using the change of variable s = 2x/h−1, where vhp ∈ Vhp we obtain ˜ θ˜(p) , vhp ) = h−1 B(Ψ ˆ p + Φp , V ), B( 1 which vanishes because of the basis polynomials defined in (5.2)–(5.3). Consequently, (p) the orthogonality property is satisfied. Similar arguments hold in the case of θ˜0 . Since ωh = 2κ, the change of variable reveals that ˆ p + Φp , Ψp + Φp ). ˜ θ˜(p) , θ˜(p) ) = h−1 B(Ψ B( 1 1 ˆ p, ˜ θ˜(p) , θ˜(p) ) = h−1 [B(Ψ Furthermore, exploiting the parities of Ψ and Φ, we obtain B( 1 1 (p) (p) p p p p −1 p ˆ ˜ θ˜ , θ˜ ) = (−1) (2h) [B(Ψ ˆ , Φ )]. Similar arguments give B( , Ψp ) − Ψ ) + B(Φ 0 1 ˆ p , Φp )]. Substituting these results into (4.2) gives B(Φ cos μ ˜(p) h = (−1)p+1

(5.26)

ˆ p , Ψp ) + B(Φ ˆ p , Φp ) B(Ψ . ˆ p , Ψp ) − B(Φ ˆ p , Φp ) B(Ψ

Now using (5.4) and (5.5) together with the recurrence relation defined in (4.4) and performing ordinary manipulations gives (4.5) as required. 5.3. Proof of Theorem 4.2. We use Theorem 5.4 to prove Theorem 4.2 as follows. Proof. We first consider the case when p is even. We start with (5.19)–(5.20) and ˆ p , Φp ) in terms of E (p) and E (p) , respectively. Now, inserting ˆ p , Ψp ) and B(Φ write B(Ψ Ψ Φ these expressions into (5.26) and performing ordinary calculations, we obtain the following expression for the error in the discrete dispersion relation:  (5.27)

 ωh cos μ ˜(p) h − cos ωh sin ωh       −1 ωh ωh sin ωh  (p) (p) (p) (p) 2 2 = EΦ sin EΦ − EΨ . + EΨ cos 1+ 2 2 2ωh

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3930

MARK AINSWORTH AND HAFIZ ABDUL WAJID

2p+2

(p)

 ωh 2

(p) EΦ 2 2p+2 p! 2 (ωh) 1 p )[ (2p)! ] 2p+1 .

In particular, for small ωh  1, cos μ ˜(p) h − cos ωh ≈ (p)

(p)

2 + EΨ , where ( ωh 2 )

p! 2 (ωh) EΦ ≈ 12 [ (2p)! ] 2p+1 and EΨ ≈ − 12 (1 + Substituting these values, we get (4.6). We now consider the case when p is odd. In this case, starting with (5.21)–(5.22) and following the arguments used for the even order case, we arrive at

 (5.28)

 ωh cos μ ˜(p) h − cos ωh sin ωh       −1 sin ωh  (p) ωh ωh (p) (p) (p) 2 2 = EΨ sin . + EΦ cos 1+ EΨ − EΦ 2 2 2ωh

The proof for p odd follows exactly the same lines. Hence, it is once again evident that the leading term in the remainder is the same, regardless of the parity of the polynomial order p. It is also worth noting that unlike in [1] there were no parity dependent dominating terms in the expressions for the error. Furthermore, for small μ ˜(p) h, we can use the approximation cos μ ˜(p) h − cos ωh = −(˜ μ(p) h − ωh)ωh + · · · to obtain estimate (4.7). 5.4. Proof of Theorem 4.3. Proof. We begin with the regime where p  ωh. The series (5.14)–(5.15) defined in Lemma 5.2 gives ap = 1 −

(p + 2)! + O(κ−4 ) 8(p − 2)!κ2

and

bp =

(p + 3)! (p + 1)! − + O(κ−5 ). 2(p − 1)!κ 48(p − 3)!κ3

Now, substituting these values into (4.5) we obtain R(p) (2κ) = cos μ ˜(p) h ≈

 (−1)p+1  4 2p + 4p3 + p2 − p − 12κ2 6

after elementary calculations. Since κ = ωh/2  1, we can also obtain R(p) (2κ) = cos μ ˜(p) h ≈ (−1)p

(ωh)2 . 2

As E (p) = R(p) (2κ) − cos ωh, the error in the discrete dispersion relation for p  ωh is given by E (p) ≈ (−1)p

(ωh)2 2

∀p ∈ N such that p  ωh.

Now, for the rest of the proof we make use of the fact that (5.29)

(p)

EΨ (κ)

cos2 (ωh/2) ˜ p+3/2 (κ){1 − Q ˜ p+3/2 (κ) tan κ}−1 , = −Q ωh

˜ p+3/2 (κ) is studied which can be verified by rewriting (5.23), where the quotient Q in Theorem A.2 of the appendix. Since the quadrature rule is exact in bilinear form (5.1) for Φp , we can use [1, eq. (4.22)] with p = 2No , namely (5.30)

(p)

EΦ (κ)

sin2 (ωh/2) = −Qp+1/2 (κ){1 + Qp+1/2 (κ) cot κ}−1 , ωh

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3931

where the behavior of the quotient Qp+1/2 (κ) is studied in the appendix of [1]. First of all, consider the case when p is even with κ = ωh/2  1 fixed. The behavior of the error E (p) in different regimes will be determined by the behavior of ˜ p+3/2 (κ) or Qp+1/2 (κ). either Q We consider the preasymptotic regime where 2p + 1 < ωh − o(ωh)1/3 for both p+1/2 ˜ p+3/2 (κ). For p in this range, it is evident that in the preasymptotic (κ) and Q Q p+3/2 ˜ regime both Q (κ) and Qp+1/2 (κ) oscillate around unity but the error E (p) does not oscillate around unity because the denominator of (5.27) becomes very small. Therefore, the error E (p) starts from (ωh)2 /2 and decays by several orders of magnitude to O(1). ˜ p+3/2 (κ) have the same bounds for p in the transition Both Qp+1/2 (κ) and Q region, i.e., p in the transition region satisfies ωh− o(ωh)1/3 < 2p+ 1 < ωh+ o(ωh)1/3 . Moreover, with the spectral element method the term appearing in the denominator of (5.27) is oscillating around unity for ωh  1, which was not the case with the ˜ p+3/2 is also finite element method. In addition, from Theorem A.2, the quotient Q p+1/2 −1/3 decays algebraically at a rate O(p ), which is oscillating in this region whilst Q ˜ p+3/2 dominates Qp+1/2 in the error expression, proved in Theorem A.2 of [1]. Hence, Q and therefore, the error E (p) oscillates in the transition region. In the superexponential region, where p satisfies 2p + 1 > ωh + o(ωh)1/3 , the term appearing in the denominator of (5.27) is of O(1) for ωh  1. Hence, the error E (p) is ˜ p+3/2 and Qp+1/2 . Since in this region both dictated by the behavior of the sum of Q p+3/2 p+1/2 ˜ Q and Q decay at a superexponential rate, the error E (p) also decays at a ˜ p+3/2 and Qp+1/2 . superexponential rate as it is the sum of Q For the odd order case, the proof of the error expression (5.28) follows exactly the same arguments for the oscillatory, transition, and superexponential regions as in the even order case. ˜ m (κ). We now consider the behavior of the quotient Appendix. Behavior of Q Q (κ) for both cases, i.e., when κ  1 and κ  1. This quotient is defined by ˜m

(A.1)

˜ m (κ) = (2m − 1) Jm−2 (κ) − (2m − 3) Jm (κ) , Q (2m − 1) Ym−2 (κ) − (2m − 3) Ym (κ)

1 m = integer + , 2

and appeared in Theorem 5.4 in a form which we shall prove is equivalent to the above in the following lemma for small values of κ. ˜ m be defined as above. Then Lemma A.1. Let m = integer + 1/2, and let Q (A.2)

  2 m − 32 ! (2κ)2m−2 1 (2m − 1) + ··· Q (κ) = 2 (2m − 3) (2m − 3)! 2m − 2 ˜m

for κ  1.

Proof. Write m = n + 1/2, where n ∈ Z. For small κ, identities (8.440) and (8.399)2 of [13] give (A.3)

 κ n+1/2 2n+1 + ··· , Jn+1/2 (κ) = √ π(2n + 1)!! 2

while combining identities (8.465)1 and (8.399)3 along with identity (8.440) of [13] gives (A.4)

Yn+1/2 (κ) = −

(2n − 1)!!  κ −n−1/2 √ + ··· . 2 2n π

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3932

MARK AINSWORTH AND HAFIZ ABDUL WAJID Case: κ = 10

2

10

1

0

10

10

0

10

−1

Oscillatory Phase

−4

10

10

Exponential Decay Phase

−2

˜ m (κ)| |Q

˜ m (κ)| |Q

Exponential Decay Phase

−2

10

10

Oscillatory Phase −3

−6

10

−8

10

10

−4

−10

10

10

−5

−12

10

10

−6

10

Transition Zone (Algebraic Increase)

Transition Zone (Algebraic Decay)

−7

Transition Zone (Algebraic Decay)

Transition Zone (Algebraic Increase)

−14

10

10

Case: κ = 20

2

10

−16

2

4

6

8

10

12

14

16

18

10

20

0

5

10

Order m

(a)

15

20 Order m

25

30

35

40

(b)

˜ m (κ)| for (a) κ = 10 and (b) Fig. A.1. Graphs showing the three phases in the behavior of |Q κ = 20 as the order m is increased.

Moreover, replacing n by n − 2 in the above two identities and inserting them into (A.1), together with (A.3) and (A.4), and simplifying, we arrive at ˜ m (κ) = 1 Q 2



n n−1



(n − 1)! (2n − 2)!

2

(2κ)2n−1 + ··· , 2n − 1

which in terms of m gives the result claimed. ˜ m (κ) also decays algebraically as κ becomes small. The Lemma A.1 shows that Q m ˜ behavior of the ratio Q (κ) when the order of the Bessel functions and their arguments are both very large is shown in Figure A.1 for κ = 10 and κ = 20. In this case two distinct phases are observed depending on the order m. It is evident from both of the ˜ m (κ) initially oscillates around unity. As the graphs of Figure A.1 that the quotient Q order m passes through κ+1, there is a relatively short-lived transition zone where the ratio first increases and then decays at an algebraic rate as the order m is increased. ˜ m (κ) decays Therefore, the ratio is still oscillating in the transition zone. Finally, Q at an exponential rate. The remainder of the section describes the behavior of the ˜ m for high wave numbers, i.e., when κ = ωh/2  1. ratio Q ˜ m (κ) be defined as above and m = integer + 1 . Then, as Theorem A.2. Let Q 2 ˜ m (κ) passes through two phases: m is increased, Q ˜ m (κ) oscillates around unity but does not decay 1. For m < κ + 1 + o(κ1/3 ), Q as m is increased. ˜ m (κ) decays at a superexponential rate: 2. For m > κ + 1 + o(κ1/3 ), Q 1 Q (κ) ≈ 2 ˜m



  1 − k 2 /(m − 2)2 2√1−k2 /(m−2)2 e  1 + 1 − k 2 /(m − 2)2 1−

˜ m (κ) ≈ so that Q

1 2



κe 2(m−1)

2(m−1)

m−2  m  2 1 − 1 − k 2 /m2 2√1−k2 /m2 2 e  1 + 1 − k 2 /m2

∀m  κ + 1.

A.1. Preasymptotic regime: m < κ + 1. In the preasymptotic regime we ˜ m (κ) for uniform asymptotic expansions of Bessel functions discuss the behavior of Q with large order and large argument. Moreover, in this case the value of the order m does not exceed the argument κ + 1 of the Bessel functions. Langer’s formulae given

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3933

in section 7.13.4 of [11] provide uniform asymptotic expansions for Bessel functions of large order and large argument. Inserting these formulae into (A.1) with higher order terms dropped and performing some manipulation gives √ √   J (z ) 3 − Y (z ) − β J (z ) 3 − Y (z ) m−2 m−2 m m 1/3 1/3 1/3 1/3 m ˜ (κ) = √ √ ,  (A.5) Q J1/3 (zm−2 ) + Y1/3 (zm−2 ) 3 − β J1/3 (zm ) + Y1/3 (zm ) 3 where

  1/4 1/2  κ2 − (m − 2)2 zm (2m − 3) (A.6) β= (κ2 − m2 ) (2m − 1) zm−2

−1 2 2 with zm = m(wm − tan

wm ), wm = κ /m − 1, and zm−2 = (m − 2)(wm−2 − −1 tan wm−2 ), wm−2 = κ2 /(m − 2)2 − 1. ˜ m (κ) oscillates A.1.1. Oscillatory phase: m < κ + 1 − o(κ1/3 ). The ratio Q with a magnitude of order unity for m small relative to κ + 1 in the preasymptotic regime. Since, for small m relative to κ + 1, the arguments zm and zm−2 of the Bessel functions appearing in (A.5) are large and positive, we can use the asymptotic expansions of Bessel functions with large arguments given in (8.440)1 and (8.440)2 of [13]:      πz − 12  πz − 12 π π 1 1 cos z − νπ − sin z − νπ − Jν (z) ∼ and Yν (z) ∼ . 2 2 4 2 2 4 With the aid of the above expressions, followed by elementary manipulations, we arrive at   ⎡ π π ⎤    sec z cos z − − 1 − β m m−2 4 4 ⎦ ˜ m (κ) ≈ cot zm−2 − π ⎣   (A.7) Q π π 4 cosec zm−2 − 1 − β  sin zm − 4 4 2

2

) 1/4 (2m−3) ] (2m−1) . It is evident from Figure A.2 that the expression with β  = [ (κ (κ−(m−2) 2 −m2 ) ˜ m (κ) and represents qualitatively very accurate behavior (A.7) agrees closely with Q even for modest values of κ in the preasymptotic regime.

A.1.2. Transition zone: κ + 1 − o(κ1/3 ) < m < κ + 1. Here we consider the behavior when m approaches κ + 1 from the left. Since in this region the value of the expression (A.6) used in (A.5) is approximately 1, using the series representations for Bessel functions given in identity (8.440) of [13, eq. (A.1)] the value becomes  2  zm−2 1/3  zm 1/3 ˜ m (κ) ≈ − √1 + 3 Γ 2 (A.8) Q + ··· . 3 2 2 3 π 1/2 1/2 and ωm−2  [ κ−(m−2) , where both wm and wm−2 Furthermore, ωm  [ κ−m m/2 ] (m−2)/2 ] are of order 1. Consequently, 3 3   κ−m 2 κ − (m − 2) 2 1 2 1 2 3 3 , zm−2  (m − 2)ωm−2  . zm  mωm  3 3 (m/2)1/3 3 3 ((m − 2)/2)1/3

Substituting the values of zm and zm−2 into (A.8), we finally arrive at  2  1/6 √ 2 √ 36 1 1 m ˜ Q (κ) ≈ − √ + Γ κ−m κ−m+2 + ··· , 3 (m(m − 2)) 3 π which is increasing algebraically at a rate of O(m1/3 ) as m increases in this region.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3934

MARK AINSWORTH AND HAFIZ ABDUL WAJID Case: κ = 20

2

Case: κ = 40

5

10

10

0

10

0

10

−2

10

−5

10

−4

10

−10

10 −6

10

−15

10 −8

10

−20

10

−10

10

−25

10

−12

10

−16

10

−30

˜ m (κ)| |Q |Osc| |E|

−14

10

0

5

˜ m (κ)| |Q |Osc| |E|

10

−35

10

15

20 Order m

25

30

35

40

10

0

10

20

30

40 Order m

50

60

70

80

˜ m (κ)| (given in (A.1)) for m = 1, . . . , 2κ, | Osc | (the absolute value of Fig. A.2. Graphs of |Q the right-hand side of (A.7)) for m = 1, . . . , κ + 1, and |E| (the absolute value of the right-hand side of (A.11)) for m = κ + 1, . . . , 2κ. Values of κ = 20 and κ = 40 are shown. Observe the oscillatory ˜ m (κ)| and the good qualitative agreement provided by the | Osc | in the preasymptotic behavior of |Q ˜ m (κ)| and regime m < κ + 1 for ωh ≥ 20. Furthermore, note the quantitative agreement between |Q |E| in the asymptotic regime m > κ + 1.

A.2. Asymptotic regime: m > κ + 1. In this regime the order m of the Bessel functions exceeds the argument κ + 1. Again, using Langer’s formulae with higher order terms dropped from section 7.13.4 of [11] together with (A.1) gives   π −1 K1/3 (zm−2 ) − γK1/3 (zm ) m ˜  , (A.9) Q (κ) = − I1/3 (zm−2 ) + I−1/3 (zm−2 ) − γ I1/3 (zm ) + I−1/3 (zm ) where (A.10)

 γ=

(m − 2)2 − κ2 (m2 − κ2 )

1/4

(2m − 3) (2m − 1)



zm

1/2

zm−2

with zm = m(tanh−1 wm − wm ), wm = 1 − κ2 /m2 , and zm−2 = (m − 2)(tanh−1 wm−2 − wm−2 ), wm−2 = 1 − κ2 /(m − 2)2 . Now, writing z = 23 ξ 3/2 and combining this with (11.1.04) and (11.1.12), given in [26], gives  − 12  − 12 ξm ξm −1 and I1/3 (zm ) + I−1/3 (zm ) = Bi(ξm ) . π K1/3 (zm ) = Ai(ξm ) 3 3 Again, using formulae (11.1.07) and (11.1.16) from [26], we get   3 e−zm −3/4 3 zm −3/4 −1 ξ e ξm . π K1/3 (zm ) ≈ and I1/3 (zm ) + I−1/3 (zm ) ≈ − π 2 m π Replacing m by m − 2 in the above expressions will give us the rest of the values. Now, inserting these values into (A.9) and simplifying gives   1 − γ  e(zm−2 −zm ) e−2zm−2 m ˜ (A.11) Q (κ) ≈ − 2 1 − γ  e−(zm−2 −zm ) with (A.12)



(m − 2)2 − κ2 γ = m2 − κ 2 

1/4

(2m − 3) . (2m − 1)

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3935

A.2.1. Transition zone: κ + 1 < m < κ + 1 + o(κ1/3 ). In this region, we start with the expression 1/2

1/2

˜ m (κ) ≈ − ξm−2 Ai(ξm ) − ξm Ai(ξm−2 ) , Q 1/2 1/2 ξm−2 Bi(ξm ) − ξm Bi(ξm−2 ) which is obtained from (A.9) by using the formulae (11.1.04) and (11.1.12) from [26]. Moreover, using series representations for Airy functions, given in formulae (11.1.07) and (11.1.16) of [26], the above expression gives  2 2 31/3 1 m 1/2 1/2 ˜ Γ ξm ξm−2 , (A.13) Q (κ) ≈ − √ − π 3 3 where the ξm and the zm are related by the expression ξm = (3/2zm )2/3 and ξm−2 = (3/2zm−2)2/3 . Hence, we obtain  ξm 

2 m



1/3 (m − κ)

and ξm−2 

2 m−2

1/3 (m − 2 − κ).

Inserting these values into (A.13), we finally obtain ˜ m (κ) ≈ − √1 − 1 Γ Q 3 π

 2  1/6 √ 2 √ 36 κ−m κ−m+2 + ··· , 3 m(m − 2)

˜ m is decreasing algebraically at a rate of O(m−1/3 ) with increasing m. Therewhere Q ˜ fore, Qm is still oscillating in the transition region. Concluding with the spectral ˜ m oscillates longer compared with the ratio Qm obtained element method, the ratio Q with the finite element method given in [1, eq. (A.9)]. A.2.2. Exponential decay phase: m > κ + 1 + o(κ1/3 ). When m exceeds κ + 1 + o(κ1/3 ), it is easy to verify that γ  , which is defined in (A.12), is approximately equal to 1, and consequently, (A.11) takes the form −(zm−2 +zm ) ˜ m (κ) ≈ e . Q 2

(A.14)

Replacing zm−2 and zm by wm−2 and wm in the expressions for e−zm−2 and e−zm , respectively, we obtain

e

−zm−2



1 − wm−2 2wm−2 = e 1 + wm−2

m − 2 2

and e

−zm

m 1 − wm 2wm 2 = e 1 + wm 

after elementary manipulations. Define a monotonic decreasing function f : w → (1 − w)/(1 + w)e2w on [0, 1] with values ranging from 1 to 0. In the limiting case we find that 2     κe 2   κe f 1 − k 2 /(m − 2)2  and f 1 − k 2 /m2  . 2(m − 2) 2m So, substituting the values of e−zm−2 and e−zm into (A.14) and simplifying gives ˜ m (κ) ≈ 1 [ κe ]2(m−1) . Hence, we obtain a superexponential rate of decay when Q 2 2(m−1) m − 1 > κe/2.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

3936

MARK AINSWORTH AND HAFIZ ABDUL WAJID REFERENCES

[1] M. Ainsworth, Discrete dispersion relation for hp-version finite element approximation at high wave number, SIAM J. Numer. Anal., 42 (2004), pp. 553–575. [2] M. Ainsworth, Dispersive and dissipative behaviour of high order discontinuous Galerkin finite element methods, J. Comput. Phys., 198 (2004), pp. 106–130. [3] M. Ainsworth, Dispersive properties of high-order N´ ed´ elec/edge element approximation of the time-harmonic Maxwell equations, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 362 (2004), pp. 471–491. [4] M. Ainsworth, P. Monk, and W. Muniz, Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation, J. Sci. Comput., 27 (2006), pp. 5–40. [5] I. M. Babuˇ ska and S. A. Sauter, Is the pollution effect of the FEM avoidable for the Helmholtz equation considering high wave numbers?, SIAM J. Numer. Anal., 34 (1997), pp. 2392– 2423. [6] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods, Evolution to Complex Geometries and Applications to Fluid Dynamics, Scientific Computation, Springer, Berlin, 2007. [7] G. Cohen and M. Durufl´ e, Non spurious spectral-like element methods for Maxwell’s equations, J. Comput. Math., 25 (2007), pp. 282–304. [8] G. Cohen and P. Monk, Gauss point mass lumping schemes for Maxwell’s equations, Numer. Methods Partial Differential Equations, 14 (1998), pp. 63–88. [9] G. Cohen and P. Monk, Mur-N´ ed´ elec finite element schemes for Maxwell’s equations, Comput. Methods Appl. Mech. Engrg., 169 (1999), pp. 197–217. [10] G. C. Cohen, Higher-Order Numerical Methods for Transient Wave Equations, Scientific Computation, Springer-Verlag, Berlin, 2002. [11] A. Erd´ elyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi, Higher Transcendental Functions, Vol. I, McGraw-Hill, New York, Toronto, London, 1953. [12] A. Erd´ elyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi, Higher Transcendental Functions, Vol. II, McGraw-Hill, New York, Toronto, London, 1953. [13] I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products, 4th ed., Academic Press, New York, 1965. [14] F. Ihlenburg, Finite Element Analysis of Acoustic Scattering, Appl. Math. Sci. 132, SpringerVerlag, New York, 1998. [15] F. Ihlenburg and I. Babuˇ ska, Finite element solution of the Helmholtz equation with high wave number. I. The h-version of the FEM, Comput. Math. Appl., 30 (1995), pp. 9–37. [16] F. Ihlenburg and I. Babuˇ ska, Finite element solution of the Helmholtz equation with high wave number. II. The h-p version of the FEM, SIAM J. Numer. Anal., 34 (1997), pp. 315– 358. [17] F. Ihlenburg, I. Babuˇ ska, and S. Sauter, Reliability of finite element methods for the numerical computation of waves, Adv. Eng. Soft., 28 (1997), pp. 417–424. [18] M. Iskandarani, D. B. Haidvogel, and J. C. Levin, A three-dimensional spectral element model for the solution of the hydrostatic primitive equations, J. Comput. Phys., 186 (2003), pp. 397–425. [19] M. S. Jensen, High convergence order finite elements with lumped mass matrix, Internat. J. Numer. Methods Engrg., 39 (1996), pp. 1879–1888. [20] D. Komatitsch and J. Tromp, Spectral-element simulations of global seismic wave propagation-I. Validation, Geophys. J. Int., 149 (2002), pp. 390–412. [21] D. Komatitsch and J. Tromp, Spectral-element simulations of global seismic wave propagation-II. 3-D models, oceans, rotation, and self-gravitation, Geophys. J. Int., 150 (2002), pp. 303–318. [22] D. Komatitsch and J. P. Vilotte, The spectral-element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures, Bull. Seismol. Soc. Am., 88 (1998), pp. 368–392. [23] F. Maggio and A. Quarteroni, Acoustic wave simulation by spectral methods, East-West J. Numer. Math., 2 (1994), pp. 129–150. [24] J. M. Melenk and C. Schwab, HP FEM for reaction-diffusion equations. I. Robust exponential convergence, SIAM J. Numer. Anal., 35 (1998), pp. 1520–1557. [25] P. B. Monk and A. K. Parrott, A dispersion analysis of finite element methods for Maxwell’s equations, SIAM J. Sci. Comput., 15 (1994), pp. 916–937. [26] F. W. J. Olver, Asymptotics and Special Functions, Computer Science and Applied Mathematics, Academic Press, New York, London, 1974.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

DISPERSIVE AND DISSIPATIVE BEHAVIOR OF SEM

3937

[27] A. Ralston, A First Course in Numerical Analysis, McGraw-Hill, New York, 1965. [28] D. Stanescu, D. A. Kopriva, and M. Y. Hussaini, Dispersion analysis for discontinuous spectral element methods, J. Sci. Comput., 15 (2000), pp. 149–171. [29] L. L. Thompson and P. M. Pinsky, Complex wavenumber Fourier analysis of the p-version finite element method, Comput. Mech., 13 (1994), pp. 255–275. [30] J. Valenciano and M. A. J. Chaplain, An explicit subparametric spectral element method of lines applied to a tumour angiogenesis system of partial differential equations, Math. Models Methods Appl. Sci., 14 (2004), pp. 165–187. [31] J. L. M. van Dorsselaer and W. Hundsdorfer, Stability estimates based on numerical ranges with an application to a spectral method, BIT, 34 (1994), pp. 228–238. [32] E. Zampieri and L. F. Pavarino, Approximation of acoustic waves by explicit Newmark’s schemes and spectral element methods, J. Comput. Appl. Math., 185 (2006), pp. 308–325. [33] E. Zampieri and L. F. Pavarino, An explicit second order spectral element method for acoustic waves, Adv. Comput. Math., 25 (2006), pp. 381–401.

Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.