shown that for relatively weak singularities like that in (6.50), the corner branch points

may be irrelevant unless one is interested in very high accuracy (greater than ¬ve decimal

places). Note that the function in (6.50) is everywhere bounded; only its higher derivatives

are in¬nite at r = 0, and it is in this sense that the singularity is “weak”.

For stronger singularities, such as cracks in a solid material or obtuse angles such as

that in the oft-studied “L-shaped domain eigenvalue problem”, one needs special methods

and tricks. Changes-of-coordinate and “singularity-subtraction” are explained in Chapter

16 and “global sectorial elements” in Chapter 22. The generic recommendation is to ignore

the singularities unless either (i) one has prior knowledge that u(x, y) is discontinuous or

has other strongly pathological behavior or (ii) poor convergence and rapid variation of

the numerical solution near the corners suggests a posteriori that the solution is strongly

singular.

The vagueness of this advice is unsatisfactory. However, it is a fact of life that shock

waves and other types of discontinuities always require special treatment whether the ba-

sic algorithm is spectral, ¬nite difference, or ¬nite element. If the solution is suf¬ciently

well-behaved so that it can be computed by ¬nite differences without tricks, then it can be

computed by unmodi¬ed spectral methods, too, even with corners.

6.13 Matrix methods

Because inverting the (dense) pseudospectral matrix is usually the rate-determining step,

we devote a whole chapter (Chapter 15) to good strategies for this. A couple of preliminary

remarks are in order.

First, if the total number of basis functions is small (N < 200), use Gaussian elimination.

The LU factorization of a 100 — 100 matrix requires less than a third of a second on a

personal computer with a mega¬‚op or better execution rate, so there is no excuse for being

fancy. Even Boyd (1986c), which solves a two-dimensional nonlinear problem with 200

basis functions, used Gaussian elimination on a lowly IBM AT. For over-determined or

ill-conditioned systems, a QR factorization should be used instead.

For large, multi-dimensional problems, Gaussian elimination may be too expensive. In

this case, the iterations and multi-grid strategies described in Chapter 15 are very impor-

tant, perhaps essential. They are almost never justi¬ed for one-dimensional problems. No

one is awarded style points in numerical analysis: Keep-It-Simple,Stupid and use direct,

non-iterative Gaussian elimination or Householder™s method whenever possible.

6.14 Checking

The most reliable way to verify the correctness of a spectral code is to evaluate the residual

by an independent algorithm. The recommendation is to use ¬nite differences with a very

small grid spacing h. Observe that since the pseudospectral approximation uN (x) is a

CHAPTER 6. PSEUDOSPECTRAL METHODS FOR BVPS

122

series expansion, we can evaluate it and its derivatives for arbitrary x, not merely at values

of x on the pseudospectral interpolation grid.

Thus, to verify that uN (x) is a good approximation to the solution to uxx +q(x)u = f (x),

evaluate

[uN (x + h) ’ 2 uN (x) + uN (x ’ h)]

Rfd (x) ≡ + q(x) uN (x) ’ f (x) (6.51)

h2

where h ∼ O(1/1000). The ¬nite difference residual should be small over the whole com-

putational interval, and especially small [equal to the ¬nite differencing error] at the N

points of the pseudospectral grid.

There is one mild caveat: for high order differential equations, the residual may be

O(1) even when uN (x) is a terri¬cally good approximation to u(x). For example, suppose

a differential equation of sixth order has an exact solution which is a trigonometric cosine

polynomial of degree 30. Suppose that the pseudospectral method exactly computes all of

the Fourier coef¬cients except for a30 , which is in error by 1—10’6 because of roundoff. The

maximum pointwise error of the approximation to u(x) is everywhere less than or equal to

1 — 10’6 [recall that all cosines are bounded by 1]. However, the sixth derivative of u30 (x)

multiplies each coef¬cient of u30 (x) by n6 . The error in u30,xxxxxx (x) is 306 — 10’6 = 729!

The residual will be large compared to one even though the pseudospectral solution is

terri¬cally accurate.

It follows that calculating the residual of a numerical approximation is a very conser-

vative way to estimate its error. By “conservative”, we mean that the error in uN (x) will

usually be much smaller than the residual. If k is the order of the differential equation,

then Rf d (x) is typically O(N k ) [Fourier] or O(N 2k ) [Chebyshev] times the error in uN (x)

(Boyd, 1990c).

To verify the adequacy of a given truncation N , i. e. to estimate the error in u(x), there

are two good tactics. The reliable-but-slow method is to repeat the calculation with a differ-

ent value of N . The error in the run with larger N is smaller (probably much smaller) than

the difference between the two computations. The quick-and-dirty method is to simply

make a log-linear plot of |an | versus n. The error will be roughly the order of magnitude

of the highest computed Chebyshev coef¬cient (or the magnitude of the “envelope” of the

coef¬cients (Chapter 2) at n = N , if the an oscillate with n as they decay).

De¬nition 16 (IDIOT) Anyone who publishes a calculation without checking it against an iden-

tical computation with smaller N OR without evaluating the residual of the pseudospectral approx-

imation via ¬nite differences is an IDIOT.

The author™s apologies to those who are annoyed at being told the obvious. However,

I have spent an amazing amount of time persuading students to avoid the sins in this

de¬nition.

Gresho et al.(1993) describe an embarrassment in the published literature. The ¬‚ow

over a backward-facing step was used as a benchmark for intercomparison of codes at

a 1991 minisymposium. Except for the spectral element code, all programs gave steady,

stable ¬‚ow . Despite the disagreement, the spectral element team published their unsteady,

¬‚uctuating solutions anyway in the Journal of Fluid Mechanics. This prompted the Gresho et

al.(1993) team to reexamine the problem with multiple algorithms at very high resolution.

They found that spectral elements did indeed give unsteady ¬‚ow when the number of

elements was small and the degree N of the polynomials in each element was moderate.

However, when the resolution was increased beyond that of the Journal of Fluid Mechanics

paper, either by using more subdomains or very high N within each subdomain, the ¬‚ow

6.15. SUMMARY 123

Table 6.1: Checking Strategies

1. Geometric decrease of coef¬cients an with n. (least reliable)

2. Finite difference residual (most conservative)

3. Varying the truncation N

became steady and in good agreement with the other ¬nite difference and ¬nite element

simulations. It is easy for even experienced and knowledgeable people to be fooled!

One must watch the convergence of a numerical code as carefully as a father watching

his four year old play close to a busy road.

The quick-and-not-so reliable tactic for estimating error is to simply look at the rate of

decrease of the spectral coef¬cients (for a non-cardinal basis). If aN is only 10 or 100 times

smaller than a0 , then the calculation is certainly accurate to no more than 1 decimal place

and may be completely in error. On the other hand, if aN is O(10’8 )a0 , then the solution is

probably accurate to many decimal places.

One must be cautious, however, because this is an optimistic means of estimating error.

To solve a differential equation of order k, one may need good accuracy in approximating

the k-th derivative, and this is hard to estimate from the series for u(x) itself. (A rough rule-

of-thumb is to multiply an by nk for Fourier or Chebyshev series.) A quick inspection of the

Chebyshev coef¬cients is reliable for rejecting calculations whose accuracy is poor because

N is too small. Coef¬cient-inspection is suggestive but not reliable for certifying that a

computation is highly accurate. When inspection suggests that a run might be accurate,

repeat with different N and be sure!

6.15 Summary

I assert that a well-organized pseudospectral code to solve a differential equation is only

a little harder to write than its ¬nite difference counterpart. Table 6.2 illustrates a sample

program to support this contention. For clarity, comments are appended to the end of

many lines (although not allowed in FORTRAN) and also deleted some type conversion

statements.

Excluding comments, the output block, and the user-supplied functions which would

have to be included in a ¬nite difference code, too, the program contains only 58 lines. A

two-dimensional boundary value solver would have an extra subroutine to compute the

auxiliary indexing arrays, and the basis subroutine would be slightly longer. Nevertheless,

one cannot escape the conclusion that solving BVP™s using spectral methods is simply not

very dif¬cult.

Real world nonlinear eigenvalue solvers, of course, may be rather complicated. How-

ever, the complexity comes from path-following algorithms (Appendices C & D), subrou-

tines to initialize the continuation, blocks of code to compute complicated equation coef-

¬cients and so on. These blocks would have to be included in the corresponding ¬nite

difference or ¬nite element codes as well.

It follows that approximating the derivatives in the differential equation is almost never

the hard part of writing a scienti¬c code. Given that this is so, one might as well approxi-

mate the derivatives well instead of badly. One should solve the problem to high accuracy

so that one can then forget the numerics, banish truncation worries from the mind, and

concentrate on the physics of the solution. Pseudospectral methods do just this.

CHAPTER 6. PSEUDOSPECTRAL METHODS FOR BVPS

124

Table 6.2: A sample FORTRAN program for solving a two-point boundary value problem:

x ∈ [’1, 1]

d2 (x) uxx + d1 (x) ux + d0 (x) u = f (x)

u(’1) = ± u(1) = β

DIMENSION XI(20),APHI(20),G(20),H(20,20),UGRAPH(101)

COMMON/BASIS/PHI(20),PHIX(20),PHIXX(20)

PI = 3.14159...

C SPECIFY PARAMETERS

[ u(-1) ]

ALPHA = 1.

[ u(1) ]

BETA = -0.5

[No. of Chebyshev polynomials]

N = 22

[No. of basis functions φj (x); φj (±1) = 0]

NBASIS = N - 2

C COMPUTE THE INTERIOR COLLOCATION POINTS AND THE FORCING VECTOR G.

DO 100 I=1,NBASIS

XI(I)= COS(PI*I/(NBASIS+1))

X = XI(I)

B = ALPHA*(1-X)/2.+BETA*(1+X)/2. [Function added to v(x) so that

v(x) ≡ u(x) ’ B(x) satis¬es

homogeneous boundary conditions.]

[x-derivative of B(x)]

BX = (-ALPHA + BETA)/2.

[modi¬ed inhomogeneous term]

G(I)=F(X) - D0(X)*B - D1(X)*BX

100 CONTINUE

C COMPUTE THE SQUARE MATRIX

DO 200 I=1,NBASIS

X = XI(I)

CALL BASIS(X,NBASIS) [Compute all N basis functions & derivatives.

Results are returned via arrays in COMMON.]

DD0 = D0(X)

[These three lines avoid needless calls to

DD1 = D1(X)

D0(x), etc., inside the J loop.]

DD2 = D2(X)

DO 200 J=1,NBASIS

H(I,J)=DD2*PHIXX(J)+DD1*PHIX(J)+DD0*PHI(J)

200 CONTINUE

C CALL A LINEAR ALGEBRA ROUTINE TO SOLVE THE MATRIX EQUATION

H * APHI = G

CALL LINSOLV(H,G,APHI,NBASIS)

C The array APHI now contains the NBASIS coefficients of the

expansion for v(x) in terms of the basis functions φj (x).

C

C We may convert these coefficients into those of the ordinary

C Chebyshev series as described in Sec. 5, but this is optional.

6.15. SUMMARY 125

Make a graph of u(x) at 101 points by summing the basis function

C

C [not Chebyshev] series.

UGRAPH(1) = ALPHA

UGRAPH(101) = BETA

DO 300 I = 2,100

X = -1 + (I-1)*0.02

CALL BASIS(X,NBASIS)

Recall that u(x) = v(x) + B(x) where v(x) is given by the

C

computed series and B(x) is the R. H. S. of the next line.

C

UGRAPH(I) = (ALPHA*(1-X) + BETA*(1+X))/2

DO 350 J=1,NBASIS

350 UGRAPH(I)=UGRAPH(I) + APHI(J) * PHI(J)

[At the end of this loop, UGRAPH contains

300 CONTINUE

the desired values of u(x) . . . input them to

your favorite plotter!]

END

SUBROUTINE BASIS(X,NBASIS)

COMMON/BASIS/PHI(20),PHIX(20),PHIXX(20)

C After this call, the arrays PHI, PHIX, and PHIXX contain the

C values of the basis functions (and their first two derivatives,

C respectively, at X.

IF (ABS(X).LT.1) THEN [block below is executed only on the

interior. This IF branches to a later

block to evaluate functions at x = ±1.]

[T is the argument of the trig. functions]

T = ACOS(X)

C = COS(T)

S = SIN(T)

DO 100 I=1,NBASIS

N = I+1

[Trig. form of TN (x)]

TN = COS(N*T)

TNT = - N * SIN(N*T) [Derivative of TN with respect to t]

[Second t-derivative of TN ]

TNTT= - N*N * TN

Convert t-derivatives into x-derivatives

C

[x-derivative of N-th Chebyshev polynomial.]

TNX = - TNT / S

TNXX= TNTT/(S*S) - TNT*C/(S*S*S) [Second x-derivative]

Final step: convert TN s into the basis functions φN s.

C

We subtract 1 from the even degree T N s, x from the odd degree

C

TN and 1 from the first derivative of the odd polynomials only

C