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)

CHAPTER 22. DOMAIN DECOMPOSITION METHODS

482

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

(22.13a)

L hLj = 0 & hLj (dj’1 ) = 1 & hLj (dj ) = 0

[homogeneous solutions]

(22.13b)

L hRj = 0 & hRj (dj’1 ) = 0 & hRj (dj ) = 1

22.4. WEAK COUPLING OF ELEMENTAL SOLUTIONS 483

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

(22.14)

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

(22.16)

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

CHAPTER 22. DOMAIN DECOMPOSITION METHODS

484

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

1

2

I(v) ≡ p(x) [vx ] + q(x) v 2 ’ 2 v f (x) dx (22.20)

’1

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)

1

(r, s) ≡ (22.23)

r(x) s(x) dx

’1

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

22.6. CHOICE OF BASIS & GRID 485

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

equations

da

(22.25)

M = ...

dt

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

CHAPTER 22. DOMAIN DECOMPOSITION METHODS

486

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.

22.8. MATRIX INVERSION 487

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