. 99
( 136 .)


An important theme is that solutions on different subdomains are connected only through
the interface matching conditions, (22.5). Almost all spectral element methods exploit this
weak element-to-element coupling.
To illustrate how and why, consider a linear, second order different equation, i. e. (22.1)
and (22.2). Instead of just two subintervals as in the previous section, however, divide the
interval x ∈ [a, b] into M subdomains. Let dj denote the boundary between element (j-1)
and element j and de¬ne

Uj ≡ u(dj ) (22.11)

Figure 22.1: Schematic of the interval [a, b] split into subdomains. The element boundaries
are denoted by dj ; the values of u(x) at these boundaries are denoted by Uj . The solution
on the j-th element, uj (x), is the sum of a particular integral pj (x) plus the weighted con-
tributions of the homogeneous solutions hLj (x) and hRj (x). The schematic illustrates the
boundary conditions satis¬ed by each component.

The condition of continuity of u(x) at the domain boundaries may be enforced implicitly by
demanding that both uj’1 (x) and uj (x) equal Uj at x = dj . These de¬nitions are illustrated
in Fig. 22.1.
The solution on the j-th element, uj (x), can always be written as the sum of a particular
integral pj (x) plus two homogeneous solutions. One always has the freedom to choose
the particular integral so that it vanishes at both elemental boundaries. One may similarly
choose the homogeneous solutions so that one ” call it hLj (x) ” is equal to one at the
left domain boundary and zero at the other while hRj (x) is non-zero at the right boundary
but not at the left wall. If the differential equation is L u = f , these three components are
de¬ned by
[particular integral] (22.12)
L pj (x) = f & pj (dj’1 ) = pj (dj ) = 0

L hLj = 0 & hLj (dj’1 ) = 1 & hLj (dj ) = 0
[homogeneous solutions]
L hRj = 0 & hRj (dj’1 ) = 0 & hRj (dj ) = 1

These components are schematically shown in Fig. 22.1.
It follows that the general solution to the differential equation can be written ” without
approximation ” as

uj (x) = pj (x) + Uj’1 hLj (x) + Uj hRj (x)

This decomposition is useful because one can numerically compute all the particular in-
tegrals and homogeneous solutions intra-elementally ” that is, calculate each pj (x) or
hLj (x) or hRj (x) independently of the solutions on all the other elements. The boundary
value problems de¬ned by (22.12) and (22.13) are completely uncoupled.
Unfortunately, the elemental solutions uj (x) are not completely speci¬ed because (22.14)
contains (M + 1) unknown parameters, Uj , the values of u(x) at the domain boundaries.
The two end values are determined by the boundary conditions for the original problem:

’ ’ (22.15)
u(a) = ± U0 = ± & u(b) = β UM = β

The remaining domain boundary values of u(x) are determined by the requirement of con-
tinuity of the ¬rst derivative at each of the (M ’ 1) interior domain walls. This gives, using
primes to denote differentiation,

hLj Uj’1 + hRj ’ hL, j+1 Uj ’ hR, j+1 Uj+1 = pj+1 ’ pj j = 1, . . . , M ’ 1

where all functions are evaluated at x = dj . The reason that this matrix system is tridi-
agonal, soluble in only O(8M ) multiplications and additions, is that only four different
homogeneous solutions, two for element (j ’ 1) and two for element j, enter the conti-
nuity condition at x = dj . There are only three unknowns in each row of (22.16) because
Uj [≡ u(dj )] is the coef¬cient of both hRj (x) and hL,j+1 (x).
Thus, the cost of coupling elements together via continuity of u(x) and du/dx is trivial.
The real work is in computing the particular and homogenous solutions on each element.
This may be done by whatever method ” indeed, by whatever subroutine ” one would
have employed to solve the boundary value problem by a single Chebyshev series on the
whole interval.
To compare the global expansion with the subdomain procedure, assume that the global
series has N degrees of freedom while N/M spectral coef¬cients are used on each subdo-
main. If the cost to solve a boundary value problem with ν Chebyshev coef¬cients is O(ν 3 )
operations and O(ν 2 ) storage ” true if the pseudospectral method is combined with Gaus-
sian elimination ” then the relative costs are

O(N 3 ) O(N 2 )
Global: ops. store (22.17a)
M O([N/M ]3 ) ops. & M O([N/M ]2 ) store
Subdomain: (22.17b)

The spectral element method is cheaper by a factor of

1/M 2 operations 1/M storage (22.18)

These are enormous factors if M 1! Of course, to retain the high accuracy of spectral
elements, one must keep N/M ≥ 6 (or so). Furthermore, more ef¬cient methods for solv-
ing the boundary value problems, such as the preconditioned iterations discussed earlier,
would reduce the competitive advantages of the subdomain strategy versus the use of a
single global expansion. Nevertheless, even if (22.18) is a little optimistic, it still dramatizes
the potential of “spectral substructuring”.

The key to ef¬ciency is that most of the work is reduced to solving boundary value
problems intra-elementally, i. e. (22.12) & (22.13). If we assumed M different expansions
on each of the M elements and then lumped all the collocation, boundary, and continuity
conditions into a single matrix of dimension N , the subdomain method would be just as
expensive as using a single global expansion ” and not as accurate. An ef¬cient spectral
element code is one that solves the global problem in two stages. The ¬rst is to solve
many uncoupled smaller problems within each element. The second stage is to couple the
elemental solutions together into a single, continuous global solution.

22.5 Variational Principles
The heart of the variational formulation is to multiply the differential equation residual by
a test function and then integrate-by-parts so that a term like φi (x) φj,xx (x) becomes trans-
formed into the product of ¬rst derivatives. When the basis functions are very low order
polynomials, this trick is important. For example, a popular ¬nite element basis is com-
posed of piecewise linear polynomials, the so-called “tent” or “chapeau” functions. Since
the second derivative of these functions is identically zero, it is not possible to solve second
order equations with “tent” functions by collocation. After integration-by-parts, however,
the piecewise linear polynomials give second order accuracy.
For example, the self-adjoint boundary value problem

’ [p(x) ux ]x + q(x) u = f (x) (22.19)
u(’1) = u(1) = 0

is equivalent to minimizing the functional
I(v) ≡ p(x) [vx ] + q(x) v 2 ’ 2 v f (x) dx (22.20)

over the class of functions v(x) such that v(’1) = v(1) = 0. In practice, (22.20) is minimized
over a ¬nite dimensional space of basis functions, {φj , j = 1, . . . , N }. The coef¬cients of
the series for u(x) are determined by solving the matrix equation Ha = F where

Hij ≡ (φi,x , p φj,x ) + (φi , q φj ) (22.21)

Fi ≡ (φi , f ) (22.22)
(r, s) ≡ (22.23)
r(x) s(x) dx

Although (22.21) is justi¬ed by the calculus of variations, it is obvious that (22.21) is iden-
tical with the usual Galerkin matrix element except for use of the integration-by-parts
1 1
φi ’ [ p φj, x ]x dx ’ (22.24)
φi, x p φj, x
’1 ’1

(The boundary term, (φi (1) p(1) φj, x (1) ’ φi (’1) p(’1) φj, x (’1)), is zero because the basis
functions satisfy the boundary conditions (22.19).)
As noted above, the variational form is essential to successful use of piecewise func-
tions of ¬rst degree to solve differential equations of second order. As the degree of the
polynomials increases, however, the advantage of the variational form becomes smaller
and smaller. This is the reason that the integrated-by-parts form (22.21) is rarely used with

global spectral methods; in comparison to the standard Galerkin procedure, the variational
matrix is not worth the bother.
The domain decomposition/variational principle strategy allows considerable ¬‚exi-
bility. Delves and Hall (1979) add additional terms to the variational functional (22.20)
(“global elements”) so that the basis functions need not even be continous at the element
boundaries. As the resolution is increased, however, the discontinuities in u decrease
rapidly, and the order-of-magnitude of these jumps is never larger than that of the maxi-
mum pointwise error.
Patera (1984), who dubbed his method “spectral elements”, takes a middle ground. His
basis functions belong to C 0 , that is, are continuous. However, continuity of the ¬rst and
higher derivatives at inter-element boundaries is only approximate.
Phillips and Davies (1988), who also use the term “spectral elements”, employ basis
functions that are at least C 1 . Their basis explicitly satis¬es all boundary conditions and
also obeys the element-interface conditions that are required to specify a unique solution
to the (undiscretized) problem.
The higher the degree of continuity, the smaller the matrix which discretizes the differ-
ential equation. With Delves™ and Hall™s “global elements”, the basis has extra degrees of
freedom ” the freedom to be discontinuous on inter-element boundaries ” which are ex-
plicitly eliminated in Phillips™ and Davies™ method. Furthermore, the variational functional
is simpler when the basis functions are C 1 because the functional lacks the extra terms that
are needed to enforce continuity of u(x) if the basis functions are not continuous. Thus, in
principle, using basis functions with continuous ¬rst derivatives would be optimum. In
practice, computing and manipulating basis functions in C 1 may be very dif¬cult if the ge-
ometry is complicated, so the choice of the degree of continuity of the basis set must be left
to individual discretion. However, the C 0 basis has become the most popular (Karniadakis
and Sherwin, 1999).

22.6 Choice of Basis & Grid: Cardinal versus Orthogonal
Polynomial, Chebyshev versus Legendre,
Interior versus Extrema-and-Endpoints Grid
Phillips and Davies (1988) use ordinary Chebyshev polynomials as the intra-element basis,
but it is far more common to use the Chebyshev or Legendre cardinal functions. The reasons
are similar to those that have made the pseudospectral method popular for global spectral
methods. The issue of Tn (x) versus cardinal function is not very important for boundary
value problems if the matrix problem is solved by direct methods, but ” as true of global
algorithms ” it is of paramount important when time-integration or iteration is applied.
In older books on ¬nite elements such as Strang and Fix (1973), only a non-cardinal basis
was used. The drawback is that because the piecewise polynomials are not orthogonal,
the time derivatives for the differential equation du/dt = f (u, x, t) were multiplied, in
discretized form, by the so-called “mass matrix” M. That is, the ¬nite element method
converted a partial differential equation into the coupled system of ordinary differential

M = ...

where the ellipsis ( . . . ) denotes an unspeci¬ed right-hand side and where the elements of

M are
Mij ≡ (φi , φj ) “mass matrix” (22.26)
In a non-cardinal basis, ¬nite element methods are always implicit in the sense that one
must backsolve using the LU -factorization of the “mass matrix” at every time step.
If one shifts to a cardinal function basis and applies collocation instead of integration,
however, the new mass matrix is diagonal. This pseudospectral approach is essential for
ef¬cient, explicit time-integration, just as for global Chebyshev methods.
“Spectral elements” began with a Chebyshev basis, but Legendre polynomials have
become popular, too. One reason is that the usual drawbacks of Legendre polynomials ”
collocation points that are not known in closed form, poorer accuracy near the boundaries,
and so on ” are less important for polynomials of moderate degree than when N is large.
A more compelling reason is that the variational principle employs inner products with a
weight function of unity. Because the Chebyshev polynomials are orthogonal with respect

to a weight function of 1/ 1 ’ x2 , it is more awkward to use them instead of Legendre
polynomials while still remaining faithful to the variational principle.
Either basis, however, works just ¬ne; the difference between them seems to have more
to do with mathematical tidiness than with ef¬ciency. The variational formalism is popular,
but collocation with either Chebyshev or Legendre polynomials is successful, too.
As noted earlier, there are two standard collocation grids for Chebyshev polynomials:
(i) the “roots” or “Gauss-Chebyshev” grid, whose points are the zeros of TN (x) and (ii) the
“extrema-and-endpoints” or “Gauss-Chebyshev-Lobatto” grid, which includes the roots of
the derivative of TN (x) plus the endpoints ±1. As shown in Chapter 2, both have equally
good theoretical properties and the choice between them is usually a matter of preference.
With spectral elements, however, the “extrema-and-endpoints” or “Lobatto” grid is
much more popular than the roots grid. The reason is that one must explicitly match so-
lutions along inter-element boundaries. In the variational formalism, it is awkward to im-
pose continuity on the basis sets, apply the isoparameteric mapping described in Sec. 10,
and so on if none of the grid points coincides with the element walls. However, Schumack,
Schultz, and Boyd (1989) have found that the roots grid is just as accurate as the Lobatto
grid for computing Stokes™ ¬‚ows.
Nevertheless, the Legendre basis with Lobatto grid has become the canonical choice for
spectral elements. The cardinal functions, derivative matrix and grid points are all given
in Appendix F, Sec. 10.

22.7 Patching versus Variational Formalism
Canuto et al. (1988) show that the differences between the patching & variational ap-
proaches to spectral elements are small. Indeed, the variational principle can be expressed
in terms of patching and vice versa if the two formalisms are appropriately generalized.
For example, the variational method is equivalent to patching if the condition of conti-
nuity of du/dx at inter-element boundaries, (22.5b), is replaced by
u1,x (d) ’ u2,x (d) + w0 R1 (d) + wN2 R2 (d) (22.27)
where R1 and R2 are the differential equation residuals in elements 1 and 2 and where w0
and wN2 are the ¬rst and last weights of the Gauss-Lobatto quadrature. Since the quadra-
ture weights are O(1/Ni2 ) and the residuals are very small if the numerical solution is ac-
curate, it follows that (22.27) is only a small perturbation of (22.5b). Further, the magnitude
of the residual terms in (22.27) is always the same order of magnitude as the error in the
numerical solution: exponentially small.

Thus, whether one applies (22.5b) or (22.27) is irrelevant to accuracy. It is obvious, how-
ever, that (22.27) is harder to program and also harder to explain and justify. In practice,
one would always either apply patching in the standard form (22.5b) or apply the varia-
tional principle as explained in Sec. 5. The only use of (22.27) is conceptual: it shows that


. 99
( 136 .)