<<

. 28
( 136 .)



>>

However, articles such as Boyd (1986c) and Schultz, Lee and Boyd (1989) have already
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

<<

. 28
( 136 .)



>>