Spectral Methods

Claudio Canuto1 and Alfio Quarteroni2,3

1 Polytechnic University of Turin, Turin, Italy
2 Swiss Federal Institute of Technology in Lausanne (EPFL), Lausanne, Switzerland
3 Polytechnic University of Milan, Milan, Italy

approximations based on algebraic polynomial expansions.

1 Introduction 1 The different concepts are first explained on one-dimensional
2 Fourier Methods 1 intervals. Then we address the case of a square or a cube or
3 Algebraic Polynomial Expansion 3 a simplex, and finally the case of more complex geometrical
domains. We illustrate the case of elliptic equations, Stokes
4 Algebraic Expansions on Triangles 5
and Navier–Stokes equations, and then advection equations
5 Stokes and Navier–Stokes Equations 6 and conservation laws.
6 Advection Equations and Conservation Laws 8
7 The Spectral Element Method 11
2 FOURIER METHODS
8 The Mortar Method 12
9 The Spectral Discontinuous Galerkin Method 13 In their early stage, spectral methods were designed to
References 15 approximate the periodic solution of partial differential
equations by a truncated Fourier series. If

∑
+∞
1
u(x) = uk 𝜑k (x), 𝜑k (x) = √ eikx
1 INTRODUCTION k=−∞ 2𝜋

Spectral methods are numerical methods for solving is the unknown solution, the numerical solution is sought in
boundary-value problems for partial differential equations. the form
∑
(N∕2)−1
They have evolved over the years from their noble ancestor, uN (x) = uN,k 𝜑k (x)
the Fourier method, based on trigonometric expansions, k=−N∕2
through the more flexible Galerkin method with Gaussian
integration, all the way maintaining their most distinguished where N is an (even) integer that dictates the size of the
feature: the very high rate of convergence. approximate problem. Note that the unknowns are repre-
For the reader’s convenience, we will gradually approach sented by the Fourier coefficients {uN,k }. Potentially, this
this subject by first addressing the case of periodic prob- approximation has a tremendously high quality, since for all
lems, where the so-called Fourier methods are used. Then 0 ≤ k ≤ s, there exists a positive constant Ck,s such that
we turn to nonperiodic problems and address collocation
inf ||u − vN ||H k (0,2𝜋) ≤ Ck,s N k−s ||u||Hps (0,2𝜋) (1)
vN ∈SN
Encyclopedia of Computational Mechanics Second Edition,
Edited by Erwin Stein, René de Borst and Thomas J.R. Hughes. provided u belongs to the Sobolev space Hps (0, 2𝜋) of peri-
© 2017 John Wiley & Sons, Ltd. odic functions having s derivatives in L2 (0, 2𝜋). Here, SN =
DOI: 10.1002/9781119176817.ecm2003m span{𝜑k | − N∕2 ≤ k ≤ (N∕2) − 1} denotes the space of the
2 Spectral Methods

finite Fourier series of order N. This abstract approxima- (v, w)N = (v, w), ∀v, w ∈ SN , enables us to recover from (6)
tion property is reflected by a corresponding error estimate a collocation formulation LN ucN = f at all nodes xj , where
for the difference u − uN . Actually, in the most classical LN is obtained from L by replacing each derivative by the
approach, the spectral Fourier method consists of approxi- corresponding so-called pseudospectral derivative. This
mating a given PDE (partial differential equation), say means that for any smooth function v, Dv is replaced by
DN v = D(IN v), where
Lu(x) = f (x), x∈Ω (2)
∑
(N∕2)−1

with Ω = (0, 2𝜋), where L is a differential operator and f is IN v ∈ SN , IN v(x) = v∗k 𝜑k (x), v∗k = (v, 𝜑k )N
a given 2𝜋-periodic function, by requiring the L2 -projection k=−N∕2

of the residual upon the subspace SN to vanish, that is, (7)

find uN ∈ SN s.t. (LuN − f , 𝜑) = 0 ∀𝜑 ∈ SN (3) is the interpolant of v at the nodes {xj }.

The interpolation error satisfies, for 0 ≤ k ≤ s, s ≥ 1,
Here, (v, w) = ∫Ω vw̄ denotes the L2 (Ω) inner product.
If L is a constant coefficient operator, this yields an ||v − IN v||H k (0,2𝜋) ≤ Cs N k−s ||v||Hps (0,2𝜋) (8)
embarrassingly simple problem. As a matter of fact,
owing to the L2 -orthogonality of the functions {𝜑k }, that and so does the collocation error u − ucN . A consequence
is, (𝜑k , 𝜑m ) = 𝛿km , ∀k, m ∈ ℤ, equation (3) yields, after of (8) (when k = 1) is that the error on the pseudospectral
Fourier-transforming the residual LuN − f , the following set derivative ||Dv − DN v||L2 (0,2𝜋) decreases like a constant time
of explicit equations for the unknowns: N 1−s , provided that v ∈ Hps (0, 2𝜋) for some s ≥ 1. Actually,
one can even prove that
N N
𝜆k uN,k = f̂k , − ≤k ≤ −1 (4)
2 2 ||Dv − DN v||L2 (0,2𝜋) ≤ C(𝜂)Ne−N𝜂∕2 ∀ 0 < 𝜂 < 𝜂0
where f̂k = (f , 𝜑k ) is the kth Fourier coefficient of f , while 𝜆k provided that v is analytic in the strip |Im z| < 𝜂0 . This
is the kth eigenvalue of L. For instance, if exponential rate of convergence is often referred to as spec-
tral convergence, as it is a distinguishing feature of spectral
L = −D(𝛼D) + 𝛽D + 𝛾I (with 𝛼, 𝛽, 𝛾 constants) (5)
methods.
There is, however, a major difference between the collo-
where D denotes differentiation with respect to x and I the
cation approach and the L2 -projection approach (3). In the
identity, then 𝜆k = 𝛼k2 + i𝛽k + 𝛾.
latter, the unknowns are the frequency coefficients {uN,k } of
Moreover, in this special case, uN indeed coincides with
uN , whereas in the collocation approach one looks for the
the truncated Fourier series of order N of the exact solution u,
nodal values {uj = ucN (xj )} of ucN . These values may be inter-
thus the bound (1) (with vN = uN ) provides an error estimate.
preted as the coefficients of ucN with respect to the trigono-
However, the one that we have just described is an overly
metric Lagrange basis associated with the nodes xj ; indeed,
fortunate circumstance. Should indeed some of the coef-
observing that ucN = IN ucN , using (7) and exchanging summa-
ficients 𝛼, 𝛽, or 𝛾 be functions of x (or, even worse, of
tions over k and j, one gets
u, yielding a nonlinear equation), then convolution sums
between the unknown frequency coefficients {uN,k } and
∑
N−1
2𝜋 ∑
(N∕2)−1
∑
N−1
the Fourier coefficients of 𝛼, 𝛽, 𝛾 will arise, and the diag- ucN (x) = uj 𝜑k (xj )𝜑k (x) = uj 𝜓j (x)
onal structure of equation (4) would be lost. A variant of j=0
N k=−N∕2 j=0
the projection approach (3) can be based on evaluating
the convolution sums by discrete Fourier transform. This where 𝜓j ∈ SN satisfies 𝜓j (xm ) = 𝛿j,m , 0 ≤ j, m ≤ N − 1. A
requires introducing equally spaced nodes, xj = 𝜋j∕N, j = modal representation is used in the former case (Fourier),
0, … , N − 1, then replacing the exact integrals in (3) by whereas a nodal one is adopted in the latter (collocation).
numerical integration and the resulting scheme is The same approach can be pursued for boundary-value
problems set on multidimensional intervals Ω = (0, 2𝜋)d ,
find ucN ∈ SN s.t. (LucN − f , 𝜑)N = 0 ∀𝜑 ∈ SN (6) d = 2, 3 by tensorizing basis functions and collocation nodes.
Fourier methods represent the most classical approach in
∑
where (v, w)N = (2𝜋∕N) N−1 ̄ j ) is the Gaussian
j=0 v(xj )w(x spectral methods. The interested reader can find a compre-
approximation of the scalar product (v, w). The exactness of hensive coverage of the subject in the monographs (Gottlieb
the Gaussian approximation on SN , namely, the property that and Orszag, 1977; Boyd, 2001; Canuto et al., 2006).
 Spectral Methods 3

3 ALGEBRAIC POLYNOMIAL {uN,k } with a structured matrix for which efficient diagonal-
EXPANSION ization algorithms can be devised, a circumstance that is also
featured by the multidimensional problems that are generated
When a boundary-value problem with nonperiodic data (of by tensorization.
Dirichlet, Neumann, or mixed type) has to be solved numer- However, this is not general enough, as this structure gets
ically, the trigonometric expansion is no longer adequate to lost for a more general kind of differential operators. A more
guarantee high order of accuracy. Then, Jacobi orthogonal flexible approach (in analogy with what was done in the
polynomials are used to provide orthogonal bases for the Fourier case) consists of adopting a nodal representation
approximation space. of uN at selected Gauss–Lobatto nodes xj = cos(𝜋j∕N), j =
The finite dimensional approximation space is now the 0, … , N, then looking for a standard Galerkin approxima-
space ℙN made of algebraic polynomials of degree less than tion with integrals replaced by Gauss–Lobatto integration:
or equal to N.
The historical approach, inspired by the Fourier method, ∑
N
(u, v)N = 𝛼j u(xj )v(xj ) (11)
aimed at expanding the approximate solution with respect to
j=0
a basis of orthogonal polynomials
where 𝛼j = 𝜋∕N for j = 1, … , N − 1, 𝛼0 = 𝛼N = 𝜋∕(2N) are
∑
N
the quadrature coefficients.
uN (x) = uN,k pk (x) (9)
k=0
Should we still consider the baby Dirichlet boundary-value
problem for the operator L introduced in (5), the corre-
where uN,k now represent the unknown frequency coeffi- sponding discrete problem would read:
cients.
The matter of choice were the Chebyshev polynomials, find uN ∈ ℙN , uN (−1) = u− , uN (1) = u+ ,
pk (x) = Tk (x) = cos(k𝜃), 𝜃 = cos−1 (x), −1 ≤ x ≤ 1, owing
s.t. (𝛼u′N , v′N )N + (𝛽u′N , vN )N
to their analogy with trigonometric polynomials. Since the
Chebyshev basis does not necessarily match the boundary + (𝛾uN , vN )N = (f , vN )N ∀vN ∈ ℙ0N (12)
requirement (as Tk (1) = 1, Tk (−1) = (−1)k , ∀k ≥ 0), one
trick consists of projecting the equation residual on the where now ℙ0N = {vN ∈ ℙN | vN (−1) = vN (1) = 0}. This
reduced space ℙN−2 , enforcing the boundary conditions time, however, the expansion is made in terms of the nodal
afterward. For instance, for a Dirichlet boundary-value Lagrangian basis at Gauss–Lobatto nodes, that is, using
problem like (2), where now Ω = (−1, 1), and Dirichlet
boundary conditions u(−1) = u− , u(+1) = u+ , the solution ∑
N

(9) is required to satisfy the following equations: uN (x) = uj 𝜓j (x)

j=0

(LuN − f , Tk )𝜔 = 0, 0≤k ≤N−2

instead of (9), where 𝜓j is the unique algebraic polynomial
uN (−1) = u− , uN (+1) = u+ (10) of degree N such that 𝜓j (xi ) = 𝛿ij , ∀i, j = 0, … , N.
One may show that
This modal approach was termed the Lanczos–Tau method.
1 2 ′
The symbol (u, v)𝜔 = ∫−1 uv 𝜔 dx is the so-called weighted (−1)j+1 (1 − x )TN (x)
𝜓j (x) = ⋅ , j = 0, … , N (13)
scalar product with respect to the Chebyshev weight func- cj N 2 (x − xj )
tion 𝜔(x) = (1 − x2 )−1∕2 , −1 < x < 1. The weighted scalar
product is used, instead of the more traditional one (⋅, ⋅), in with cj = 2 if j = 0 or N, cj = 1 otherwise.
order to take advantage (to the highest possible extent) of the The same approximation framework can be set up by
Chebyshev orthogonality, replacing the Chebyshev polynomials with the Legendre
polynomials {Lk , k = 0, 1, … }, which are orthogonal with
(Tk , Tm )𝜔 = 0 if k≠m respect to the traditional L2 -scalar product (otherwise said
(T0 , T0 )𝜔 = 𝜋 with respect to the weight function 𝜔 = 1).
The approximate problem still reads like (12); however,
𝜋
(Tk , Tk )𝜔 = ∀k ≥ 1 this time the nodes {xj } and the coefficients {𝛼j } are those
2
of the (Legendre) Gauss–Lobatto integration.
When L has constant coefficients, the Lanczos–Tau problem The approach described above is named G-NI (Galerkin
(10) yields an algebraic system for the frequency coefficients with numerical integration). A similar G-NI approach can
4 Spectral Methods

be pursued in several dimensions. For instance, consider a on the space V × V, that is,
second-order elliptic boundary-value problem
{ ∃𝛼 ∗ > 0 independent of N s.t.
Lu = f in Ω = (−1, 1)d , d = 2, 3
(14) aN (vN , vN ) ≥ 𝛼 ∗ ||vN ||2H 1 (Ω) , ∀vN ∈ VN (19)
u = 0 on 𝜕Ω
This is the case for the problem at hand if, for example, 𝜷 is
together with its weak form
constant and 𝛾 is non-negative.
The convergence analysis of the G-NI approximation can
find u ∈ V = H01 (Ω) ∶ a(u, v) = (f , v) ∀v ∈ V (15)
be carried out by invoking the Strang lemma for generalized
where the bilinear form a ∶ V × V → ℝ is associated with the Galerkin approximation. Precisely, the following error esti-
operator L; for instance, if mate holds:
[
( )
Lu = −∇ ⋅ (𝛼∇u) + 𝜷 ⋅ ∇u + 𝛾u, with 𝛼 ≥ 𝛼0 > 0 M
||u − uN || ≤ inf 1 + ∗ ||u − wN ||
(16) wN ∈VN 𝛼
]
then a(u, v) = ∫Ω (𝛼∇u ⋅ ∇v + 𝜷 ⋅ ∇uv + 𝛾uv). 1 a(wN , vN ) − aN (wN , vN )
+ ∗ sup
Upon introducing the tensorized Legendre–Gauss–Lobatto 𝛼 vN ∈VN ∖{0} ||vN ||
quadrature nodes and coefficients xk = (xk1 , … , xkd ) and (f , vN ) − (f , vN )N
1
𝛼k = 𝛼k1 , … , 𝛼kd , (ki = 0, … , N), the Legendre–Galerkin + sup
𝛼∗ vN ∈VN ∖{0} ||vN ||
approximation of (15) with numerical integration (G-NI)
becomes
where || ⋅ || is the norm of H 1 (Ω) and M is the continuity
find uN ∈ VN = ℙ0N ∶ constant of the bilinear form a(⋅, ⋅).
Three sources contribute to the approximation error:
aN (uN , vN ) = (f , vN )N ∀vN ∈ VN (17)
• the best approximation error, which can be immediately
where ℙN is now the set of polynomials of degree ≤ N bounded by taking wN = IN−1 u:
with respect to each of the independent variables, and ℙ0N
is its subspace made of those polynomials vanishing at 𝜕Ω.
inf ||u − wN || ≤ ||u − IN−1 u||
Moreover, ∑ wN ∈VN
(u, v)N = 𝛼k u(xk )v(xk ) (18)
k • the error on the numerical quadrature, which can be
bounded as follows:
is the Gauss–Lobatto quadrature formula that approximates
the scalar product (u, v), while aN is the discrete bilinear (f , vN ) − (f , vN )N
form that is obtained from a by replacing each scalar product sup
vN ∈VN ∖{0} ||vN ||
(⋅, ⋅) with (⋅, ⋅)N . Owing to the property that the quadrature ( )
formula (18) has degree of exactness 2N − 1, the Galerkin ≤ C2 ||f − IN f ||L2 (Ω) + ||f − PN−1 f ||L2 (Ω)
numerically integrated problem (17) can still be interpreted
as a collocation method. Indeed, it follows from (17) that where PN−1 f is the truncated Legendre series of f of order
LN uN = f at all internal nodes xk , where LN is the approx- N − 1;
imation of L obtained by replacing each exact derivative by • the error generated by the approximation of the bilinear
the pseudo-spectral derivative, that is, the derivative of the form, on its hand, is less immediate to estimate. However,
interpolant IN at the Gauss–Lobatto nodes. The interpola- having chosen wN = IN−1 u, which is a polynomial of
tion operator IN is defined as follows: IN v ∈ ℙN for all v ∈ degree N − 1, using the degree of exactness of the
̄ such that IN v(xk ) = v(xk ) for all k. Then, the operator
C0 (Ω), quadrature formula and assuming that the coefficients
approximating (16) is of the operator are constant, one easily checks that
a(wN , vN ) − aN (wN , vN ) = 0, that is, this error is actu-
LN uN = −∇ ⋅ [IN (𝛼∇uN )] + 𝜷 ⋅ ∇uN + 𝛾uN ally null. If the coefficients are nonconstant, one can
control it in terms of the interpolation error measured in
Existence and uniqueness of the solution of (18) follow H 1 (Ω). More details can be found in Canuto et al. (2006),
from the assumption that aN (⋅, ⋅) is a uniformly coercive form Example 5 in Section 6.4.3.
 Spectral Methods 5

We can conclude by taking advantage of the optimality of mapping. On the other hand, triangles, tetrahedra, prisms,
the truncation error in the L2 -norm and that of the interpola- and similar figures allow one to handle complex geometries
tion error in both the L2 - and H 1 -norm: in a more flexible way. So, a natural question arises: can one
match the advantages of a tensor product structure with those
∀f ∈ H r (Ω), r ≥ 0, ||f − PN f ||L2 (Ω) ≤ C3 N −r ||f ||H r (Ω) of a triangular geometry?
A positive answer to this question was given by Dubiner
∀g ∈ H s (Ω), s > 1, N||g − IN g||L2 (Ω) + ||g − IN g||H 1 (Ω)
(1991), who advocated the use of collapsed Cartesian coor-
≤ C4 N 1−s ||g||H s (Ω) dinate systems and warped tensor products for spectral
discretizations. The method was further developed by Karni-
Thus, we obtain that adakis and Sherwin (1999). We now describe this approach
in 2D, pointing to the latter reference for the 3D extensions.
||u − uN || ≤ C5 (N −r ||f ||H r (Ω) + N 1−s ||u||H s (Ω) ) Let us introduce the reference triangle  = {(x1 , x2 ) ∈ ℝ2 ∶
−1 < x1 , x2 ; x1 + x2 < 0}, as well as the reference square
provided u and f have the requested regularity.  = {(𝜉1 , 𝜉2 ) ∈ ℝ2 ∶ −1 < 𝜉1 , 𝜉2 < 1}. The mapping
A few comments on the implementation of the method are
in order. The algebraic system associated with (17) reads 1 + x1
(x1 , x2 ) → (𝜉1 , 𝜉2 ), 𝜉1 = 2 −1
Au = f, where aij = aN (𝜓j , 𝜓i ), fi = (f , 𝜓i )N , u = (ui ), ui = 1 − x2
uN (xi ), and {𝜓j } denote the Lagrangian basis functions of
𝜉2 = x2 (20)
SNd associated with all the nodal points {xj } not sitting on
𝜕Ω. The matrix A, which is nonsingular whenever (19) is
is a bijection between  and . Its inverse is given by
fulfilled, is ill-conditioned: indeed, there exist two constants
C1 , C2 such that 1
(𝜉1 , 𝜉2 ) → (x1 , x2 ), x1 = (1 + 𝜉1 )(1 − 𝜉2 )
2
C1 N 3 ≤ cond(A) ≤ C2 N 3
x2 = 𝜉2
where cond(A) is the (spectral) condition number of A. The
Note that the mapping (x1 , x2 ) → (𝜉1 , 𝜉2 ) sends the ray in 
use of a preconditioned iterative procedure (e.g., the conju-
issuing from the upper vertex (−1, 1) and passing through
gate gradient when 𝜷 = 𝟎, or a Krylov iteration otherwise) is
the point (x1 , −1) into the vertical segment in  of equation
mandatory. A possible preconditioner is given by the diag-
𝜉1 = x1 . Consequently, the transformation becomes singular
onal of A. This yields a preconditioned system whose condi-
at the upper vertex, although it stays bounded therein.
tion number behaves like a constant times N 2 . A more drastic
The Jacobian of the inverse transformation is given by
improvement would be achieved by taking as a precon-
[𝜕(x1 , x2 )∕𝜕(𝜉1 , 𝜉2 )] = (1∕2)(1 − 𝜉2 ). We term (𝜉1 , 𝜉2 ) the
ditioner the matrix associated with the (piecewise-linear)
collapsed Cartesian coordinates of the point on the triangle
finite element discretization of the operator (16) at the
whose regular Cartesian coordinates are (x1 , x2 ).
same Legendre–Gauss–Lobatto nodes. This is an optimal
Denote by {P(𝛼,𝛽) k
(𝜉)} the family of Jacobi polynomials of
preconditioner as the condition number of the preconditioned
increasing degree k ≥ 0, which forms an orthogonal system
system becomes independent of N.
with respect to the measure (1 − 𝜉)𝛼 (1 + 𝜉)𝛽 d𝜉 in (−1, 1)
Spectral methods based on algebraic polynomials have
(note that Pk(0,0) (𝜉) is the Legendre polynomial Lk (𝜉) intro-
been discussed and analyzed in Canuto et al. (1988, 2006),
duced in the previous section). For k = (k1 , k2 ), we set
Bernardi and Maday (1992, 1997), Guo (1998), and Shen
et al. (2011) (see Interpolation and Quasi-Interpolation
Ψk1 (𝜉1 ) ∶= Pk(0,0) (𝜉1 ) (21)
in h- and hp-Version Finite Element Spaces). 1

(2k1 +1,0)
Ψk1 ,k2 (𝜉2 ) ∶= (1 − 𝜉2 )k1 Pk (𝜉2 )
2

4 ALGEBRAIC EXPANSIONS ON and we define the warped tensor product function on 

TRIANGLES
Φk (𝜉1 , 𝜉2 ) ∶= Ψk1 (𝜉1 )Ψk1 ,k2 (𝜉2 ) (22)
Spectral methods for multidimensional problems rely their
efficiency on the tensor product structure of the expan- which is a polynomial of degree k1 in 𝜉1 and k1 + k2 in
sions they use. This feature naturally suggests the setting of 𝜉2 . By applying the mapping (20), one obtains the function
the methods on patches of Cartesian products of intervals, 𝜑k (x1 , x2 ) ∶= Φk (𝜉1 , 𝜉2 ) defined on  . It is easily seen that
such as squares or cubes, possibly after applying a smooth 𝜑k is a polynomial of global degree k1 + k2 in the variables
6 Spectral Methods

x1 , x2 . Furthermore, owing to the orthogonality of Jacobi associated with the quadrature nodes. This means that the
polynomials, one has for k ≠ h G-NI method based on the quadrature formula described
above cannot be equivalent to a collocation method at the
𝜑k (x1 , x2 )𝜑h (x1 , x2 ) dx1 dx2 quadrature points.
∫ Finally, we observe that the G-NI mass and stiffness
1
1 matrices on  can be efficiently built by exploiting the tensor
= P(0,0) (𝜉1 )Ph(0,0) (𝜉1 ) d𝜉1 product structure of both the basis functions and the quadra-
2 ∫−1 k1 1
ture points through the sum-factorization technique.
1
(2k1 +1,0) (2h1 +1,0)
× Pk (𝜉2 )Ph (𝜉2 )(1 − 𝜉2 )k1 +h1 +1 d𝜉2 = 0
∫−1 2 2

5 STOKES AND NAVIER–STOKES

We conclude that the set {𝜑k ∶ 0 ≤ k1 , k2 ; k1 + k2 ≤ N} is an
orthogonal basis of modal type of the space N ( ) of the
EQUATIONS
polynomials of total degree ≤ N in the variables x1 , x2 .
Spectral methods are very popular among the community
While orthogonality simplifies the structure of mass and
of fluid-dynamicists. Owing to their excellent approximation
stiffness matrices, it makes the enforcement of boundary
properties, spectral methods can, in fact, provide very accu-
conditions, or matching conditions between elements,
rate simulations of complex flow patterns. However, special
uneasy. To overcome this difficulty, it is possible to modify
care is needed for the treatment of the incompressibility
the previous construction by building a new modal basis,
constraint. With the aim of simplification, let us first address
say {𝜑m k
}, made of boundary functions (3 vertex func-
the linear Stokes equations
tions plus 3(N − 1) edge functions) and internal functions
(bubbles). Each basis function retains the same “warped
⎧−𝜈Δu + grad p = f in Ω
tensor product” structure as above. Indeed, it is enough to ⎪
replace in one dimension the Jacobi basis P(𝛼,0) k
(𝜉) (with ⎨div u = 0 in Ω (23)
𝛼 = 0 or 2k + 1) with the modified basis given by the two ⎪u = 𝟎 on 𝜕Ω
⎩
boundary functions (1 + 𝜉)∕2 and (1 − 𝜉)∕2 and the N − 1
bubbles (1 + 𝜉)(1 − 𝜉)P(𝛼,1)
k−1
(𝜉), k = 1, … , N − 1. These where Ω = (−1, 1)d for d = 2, 3, whereas 𝜈 > 0 is the kine-
univariate functions are then combined as in (22) to form the matic fluid viscosity, u the fluid velocity, p the fluid pressure,
two-dimensional basis. and f the volume forces.
With such basis on hand, one can discretize a boundary- A natural spectral G-NI discretization reads as follows: find
value problem by the Galerkin method with numerical inte- uN ∈ VN , pN ∈ QN such that
gration (G-NI). To this end, one needs a high-precision
quadrature formula on  . Since ⎧((𝜈∇uN , ∇vN ))N − (pN , div vN )N
⎪
⎨= ((f, vN ))N ∀vN ∈ VN (24)
f (x1 , x2 ) dx1 dx2 ⎪−(q , div v ) = 0
∫ ⎩ N N N ∀qN ∈ QN
1 1
1 where (⋅, ⋅)N is the discrete Gauss–Lobatto scalar product
= d𝜉1 F(𝜉1 , 𝜉2 )(1 − 𝜉2 ) d𝜉2
2 ∫−1 ∫−1 (18), whereas ((⋅, ⋅))N denotes its generalization to the case
of vector functions. Moreover, VN = (ℙ0N )d while QN is a
it is natural to use a tensor product Gaussian formula in  for
polynomial space that needs to be chosen conveniently so
the measure d𝜉1 (1 − 𝜉2 ) d𝜉2 . This is obtained by tensorizing
as to satisfy the following Brezzi condition:
an (N + 1)-point Gauss–Lobatto formula for the measure
d𝜉1 with an N-point Gauss–Radau formula for the measure
∃𝛽N > 0 ∶ ∀qN ∈ QN , ∃vN ∈ VN s.t.
(1 − 𝜉2 ) d𝜉2 with 𝜉2 = −1 as integration knot (excluding the
singular point 𝜉2 = 1 from the integration knots makes the (qN , div vN )N ≥ 𝛽N ||qN ||L2 (Ω) ||vN ||H 1 (Ω) (25)
construction of the matrices easier). The resulting formula
is exact for all polynomials in  of degree ≤ 2N − 1 in The violation of this condition, that is, the existence of
each variable 𝜉1 , 𝜉2 ; hence, in particular, it is exact for all nonconstant pressures qN ∈ QN such that (qN , div vN )N =
polynomials in  of total degree ≤ 2N − 1 in the variables 0 ∀vN ∈ VN , implies the existence of spurious pressure
x1 , x2 . Note, however, that the number of quadrature nodes in modes, which pollute the computed pressure pN .
 is N(N + 1), whereas the dimension of N ( ) is (1∕2)(N + The largest constant 𝛽N , called the inf-sup constant,
1)(N + 2); thus, no basis in N ( ) can be the Lagrange basis depends on the way QN is chosen, and has a special role
 Spectral Methods 7

in the analysis of the spectral approximation (24). Two The pressure matrix S has (N − 1)2 rows and columns. It is
choices are commonly proposed in practice. The first one symmetric; moreover, it is positive definite iff Ker DT = 0, a
is QN = ℙN−2 ∩ L02 (Ω), that is, the space of polynomials of condition that is equivalent to (25).
degree N − 2 in each variable with zero average in Ω. In that If we consider the generalized eigenvalue problem Sw =
case, 𝛽N ≃ CN (1−d)∕2 . 𝜆Mw, where M is the pressure mass matrix (𝜓j , 𝜓i )N , ({𝜓j }
An alternative approach consists of choosing QN = ℙ[𝜆N] ∩ being the Lagrangian polynomials (of degree ≤ N − 2) asso-
ℙN−2 ∩ L02 (Ω), for some 𝜆 ∶ 0 < 𝜆 < 1, where [𝜆N] denotes ciated with the interior Gauss–Lobatto nodes), then the
the largest integer ≤ 𝜆N; in this case, 𝛽N ≥ 𝛽 > 0. maximum generalized eigenvalue 𝜆max is uniformly bounded
The latter approach allows one to derive uniform stability (from above) by the coercivity constant 𝛼 of the discrete
and optimal error bounds for the approximate solution. In bilinear form (∇uN , ∇vN )N (we can assume 𝛼 = 1 in the case
general, this occurs when 𝛽N is uniformly bounded from at hand), whereas the minimum one 𝜆min is proportional to
below as N increases. In fact, using (25), one can obtain that 𝛽N2 . As a consequence, the condition number of the matrix
M −1 S is ∼ 𝛽N−2 , thus ∼ N d−1 in the case of the ℙN − ℙN−2
𝜈||∇uN ||(L2 (Ω))d×d + 𝛽N ||pN ||L2 (Ω) ≤ C||f||(C0 (Ω))
̄ d discretization.
Since S is close to M (the discrete variational equivalent of
where C is a constant independent of N. the identity operator), M can serve as a preconditioner for a
On the contrary, under the assumption (25), the error esti-
conjugate gradient solution of (27). The corresponding PCG
mate on the velocity field is optimal, whereas the error on
(preconditioned conjugate gradient) method will converge
the pressure undergoes a loss of accuracy of order 𝛽N−1 .
in O(N 1∕2 ) iterations for 2D problems and in O(N) for 3D
For instance, in the case where QN = ℙN−2 ∩ L02 (Ω), the
ones. In practice, however, the convergence is faster, as the
following error bound can be proven, provided the assumed
previous estimate on the asymptotic behavior of 𝛽N is too
regularity for the exact solution u, p, and the forcing term f
pessimistic.
holds for suitable values of s ≥ 1 and t ≥ 0:
Several kinds of generalizations are in order.
First of all, we mention that the Stokes system (23) could be
||u − uN ||(H 1 (Ω))d + N (1−d)∕2 ||p − pN ||L2 (Ω)
reduced to a single (vector) equation by L2 -projection upon
≤ CN 1−s (||u||(H s (Ω))d + ||p||H s−1 (Ω) ) + N −t ||f||(H t (Ω))d the divergence-free subspace Vdiv = {v ∈ (H01 (Ω))d | div v =
0 in Ω}:
Note that the (N + 1)2 Gauss–Lobatto nodes are used to inter-
polate the discrete velocity components, while the subset
find u ∈ Vdiv ∶ 𝜈∇u ⋅ ∇v = f ⋅ v ∀v ∈ Vdiv
made of the (N − 1)2 interior Gauss–Lobatto nodes can be ∫Ω ∫Ω
used to interpolate the discrete pressure. Alternatively, one
could use a staggered grid made of the (N − 1)2 Gauss nodes Since this is a well-posed elliptic problem, a unique velocity
for the pressure [and change in (24) the discrete integrals field can be obtained and, afterward, a unique pressure
(pN , div vN )N and (qN , div uN )N accordingly]. This, however, p can be recovered in L02 (Ω). The simple structure of the
would require interpolation between meshes, as, in this case, reduced problem calls for a Galerkin (or G-NI) discretiza-
velocity and pressure feature nodal representations with tion. However, a computer implementation is far from
respect to different sets of nodes. being trivial, as one should construct a set of polynomial
The algebraic formulation of the discrete Stokes problem basis functions that are inherently divergence-free. This
(24) yields the classical block structure matrix form task has been successfully accomplished only for some
[ ] [ ] [ ] specific boundary-value problems, for instance, when Ω is
A DT u f a cylindrical domain and Fourier expansion in the angular
= (26) direction is combined with an expansion in terms of Cheby-
D 0 p 𝟎
shev polynomials in both the longitudinal and the radial
where we have used test functions vN based on the direction. A similar idea is behind the approach by Batcho
Lagrangian polynomials of degree N of the velocity approx- and Karniadakis (1994) to generate eigenfunctions of a
imation, and test functions qN based on the Lagrangian generalized Stokes operator and use them as polynomial
polynomials of degree N − 2 (at the interior nodes) of the divergence-free functions. Recovering the pressure function
pressure approximation. numerically is also quite involved.
Upon eliminating (although only formally!) the u vector, A different kind of generalization consists of using
one obtains from (26) the reduced pressure system equal-order interpolation ℙN − ℙN for both discrete velocity
and pressure fields. However, this choice would give rise
Sp = g, with S = DA−1 DT and g = DA−1 f (27) to a couple of subspaces, VN and QN , which violate the
8 Spectral Methods

Brezzi condition (25), yielding spurious pressure modes that where M is now the velocity mass matrix, while A, DT , and
swamp the physically relevant pressure. In line with what D are the matrices introduced before. To increase time accu-
is nowadays common practice in the finite element commu- racy, a BDF3 discretization is coupled with an extrapolation
nity, Canuto and van Kemenade (1996) have proposed of the nonlinear term. This gives (Karniadakis et al., 1991)
and analyzed a stabilization by bubble functions. The ( )
idea consists in adding to (ℙ0N )d a supplementary space 11
M + A un+1 + DT pn+1
spanned by local polynomial functions having support 6Δt
( )
in one small element called cell. A cell is a Cartesian 1 3 1
= M 3un − un−1 + un−2 + Mf n+1
product of intervals, each one having two consecutive Δt 2 3
Gauss–Lobatto nodes as endpoints. The new velocity space −(3C(un ) − 3C(un−1 ) + C(un−2 ))
is now given by VN = (ℙ0N )d ⊕ BdN , where BdN denotes the
Dun+1 = 0
space of bubble functions, while the pressure space is simply
QN = ℙN ∩ L02 (Ω). The continuity equation is modified by This scheme is third-order accurate with respect to Δt.
the presence of the additional term which plays the role of a The previous schemes are examples of coupled methods,
stabilizing term to damp the oscillatory pressure modes. in which at each time step one solves simultaneously for
The Navier–Stokes equations velocity and pressure. As an alternative, one can resort
to the more efficient yet possibly less accurate segregated
⎧𝜕t u − 𝜈Δu + C(u) + grad p = f in Ω methods, which split the Navier-Stokes system into a cascade
⎪ of subproblems involving either the velocity or the pressure
⎨div u = 0 in Ω (28)
as the unknown. The progenitor of this class of schemes is
⎪u = 𝟎 on 𝜕Ω
⎩ the celebrated Chorin-Temam method (1968). A thorough
discussion of various alternatives can be found in Guer-
differ from (23) because of the presence of the acceleration mond et al. (2006). A related approach is based on the alge-
term 𝜕t u and the convective term C(u). The latter can take the braic factorization of the discrete system obtained from the
standard convective form u ⋅ ∇u; however, other expressions Navier-Stokes system, see Quarteroni et al. (1999, 2000),
are used as well, such as the conservative form div(uu) or Gervasio et al. (2006), Gervasio and Saleri (2006), Gervasio
the skew-symmetric form (1∕2)(u ⋅ ∇u + div(uu)). The three (2008).
forms are all equivalent for the continuous equations (with An extensive coverage of the spectral method for
homogeneous Dirichlet boundary conditions) because of the Navier–Stokes equations can be found in the books Canuto
et al. (1988, 2006, 2007), Deville et al. (2002), and Peyret
incompressibility condition. However, this is no longer true
(2002). For their analysis, see also Bernardi and Maday
at the discrete level. Indeed, the G-NI spectral discretization
(1992, 1997) (see Spectral Element and hp Methods;
of the Navier–Stokes equations has different stability prop-
Mixed Finite Element Methods).
erties depending upon which form of C(u) is employed.
For the time discretization of (28), fully implicit methods
would produce a nonsymmetric, nonlinear system. To
avoid that, the convective term must be treated explic- 6 ADVECTION EQUATIONS AND
itly. One way is to combine backward-difference (BDF) CONSERVATION LAWS
discretization of linear terms with Adams–Bashforth (AB)
discretization of the convective one. A classical recipe is the In order to illustrate spectral approximations to hyperbolic
so-called BDF2/AB3 scheme, that is, the combination of the problems, we consider the linear and nonlinear 1D model
equations ut + aux = 0 and ut + f (u)x = 0, supplemented by
second-order BDF discretization with the third-order AB
initial and appropriate boundary conditions. In addition to
discretization:
the standard issues related to spectral discretizations (effi-
( ) cient implementation, imposition of boundary conditions,
3
M + A un+1 + DT pn+1 stability, accuracy for smooth solutions), here we face a new
2Δt problem. Indeed, the equation may propagate singularities
( )
1 1
= M 2un − un−1 + Mf n+1 along characteristics, or even (in the nonlinear case) generate
Δt 2 singularities from smooth initial data. So, the question arises:
( )
23 4 5 what is the interest of using high-order methods in such
−2 C(un ) − C(un−1 ) + C(un−2 )
12 3 12 cases? We will answer this question in the second part of the
Dun+1 = 0 present section.
 Spectral Methods 9

For periodic problems, say in (0, 2𝜋), the Fourier-Galerkin term in the skew-symmetric form
method is the conceptually simplest choice: find
uN = uN (t) ∈ SN such that 1 1 1
aux = (au)x + aux − ax u (29)
2 2 2
(uN,t + auN,x , 𝜑) = 0 or (uN,t + f (uN )x , 𝜑) = 0, and applying pseudospectral derivatives, that is, the deriva-
∀𝜑 ∈ SN tive of the interpolant, one recovers the same stability esti-
mates as for the pure Galerkin method (in practice, such an
expensive form is rarely necessary). Again, similar consider-
Taking 𝜑 = uN , integrating by parts, and using periodicity,
ations apply in the nonlinear case as well.
one obtains (d∕dt)||uN (t)||2L2 (0,2𝜋) ≤ K||uN (t)||2L2 (0,2𝜋) (with
We now turn to the discretization of nonperiodic prob-
K = max[0,2𝜋] |ax |) for the linear advection equation, and
lems, in the framework of Legendre methods. The advec-
(d∕dt)||uN (t)||2L2 (0,2𝜋) = [F(uN (t))]2𝜋 = 0 (where F denotes
0 tion equation is well-posed, provided we prescribe the solu-
any primitive of f ) for the conservation law. This proves the
tion, say u(xb ) = gb , at the inflow points xb ∈ B− , where
L2 -stability of the approximation.
B± = {xb ∈ { − 1, 1} ∶ (±1)a(xb )nb > 0} with nb = xb . The
In terms of Fourier coefficients, the Galerkin method for
most obvious way to account for the boundary conditions
the advection equation is equivalent to the set of ordinary
is to enforce them exactly (or strongly) in the discrete solu-
differential equations
tion: uN ∈ ℙN satisfies uN (xb ) = gb , ∀xb ∈ B− . The corre-
sponding Galerkin method is L2 -stable. Indeed, assuming for
N N
(uN,k )′ + (auN,x )k = 0, − ≤k ≤ −1 simplicity gb = 0, we take uN itself as the test function and
2 2
after integration by parts we get
Setting for simplicity b = uN,x , we have (ab)k =
∑(N∕2)−1 1d 1
||uN ||2L2 (−1,1) − K||uN ||2L2 (−1,1)
a b . This is a family of convolution sums,
h=−N∕2 k−h h 2 dt 2
which can be computed in O(N 2 ) operations. A more ∑
+ a(xb )nb u2N (xb ) ≤ 0
efficient scheme consists of transforming back a and b in
xb ∈B+
physical space, taking the pointwise product at the nodes
xj = 𝜋j∕N, j = 0, … , N − 1, and returning to Fourier space. whence stability easily follows. A similar result holds for
Using the FFT (fast Fourier transform), the full process the Galerkin method with numerical integration (G-NI) at
costs O(N log N) operations. This is the pseudospectral the Legendre–Gauss–Lobatto points, provided we use the
evaluation of convolutions sums. There is an error involved, skew-symmetric form (29) of the convective term. The G-NI
since one replaces the exact projection PN (ab) of ab upon scheme is equivalent to enforcing the equation at the internal
SN by its interpolant IN (ab) at the nodes. Such error, termed nodes and at the noninflow boundary points (xb ∉ B− ).
the aliasing error, is negligible if N is so large that the A more flexible way to handle the boundary conditions,
essential features of u are resolved. Otherwise, appropriate useful, for example, in domain decomposition and for
de-aliasing techniques can be applied, such as increasing the systems of equations, is to enforce them in a weak sense.
number of interpolation nodes. The rationale is that if stability holds then accuracy is
This process applies to the conservation law as well, assured, provided the boundary conditions are matched to
provided the nonlinearity is polynomial (as for Burgers’ within the same consistency error as for the equation in the
equation, where f (uN ) = (1∕2)u2N , or for the convective term interior. Thus, we seek uN = uN (t) ∈ ℙN satisfying, for all
uN ∇uN in the Navier–Stokes equations). It can be extended vN ∈ ℙN ,
to the nonperiodic case by using the Chebyshev nodes xj =
∑
cos 𝜋j∕N, j = 0, … , N. (uN,t , vN )N − (uN , avN,x )N + a(xb )nb uN (xb )vN (xb )
The Fourier–Galerkin method with the pseudospectral xb ∈B+
evaluation of convolutions sums is nothing but the Galerkin ∑
method with numerical integration described in (6), or equiv- = |a(xb )nb |gb vN (xb ) (30)
xb ∈B−
alently, the collocation method at the quadrature points
This G-NI formulation follows by integrating by parts the
ucN,t (xj ) + a(xj )ucN,x (xj ) = 0, j = 0, … , N − 1 convective term (for simplicity, we assume a to be constant:
otherwise, we use the skew-symmetric form (29)). Choosing
Unless a(x) ≥ 𝛼 > 0 for all x, this scheme is (weakly) as vN the Lagrange characteristic polynomial at each quadra-
unstable due to the aliasing error. Writing the convective ture node, we see that the advection equation is enforced at
10 Spectral Methods

all internal and noninflow nodes, whereas at the inflow nodes |k| ≤ m, Qm,k = 1 − (m∕|k|)(2s−1)∕𝜗 if |k| > m. Thus, the
we have sth order artificial viscosity is applied only to sufficiently
high-frequency modes. For s= 1, one can prove that the
1
uN,t (xb ) + auN,x (xb ) + |a(xb )nb |(uN (xb ) − gb ) = 0 solution is bounded in L∞ (0, 2𝜋), it satisfies the estimate
wb √
||uN (t)||L2 (0,2𝜋) + 𝜖N ||uN,x (t)||L2 (0,2𝜋) ≤ C||uN (0)||L2 (0,2𝜋) ,
and it converges to the correct entropy solution.
Since 1∕wb ∼ cN 2 as N → +∞, this shows that the boundary
A computationally simpler and widely used road to
condition is indeed enforced by a penalty method. The
stability of the scheme (30) immediately follows by taking stabilization consists of filtering the spectral solution when
vN = uN . Stability is actually guaranteed even if we multiply advancing in time,
each boundary term in (30) by any constant 𝜏b ≥ 1∕2, thus
enhancing the flexibility of the penalty method. ∑
(N∕2)−1 ( )
2k
Spectral methods for linear advection equations are uN (t) → N uN (t) = 𝜎 uN,k (t) exp(ikx)
k=−N∕2
N
addressed by Gottlieb and Hesthaven (2001), Funaro (1997),
and Fornberg (1996).
From the computational point of view, a simple yet efficient where 𝜎 = 𝜎(𝜂) is a smooth, even function satisfying 𝜎(0) =
device to stabilize spectral collocation methods for advec- 1, 𝜎 (j) (0) = 0 for all j with 1 ≤ j ≤ s (for a suitable s), mono-
tion or advection-dominated equations has been proposed tonically decreasing for 𝜂 > 0 and vanishing (or being expo-
by Fischer and Mullen (2001): it consists in applying after nentially small) for 𝜂 > 1. A popular choice is the exponen-
each time step a filtering operator of the form 𝛼IN−1 + tial filter 𝜎(𝜂) = exp(−𝛼𝜂 2s ). Interestingly, the effect of the
(1 − 𝛼)Id, where IN−1 is the interpolation operator at the spectral viscosity correction described above can be closely
Legendre-Gauss-Lobatto grid of order N − 1 (if the similar mimicked by applying the exponential filter with 𝜎(2k∕N) =
grid of order N is used to collocate the equation) and 𝛼 ∈ exp(−𝜖N Qm,k k2 ).
(0, 1). The effect is a dumping in the spurious high-frequency If the solution of the conservation law is piecewise analytic
oscillations. but discontinuous, its truncation PN u or its interpolation IN u
Let us now consider the nonlinear conservation law are highly oscillatory around the singularities, and converge
ut + f (u)x = 0. The stability (and convergence) of spectral slowly (O(N −1 )) to u away from them. However, they contain
discretizations is a much more delicate issue than that for the enough information to allow the reconstruction of the exact
linear advection equation. Indeed, the equation may develop solution with exponential accuracy, away from the singu-
singular solutions at a finite time, which correspond to the larities, by a postprocessing as described below. It follows
accumulation of energy in the high-frequency modes or, that the crucial feature of the discretization scheme is the
equivalently, to the onset of oscillations around discontinu- capability of producing an approximation uN , which is spec-
ities (Gibbs phenomenon). The nonlinear mechanism may trally close to PN u or to IN u. This is precisely what is
amplify the high-frequency components, leading to destruc- obtained by the spectral viscosity method or by the equiv-
tive instabilities (in stronger norms than L2 ). On the other alent filtering procedure. Given PN u (similar considerations
hand, oscillations should not be brutally suppressed: they are apply to IN u), the postprocessing reconstruction may be
inherent to the high-order representation of discontinuous local or global. In the former case, a spectrally accurate
functions, and they may hide the correct information that approximation of u at a point x0 of analyticity is given by
allows the reconstruction of the exact solution. Thus, a good u∗N (x0 ) = ∫ℝ K𝜈 (x0 , y)𝜚(x0 − y)PN u(y) dy, where 𝜈 = [𝛽N] for
spectral discretization should guarantee enough stability some 𝛽 ∈ (0, 1), K𝜈 (x, y) is, for each x, a 𝜈-degree polynomial
while preserving enough accuracy. Furthermore, the discrete approximation of the delta at x (e.g., for Fourier, K𝜈 (x, y) =
∑𝜈∕2
solution should converge to the physically relevant exact 1 + k=0 cos(x − y) is the Dirichlet kernel), whereas 𝜚(𝜂)
solution by fulfilling an appropriate entropy condition. is a C∞ -localizer around 𝜂 = 0. In the latter case, a spec-
The mathematically most rigorous discretization that trally accurate approximation of u on an interval [a, b] of
matches these requirements is the spectral viscosity method analyticity is given (Gottlieb and Shu, 1997) by the orthog-
(Tadmor, 1998). In the Fourier–Galerkin context, it amounts onal projection of PN u upon ℙ𝜈 ([a, b]) (again 𝜈 = [𝛽N]) with
b
to considering the modified equation respect to the weighted inner product ∫a u(x)v(x)𝜔𝜈 (x) dx,
with 𝜔𝜈 (x) = ((x − a)(b − x))𝜈−1∕2 , which varies with N. The
uN,t + (PN f (uN ))x = 𝜖N (−1)s Dsx (Qm Dsx uN ) projection is computed via the Gegenbauer polynomials (i.e.,
(𝜈−1∕2,𝜈−1∕2)
the Jacobi polynomials {Pk }) translated and scaled
where 𝜖N ∼ cN 1−2s , m = mN ∼ N 𝜗 for some 𝜗 < 1 − 1∕(2s), to [a, b]. The reader can refer, for example, to Gottlieb and
and the Fourier coefficients of Qm satisfy Qm,k = 0 if Shu (1997), Gottlieb and Tadmor (1984), and Tadmor (1998).
 Spectral Methods 11

7 THE SPECTRAL ELEMENT METHOD FEMs of p-type are defined in terms of the Legendre poly-
̂ Precisely,
nomials Lk (𝜉) of degree k (k = 2, … , p), 𝜉 ∈ Ω.
The spectral element method (SEM) represents another ̂
the p + 1 modal basis functions on Ω are defined by
example of the Galerkin method. The finite dimensional
space is now made of piecewise algebraic polynomials of 1−𝜉 1+𝜉
𝜑1 (𝜉) = , 𝜑p+1 (𝜉) = ,
high degree on each element of a fixed partition of the 2 2
√ 𝜉
computational domain. For a one-dimensional problem, 2k − 1 1
such as, for example (2), we split Ω = (a, b) into a set of M 𝜑k (𝜉) = L (s) ds = √
2 ∫−1 k−1 2(2k − 1)
disjoint intervals Ωe , e = 1, … , M, whose end points are
a = x̄ 0 < x̄ 1 < · · · < x̄ M = b. Then we set ×(Lk (𝜉) − Lk−2 (𝜉)), k = 2, … , p

̄ v|Ω ∈ ℙN
VN,M = {v ∈ C0 (Ω)| ∀e = 1, … , M The first two terms ensure C0 continuity of the trial func-
e
tions.
v(a) = v(b) = 0}
For the algebraic realization of SEM, nodal basis func-
tions are those introduced in (13). Being associated with
The approximation of (2) by SEM reads the special set of Legendre–Gauss–Lobatto nodes, they can
be used to generate shape functions. Moreover, LGL inte-
find uN,M ∈ VN,M ∶ a(uN,M , v) = (f , v) ∀v ∈ VN,M (31) gration formulae are used to approximate the integrals that
provide the entries of the problem’s arrays, namely the stiff-
This approach shares the same structure as the p-version ness matrix and the right-hand side vector. This is reflected
of the finite element method (FEM). As in the latter, the by replacing (31) with the more interesting SEM-NI version:
number M of subintervals is frozen, while the local poly-
nomial degree (i.e., indicated by N in the SEM context and ∑
M
∑
M
by p in the FEM context) is increased to improve accuracy. find uN,M ∈ VN,M ∶ aN,Ωe (uN,M , v) = (f , v)N,Ωe
More precisely, if h = (b − a)∕M denotes the constant length e=1 e=1
of each subinterval, one has ∀v ∈ VN,M (33)
N
||u′ − (ΠN,M u)′ ||L2 (a,b) + ||u − ΠN,M u||L2 (a,b) where (u, v)N,M is the Legendre–Gauss–Lobatto inner
h ∑N
product (11) in Ωe , (u, v)N,Ωe = j=0 𝛼je u(xje )v(xje ), with
≤ C(s)hmin(N,s) N −s ||u′ ||H s (a,b) , s ≥ 0 (32)
𝛼je = 𝛼j [(b − a)∕2], xje is the correspondent of xj in Ωe .
Moreover, aN,Ωe (u, v) is the elemental bilinear form.
where ΠN,M is the SEM interpolant. If u is arbitrarily smooth
Still considering the case of the differential operator (5) as
(s large), it is advantageous to keep h fixed and let N → ∞.
an instance, we end up with the following form:
Should the different degree of smoothness suggest the use of
a nonuniform polynomial degree, another upper bound for
the left-hand side of (32) is aN,Ωe (u, v) = (𝛼u′ , v′ )N,Ωe + (𝛽u′ , v)N,Ωe + (𝛾u, v)N,Ωe

∑
M in analogy with the left-hand side of (12).
Ce hmin(Ne ,se ) Ne e ||u′ ||H se (Ωe ) ,
−s
se ≥ 1, The multidimensional case can be addressed by first intro-
e=1 ducing the tensorized basis functions on the parental element
∀e = 1, … , M Ω̂ = (−1, 1)d (d = 2, 3), then mapping basis functions and
nodal points on every element Ωe (now a quadrilateral or
where Ne is the polynomial degree used in the eth element parallelepipedal structure, possibly with curved edges or
Ωe , and H se +1 (Ωe ) is the local smoothness of u in Ωe . surfaces). The functional structure of our problem remains
SEM was first introduced by Patera (1984) for Chebyshev formally the same as in (33), and the kind of error estimate
expansions, then generalized to the Legendre case by Y. that can be achieved is similar. Obviously, this time VN,M is
Maday and A. Patera. Both approaches (SEM and p-version made of globally continuous functions that satisfy homoge-
of FEM) make use of a parental element, say Ω ̂ = (−1, 1), neous Dirichlet boundary data (if any). They are obtained by
on which the basis functions are constructed. However, the joining the elemental functions that are the mapping of the
main difference lies in the way the basis functions are chosen nodal basis functions according to the transformation Te ∶
(and, therefore, in the structure of the corresponding stiffness Ω̂ → Ωe that maps the parental element Ω ̂ into the current
matrix). element Ωe .
12 Spectral Methods

We refer to the seminal chapter by Patera (1984) and to the The mathematical rationale behind the choice of the
books by Bernardi and Maday (1997), Deville et al. (2002), matching condition (36) (rather than a more “natural” condi-
Karniadakis and Sherwin (1999), and Schwab (1998). tion of pointwise continuity at one set of grid nodes on Γ)
becomes clear from the convergence analysis for problem
(34). With this aim, we introduce
8 THE MORTAR METHOD ( )1∕2
||v||∗ ∶= ||v||20,Ω + ||∇v|Ω1 ||20,Ω + ||∇v|Ω2 ||20,Ω (37)
1 2
This method has been introduced by Bernardi et al. (1994)
with the aim of allowing spectral elements having different which is a norm (the “graph” norm) for the Hilbert space
polynomial degrees or being geometrically nonconforming,
and also to allow the coupling of the spectral (element) H∗ ∶= {v ∈ L2 (Ω) | v|Ω1 ∈ H 1 (Ω1 ), v|Ω2 ∈ H 1 (Ω2 )} (38)
method with the finite element method. Its generality,
however, goes beyond these two specific examples. Consider, Owing to the Poincaré inequality, we have that
for the sake of illustration, the Poisson problem with homo-
∑
2
geneous Dirichlet conditions. The idea is to approximate its |∇v𝛿 |2 ≥ 𝛼∗ ||v𝛿 ||2∗ ∀v𝛿 ∈ V𝛿 (39)
weak form (13) by the following discrete problem: i=1
∫ Ωi

∑
M
∑
M from where the discrete problem (34) admits a unique solu-
find u𝛿 ∈ V𝛿 ∶ ∇u𝛿 ⋅ ∇v𝛿 = f v𝛿 ∀v𝛿 ∈ V𝛿 tion by a straightforward application of the Lax–Milgram
i=1
∫ Ωi i=1
∫ Ωi lemma. Furthermore, one can prove the following inequality
(34) for the norm of the error u − u𝛿 :
( )
Here, 𝛿 > 0 is a parameter describing the quality of the 1
||u − u𝛿 ||∗ ≤ 1 + inf ||u − v𝛿 ||∗
discretization, and V𝛿 is a finite dimensional space that 𝛼∗ v𝛿 ∈V𝛿
approximates H01 (Ω) without being contained into C0 (Ω).̄
| 𝜕u |
|∫Γ 𝜕n [w𝛿 ]Γ |
More precisely, V𝛿 is a subspace of the following space:
+
1
sup | | (40)
𝛼∗ w𝛿 ∈V𝛿 ||w𝛿 ||∗
Y𝛿 ∶= {v𝛿 ∈ L2 (Ω) | v𝛿|Ωi ∈ Yi,𝛿 , i = 1, … , M} (35)
where [w𝛿 ]Γ ∶= w(1) 𝛿|Γ
− w(2)
𝛿|Γ
denotes the jump across Γ of a
where, for each i = 1, … , M, Yi,𝛿 is a finite dimensional function w𝛿 ∈ V𝛿 . The approximation error of (34) is there-
subspace of H 1 (Ωi ): it can be either a finite element space, fore bounded (up to a multiplicative constant) by the best
or a polynomial spectral (elements) space. In any case, no approximation error (i.e., the distance between the exact
requirement of compatibility is made for the restriction of solution u and the finite dimensional space V𝛿 ) plus an extra
the functions of Y𝛿 on the element interface Γ. error involving interface jumps. The latter would not appear
Heuristically, the space V𝛿 will be made up of functions in the framework of classical Galerkin approximation (like
belonging to Y𝛿 that satisfy some kind of matching across the SEM), and is the price to pay for the violation of the
Γ. Precisely, assuming for simplicity that there are only two conforming property; that is, for the fact that V𝛿 ⊄ H01 (Ω).
elements, if v𝛿 ∈ V𝛿 and if v(1)
𝛿
∈ Y1,𝛿 , v(2)
𝛿
∈ Y2,𝛿 denote its The error estimate (40) is optimal if each one of the two
restrictions to Ω1 and Ω2 , respectively, then the following terms on the right can be bounded by the norm of local
integral matching conditions should be satisfied for a certain errors arising from the approximations in Ω1 and Ω2 , without
fixed index i (=1 or 2): the presence of terms that combine them in a multiplicative
fashion. In this way, we can take advantage of the local
(v(1) − v(2) )𝜇𝛿(i) = 0 ∀𝜇𝛿(i) ∈ Λ(i) (36) regularity of the exact solution as well as the approximation
∫Γ 𝛿 𝛿 𝛿
properties enjoyed by the local subspaces Yi,𝛿 of H 1 (Ωi ).
To generate a nodal basis for the finite dimensional space
where Λ(i)𝛿
denotes the restriction to Γ of the functions of Yi,𝛿 . V𝛿 , we can proceed as follows. For i = 1, 2, let us denote by
If we take i = 2 in (36), this amounts to letting Ω1 play i the set of nodes in the interior of Ωi , and by Γ(i) the set
the role of master and Ω2 that of slave, and (36) has to be of nodes on Γ, whose cardinality will be indicated by Ni and
intended as the way of generating the value of v(2) 𝛿|Γ
once v(1)
𝛿|Γ
NΓ(i) , respectively. Note that, in general, Γ(1) and Γ(2) can
is available. The alternative way, that is, taking i = 1 in (36) is be totally unrelated.
also admissible. Depending upon the choice of index i made Now, denote by {𝜑(1) k′
}, k′ = 1, … , N1 , the Lagrange func-
in (36), the method will produce different solutions. tions associated with the nodes of 1 ; since they vanish on
 Spectral Methods 13

Γ, they can be extended by 0 in Ω ̄ 2 . These extended functions The discrete problem (34) can also be reformulated as a
are denoted by {𝜑̃ (1)
k′
}, and can be taken as a first set of basis saddle point problem. Its algebraic counterpart is similar to
functions for V𝛿 . (26), that is, it can be formulated as
Symmetrically, we can generate as many basis functions [ ] [ ] [ ]
for V𝛿 as the number of nodes of 2 by extending by 0 in Ω ̄1 A BT u f
=
the Lagrange functions associated with these nodes. These B 0 𝛌 𝟎
new functions are denoted by {𝜑̃ (2)k′′
}, k′′ = 1, … , N2 . where u collects the nodal values of u𝛿 in each subdomain,
Finally, always supposing that Ω1 is the master domain and and 𝛌 collects the nodal values on Γ of the Lagrange multi-
Ω2 its slave, for every Lagrange function {𝜑(1) } in Ω̄ 1, m =
(1)
m,Γ plier 𝜆𝛿 associated with the “constraint” (36). The matrix
1, … , NΓ , we obtain a basis function {𝜑̃ m,Γ } as follows: A is block-diagonal (with one block per subdomain Ωi ),
{ each block corresponding to a problem for the Laplace oper-
𝜑(1) ̄1
in Ω ator with a Dirichlet boundary condition on 𝜕Ωi ∩ 𝜕Ω and
𝜑̃ m,Γ ∶= m,Γ
𝜑̃ (2) ̄2
in Ω a Neumann boundary condition on 𝜕Ωi ∖𝜕Ω. The matrix B
m,Γ
is associated with the bilinear form on the left-hand side of
where (36), that is, it accounts for the interaction between function
(2)
NΓ
∑ jumps and Lagrange multiplier on Γ.
𝜑̃ (2)
m,Γ
∶= 𝜉j 𝜑(2)
j,Γ After elimination of the degrees of freedom internal to the
j=1 subdomains, the method leads to the reduced linear system
(still of a saddle point type)
𝜑(2) are the Lagrange functions in Ω ̄ 2 associated with the
j,Γ [ ] [ ] [ ]
(2)
nodes of Γ , and 𝜉j are unknown coefficients that should S CT uΓ g
= Γ
be determined through the fulfillment of the matching C 0 𝛌 𝟎
equations (36). Precisely, they must satisfy
where the matrix S is block-diagonal, C is a jump operator,
⎛ NΓ
(2)
⎞ uΓ is the set of all nodal values at subdomain interfaces, and
∑
⎜ 𝜉j 𝜑(2) − 𝜑(1) ⎟ 𝜑(2) = 0 ∀l = 1, … , NΓ(2) (41) gΓ is a suitable right-hand side.
∫Γ ⎜ j=1 j,Γ m,Γ ⎟ l,Γ
This system can be regarded as an extension of the Schur
⎝ ⎠
complement system to nonconforming approximation (the
A basis for V𝛿 is therefore provided by the set of all Lagrange multiplier 𝛌 indeed accounts for the nonmatching
functions {𝜑̃ (1)
k′
}, k′ = 1, … , N1 , {𝜑̃ (2)
k′′
}, k′′ = 1, … , N2 , and discretization at subdomain interfaces). In fact, the ith block
{𝜑̃ m,Γ }, m = 1, … , NΓ .(1) of S is the analog of Σi,h , and corresponds to a discretized
Steklov–Poincaré operator on the subdomain Ωi .
Remark. In the mortar method, the interface matching is
achieved through a L2 -interface projection, or, equivalently,
by equating first-order moments, thus involving computation 9 THE SPECTRAL DISCONTINUOUS
of interface integrals. In particular, from equation (36), we GALERKIN METHOD
have to evaluate two different kinds of integrals (take, for
instance, i = 2): An alternative approach to the mortar method is offered
by the discontinuous Galerkin (DG) spectral method. DG
I12 ∶= v(1) 𝜇𝛿(2) , I22 ∶= v(2) 𝜇𝛿(2) methods have gained tremendous popularity in the finite
∫Γ 𝛿 ∫Γ 𝛿
element community in the past three decades (Riviere, 2008;
Arnold et al., 2002; Wohlmuth, 2001) and more recently
The computation of I22 raises no special difficulties, because in the spectral context (Hesthaven and Warburton 2008).
both functions v(2)
𝛿
and 𝜇𝛿(2) live on the same mesh, the The main idea is to allow for the use of piecewise polyno-
one inherited from Ω2 . On the contrary, v(1) 𝛿
and 𝜇𝛿(2) are mials that are discontinuous at the element (or subdomain)
functions defined on different domains, and the computa- interfaces at the cost of inserting into the discrete varia-
tion of integrals like I12 requires proper quadrature rules. tional formulation some extra interface terms. The latter
This process needs to be done with special care, especially typically depend on either jumps and averages of interface
for three-dimensional problems, for which subdomain inter- quantities (the desired solution, the test functions, and their
faces are made up of faces, edges, and vertices; otherwise normal fluxes/stresses) and fulfill a multiple purpose: besides
the overall accuracy of the mortar approximation could be restoring the consistency with the original PDE, they need to
compromised. ensure stability of the resulting numerical approximation.
14 Spectral Methods

For the sake of simplicity, let us consider the Poisson follows: find u𝛿 ∈ V𝛿 ∶
equation (14) (where we set 𝛼 = 1, 𝜷 = 𝟎, 𝛾 = 0 and d = 2)
(
and its corresponding variational formulation (15). ∑
M
∑
Let us introduce the space V𝛿 = {v𝛿 ∈ Y𝛿 ∶ v𝛿 = 𝛼∇u𝛿 ⋅ ∇v𝛿 − [v𝛿 ] ⋅ {{𝛼∇u𝛿 }}
∫ Ωm e∈𝛿)
∫e
0 on 𝜕Ω}, with Y𝛿 being the same as in (35). We denote by m=1

𝛿 the collection of all internal edges (the element interfaces) +𝜏 [u𝛿 ] ⋅ {{𝛼∇v𝛿 }}
∫e
and, on every edge e ∈ 𝛿 , by [v] = v+ n+ + v− n− the jump
across e (v+ and v− being the values of v from the two ∑ ∑M
+ 𝜌|e|−1 [u𝛿 ] ⋅ [v𝛿 ] − 𝜷u𝛿 ⋅ ∇v𝛿 (45)
sides of the edge and n+ , n− the outward normals to the e∈𝛿
∫e ∫
m=1 Ωm
edge), while {{∇v}} = 12 ((∇v)+ + (∇v)− ) is the average of the ∑ ∑
M
gradient on the edge. The DG approximation to the Poisson + {{𝜷u𝛿 }}𝜷 ⋅ [v𝛿 ] + 𝛾u𝛿 v𝛿
∫e ∫Ωm
problem with homogeneous Dirichlet conditions reads: find e∈𝛿 m=1
u𝛿 ∈ V𝛿 ∶ ∑
M
= f v𝛿 ∀v𝛿 ∈ V𝛿
∫ Ωm
∑ M
∑ m=1
∇u𝛿 ⋅ ∇v𝛿 − [v𝛿 ] ⋅ {{∇u𝛿 }}
m=1
∫ Ωm e∈𝛿
∫e where we have used the notation
∑
−𝜏 [u𝛿 ] ⋅ {{∇v𝛿 }} (42) ⎧𝜷u+ if 𝜷 ⋅ n+ > 0
e∈𝛿
∫e ⎪ 𝛿
{{𝜷u𝛿 }}𝜷 = ⎨𝜷u−𝛿 if 𝜷 ⋅ n+ < 0
∑ ∑M
⎪ 𝜷u𝛿 if 𝜷 ⋅ n+ = 0
+ 𝜌|e|−1 [u𝛿 ] ⋅ [v𝛿 ] = f v𝛿 ∀v𝛿 ∈ V𝛿 ⎩
e∈𝛿
∫e ∫
m=1 Ωm
The DG spectral methods are especially suited when
where |e| denotes the length of e, 𝜌 a suitable constant treating pure hyperbolic equations (and systems, including
depending on the local polynomial degrees, and 𝜏 is a scalar and vector conservation laws), thanks to their intrinsic
suitable constant. Equation (42) is called interior penalty capability of approximating discontinuous solutions. We
(IP) DG method. For 𝜏 = 1 the formulation (42) is named consider for the sake of simplicity the one-dimensional
SIPG (symmetric interior penalty Galerkin) (S standing for problem
Symmetric, G for Galerkin). Other special (non-symmetric)
cases are those of NIPG (𝜏 = −1) and IIPG (incomplete inte- ⎧ut + (𝛽u)x + 𝛾u = f , a < x < b, t > 0
rior penalty Galerkin) (𝜏 = 0) (N for Non, I for Incomplete). ⎪
NIPG is stable for every 𝜌 > 0 whereas both SIPG and IIPG ⎨u(a, t) = g(t), t>0 (46)
⎪u(x, 0) = u (x), a<x<b
require a sufficiently large 𝜌 for stability. If we assume that ⎩ 0
the same polynomial degree, say N, is used on all of the
elements, 𝜌 is chosen proportional to N 2 and h denotes the where 𝛽(x) > 0 and 𝛾(x) ≥ 0 are two given functions and g(t),
element diameter, then the DG solution converges to the u0 (x) are two given boundary and initial data, respectively.
exact one according to the error estimate Still using the partition of [a, b] introduced in Section 7
( )s above, and introducing the notation
h
|||u − u𝛿 ||| ≤ C N 1∕2 |u|H s+1 (Ω) , s≥1 (43)
N W𝛿 = {v ∈ L2 (a, b) ∶ v|Ωm ∈ ℙNm , m = 1, … , M}
where
the DG approximation of (46) reads: for all t > 0, find u𝛿 =
( )1∕2 u𝛿 (t) ∈ W𝛿 such that
∑ M
∑
|||v||| = |∇v|2 + 𝜌|e|−1 [ve ]2 (44)
∫
m=1 Ωm e∈𝛿
∫e ∑
M
{(u𝛿,t , v𝛿 )m − (𝛽u𝛿 , v𝛿,x )m + (𝛾u𝛿 , v𝛿 )m }
If instead of the Poisson problem we would consider the m=1
diffusion advection equation ∑
M−1
− u−𝛿 [𝛽v𝛿 ]x̄ m + 𝛽 − u−𝛿 v−𝛿|x=b
−div(𝛼∇u − 𝜷u) + 𝛾u = f in Ω m=1

∑M
with 𝛼 ≥ 𝛼0 > 0 still with homogeneous Dirichlet boundary = (f , v𝛿 )m + 𝛽 + gv+𝛿|x=a ∀v𝛿 ∈ W𝛿 (47)
conditions on 𝜕Ω, (42) should be modified accordingly, as m=1
 Spectral Methods 15

where at every node x̄ m we set v± = lim± v(̄xm + h), [v]x̄ m = Ideally, Hm would be the exact solution to the Riemann
h→0
x̄ problem defined by the conservation law with discontin-
(v+ − v− )(̄x
m ) and (u, v)m = uv.∫x̄ m
m−1 uous data at the node x̄ m , that is, uL = u−𝛿 (̄xm , t) and uR =
Since the test function v𝛿 is also discontinuous, we can
u+𝛿 (̄xm , t). For nonlinear conservation laws, however, the
reduce (47) to a set of M relations that hold on the M
Riemann problem is usually only solved approximately
elements. More precisely, setting u𝛿(m) = u𝛿 |Ωm whenever a
(without polluting the solution accuracy), see for example
potential confusion exists, we find for every m = 1, … , M
Toro (2009), and we denote the resulting numerical flux by
(u𝛿,t , vNm )m − (𝛽u𝛿 , vNm )m + 𝛽 − u𝛿(m) vNm|̄x
m
Hm (u𝛿 (t)) = H(u−𝛿 (̄xm , t), u+𝛿 (̄xm , t)) (52)
−𝛽 + u𝛿(m−1) vNm|̄x = (f , vNm )m ∀vNm ∈ ℙNm (48)
Furthermore, we set u−𝛿 (̄x0 , t) = g(t), which is the inflow
m−1

with u(0) = g, and u𝛿 |t=0 = u0𝛿 ∈ W𝛿 (a suitable approxima- datum; on the other hand, the boundary flux HM (u𝛿 ) has to
𝛿|̄x0
tion of u0 ). be properly designed to avoid spurious reflections at x = b.
When exact integration elementwise is replaced by numer- Of all the possible types of numerical fluxes, the pref-
ical Gauss-Lobatto integration, (48) can be given a nodal erence is for one that is monotone and yields a monotone
interpretation as a generalized collocation method. The scheme when Nm = 0 (i. e., when piecewise constant poly-
above DG spectral method can be easily extended to the nomials are used), namely, a scheme for which the constant
case of linear systems. value u𝛿(m) at the new time-level tn+1 depends on u𝛿(m−1) , u𝛿(m) ,
Time-discretization can be implemented by high-order and u𝛿(m+1) at the previous time-level tn through an increasing
explicit Runge–Kutta (RK) methods, in which case the monotone function. This preference rests on the funda-
time-step must undergo a stability restriction of the form mental result, attributed to Kuznetsov (1976), that monotone
Δt ≤ C(q)h∕k(N), where q is the RK order and k(N) is an schemes are stable and converge to the solution with nonde-
affine function of the polynomial degree, see Chapter 8 in creasing entropy. Notable examples are the Lax-Friedrichs
Canuto et al. (2007) and Hesthaven and Warburton (2008) flux, the Roe flux, and the Engquist-Osher flux (LeVeque,
for details and analysis. 2002; Toro, 2009). The Lax–Friedrichs flux and the Roe flux
We consider now the case of a nonlinear scalar conserva- are the ones most widely used in applications of SDGM.
tion law For a system of nonlinear conservation laws
𝜕u 𝜕 𝜕q 𝜕
+ (u) = 0 , a<x <b, t>0 (49) +  (q) = 𝟎 , a<x <b, t>0 (53)
𝜕t 𝜕x 𝜕t 𝜕x
where  (u) is a nonlinear function of u, with the initial where q(x, t) ∈ ℝq , q ≥ 2, and  (q) is a nonlinear vector
condition u|t=0 = u0 and boundary condition u|x=a = g (we function called the flux, the derivation of SDGM follows in
are assuming that x = a is the inflow point for all t > 0). We a similar way; however, the form of numerical fluxes is now
discretize this equation by the SDGM (spectral discontinuous more involved (Canuto et al., 2006, 2007).
Galerkin method) as follows, using the same notation as in
the previous section. For every t > 0, we look for u𝛿 (t) ∈ W𝛿 Remark. All results cited in this note can be recovered
such that, for all m = 1, … , M, from the books and general chapters that are quoted in the
( ) ( ) references.
𝜕u𝛿 𝜕v
, v𝛿 −  (u𝛿 ), 𝛿 + Hm (u𝛿 )v𝛿 (̄xm )
𝜕t m 𝜕x m
−Hm−1 (u𝛿 )v𝛿 (̄xm−1 ) = 0 ∀v𝛿 ∈ ℙNm (50)

The function Hm denotes the numerical flux at the node x̄ m . REFERENCES

This formulation may be derived from the weak form of
the exact problem (49) in which we replace the exact flux by Arnold DA, Brezzi F, Cockburn B and Marini LD. Unified analysis
the numerical flux. Note that counter-integrating by parts the of discontinuous Galerkin methods for elliptic problems. SIAM J.
second term in (50) yields, for m = 0, … , M − 1, Numer. Anal. 2001/02; 39(5):1749–1779.
( ) Batcho PF and Karniadakis GE. Generalized Stokes eigenfunctions:
𝜕u𝛿 𝜕 a new trial basis for the solution of incompressible Navier–Stokes
+ (u𝛿 ), v𝛿 + (Hm (u𝛿 ) −  (u𝛿(m) ))v𝛿 |x̄ m
𝜕t 𝜕x m equations. J. Comput. Phys. 1994; 115:121–146.
Bernardi C and Maday Y. Approximations Spectrales de Problèmes
+ ( (u𝛿(m) ) − Hm−1 (u𝛿 ))v𝛿 |x̄ m−1 = 0 (51) aux Limites Elliptiques. Springer-Verlag: Paris, 1992.
16 Spectral Methods

Bernardi C and Maday Y. Spectral methods. In Handbook of Numer- in CFD, Murman EM and Abarbanel SS (eds). Birkhäuser: Boston,
ical Analysis, Techniques of Scientific Computing, vol. V, Ciarlet PJ MA, 1984; 357–375.
and Lions JL (eds). North Holland: Amsterdam, 1997; 209–486. Guermond JL, Minev P and Shen J. An overview of projection
Bernardi C, Maday Y and Patera AT. A new nonconforming methods for viscous incompressible flows. Comput. Methods Appl.
approach to domain decomposition: the mortar element method. Mech. Eng. 2006; 195:6011–6045.
In Nonlinear Partial Differential Equations and their Applications, Guo B-Y. Spectral Methods and their Applications. World Scientific
Collège de France Seminar, vol. XI, Brezis H and Lions JL (eds). Publishing: Singapore, 1998.
Longman Group: Harlow, 1994; 13–51.
Hesthaven JS and Warburton T. Nodal Discontinuous Galerkin
Boyd JP. Chebyshev and Fourier Spectral Methods. Dover: New Methods. Algorithms, Analysis, and Applications, Texts in Applied
York, 2001. Mathematics, vol. 54. Springer-Verlag: New York, 2008.
Canuto C and van Kemenade V. Bubble stabilized spectral methods Karniadakis GE, Israeli M and Orszag SA. High-order splitting
for the incompressible Navier–Stokes equations. Comput. Methods methods for incompressible Navier–Stokes equations. J. Comput.
Appl. Mech. Eng. 1996; 135:35–61. Phys. 1991; 97:414–443.
Canuto C, Hussaini MY, Quarteroni A and Zang TA. Spectral Karniadakis GE and Sherwin SJ. Spectral hp Element Method for
Methods in Fluid Dynamics. Springer-Verlag: Berlin, Heidelberg, CFD. Oxford University Press: New York, 1999.
1988.
Kuznetsov NN. The accuracy of certain approximate methods for
Canuto C, Hussaini YM, Quarteroni A and Zang TA. Spectral the computation of weak solutions of a first order quasilinear
Methods. Fundamentals in Single Domains. Springer-Verlag: equation. (Russian) Ž. Vyčisl. Mat. i Mat. Fiz. 1976; 16:1489–1502
Heidelberg, 2006. (in Russian).
Canuto C, Hussaini YM, Quarteroni A and Zang TA. Spectral LeVeque RJ. Finite Volume Methods for Hyperbolic Problems.
Methods. Evolution to Complex Geometries and Applications to Cambridge University Press: New York, 2002.
Fluid Dynamics. Springer-Verlag: Heidelberg, 2007.
Patera AT. A spectral element method for fluid dynamics: laminar
Deville MO, Fischer PF and Mund EH. High-Order Methods flow in a channel expansion. J. Comput. Phys. 1984; 54:468–488.
for Incompressible Fluid Flow. Cambridge University Press:
Peyret R. Spectral Methods for Incompressible Viscous Flow.
Cambridge, 2002.
Springer-Verlag: New York, 2002.
Dubiner M. Spectral methods on triangles and other domains. J. Sci.
Quarteroni A, Saleri F and Veneziani A. Analysis of the Yosida
Comput. 1991; 6:345–390.
method for the incompressible Navier-Stokes equations. J. Math.
Fischer PF. and Mullen JS. Filter-based stabilization of spectral Pures Appl. 1999; 78(5):473–503.
element methods. C.R. Acad. Sci. Sér. I 2001; 332:265–270.
Quarteroni A, Saleri F and Veneziani A. Factorization methods for
Fornberg BA. Practical Guide to Pseudospectral Methods. the numerical approximation of Navier-Stokes equations. Comput.
Cambridge University Press: Cambridge, 1996. Methods Appl. Mech. Eng. 2000; 188(1-3):505–526.
Funaro D. Spectral Elements for Transport-Dominated Equations. Rivière B. Discontinuous Galerkin Methods for Solving Elliptic and
Springer-Verlag: Berlin, Heidelberg, 1997. Parabolic Equations: Theory and Implementation, Frontiers in
Gervasio P. Convergence analysis of high order algebraic fractional Applied Mathematics, vol. 35. Society for Industrial and Applied
step schemes for time-dependent Stokes equations. SIAM J. Numer. Mathematics (SIAM): Philadelphia, PA, 2008.
Anal. 2008; 46:1682–1703. Schwab Ch. p- and hp- Finite Element Methods. Oxford Science
Gervasio P and Saleri F. Algebraic fractional-step schemes for Publications: New York, 1998.
time-dependent incompressible Navier-Stokes equations. J. Sci. Shen J, Tang T and Wang LL. Spectral Methods – Algorithms,
Comput. 2006; 27(1-3):257–269. Analysis and Applications. Springer-Verlag: Heidelberg, 2011.
Gervasio P, Saleri F and Veneziani A. Algebraic fractional Tadmor E. Approximate solutions of nonlinear conservation laws.
step schemes with spectral methods for the incompressible In Advanced Numerical Approximations of Nonlinear Hyperbolic
Navier-Stokes equations. J. Comput. Phys. 2006; 214(1):347–365. Equations, Quarteroni A (eds). Springer-Verlag: Berlin, 1998;
Gottlieb D and Hesthaven JS. Spectral methods for hyperbolic prob- 1–150.
lems. J. Comput. Appl. Math. 2001; 128:83–131. Toro EF. Riemann Solvers and Numerical Methods for Fluid
Gottlieb D and Orszag SA. Numerical Analysis for Spectral Methods: Dynamics. Springer-Verlag: Berlin, 2009.
Theory and Applications. SIAM-CBMS: Philadelphia, PA, 1977. Wohlmuth B. Discretization Methods and Iterative Solvers Based on
Gottlieb D and Shu C-W. On the Gibbs phenomenon and its resolu- Domain Decomposition. Springer-Verlag: Heidelberg, 2001.
tion. SIAM Rev. 1997; 39:644–668.
Gottlieb D and Tadmor E. Recovering pointwise values of discontin-
uous data with spectral accuracy. In Progress and Supercomputing