. 90
( 136 .)


18.30 Spherical Harmonics Besieged
Although spherical harmonics are the discretization of choice for the ECMWF forecasting
model, the NCAR Community Climate model, and Glatzmaier™s stellar convection and
mantle ¬‚ow simulations at the turn of the milllenium, this spectral basis is under attack on
two fronts. The ¬rst assault is from massively parallel (MP) computing. The worry is that
global basis sets of any kind, and spherical harmonics with their slow Legendre transforms
in particular, will become lamentably inef¬cient when the number of processors is huge.
As a result, a number of national weather services have switched to ¬nite element or other
local, low order algorithms. Some American research groups, such as those at NASA,
employ ¬nite difference codes almost exclusively.

Table 18.12: Comparisons of Finite Difference with Spectral Models: Bibliography

References Comments
Doron&Hollingsworth Weather forecasting
Hoskins&Simmons(1975) Comparison of spherical harmonic/differences in z
with 3D ¬nite difference weather forecasting code
Simmons&Hoskins(1975) Spectral vs. ¬nite-differences: growing baroclinic wave
Merilees et al.(1977) FD with polar ¬ltering compared with pseudospectral
Browning&Hack ¬nite differences of 2d, 4th, and 6th order
&Swarztrauber(1988) compared with spherical harmonics
Held&Suarez (1994) 2d/4th ¬nite difference; climate model
Gustafsson&McDonald(1996) HIRLAM grid point vs. spectral semi-Lagrangian;
models are similar in accuracy and ef¬ciency
Fornberg&Merrill(1997) Double Fourier vs. spherical harmonics vs. 2d & 4th
order differences for a very smooth solution
Spotz&Taylor Comparisons of compact 4th order differences with
& Swarztrauber(1998) spherical harmonics and ¬ltered double Fourier

So far, however, the MP/low order assault has not yet carried the fort although some
of the outworks have fallen. In the ¬rst place, ingenious coding and the inherent ef¬ciency
of matrix multiplication has made spectral codes run fast even on machines with moderate
parallelism (forty-eight processors). (See, for example, Drake et al.(1995) and other articles
in a special issue of Parallel Computing, vol. 21, number 10.) Because the architecture of
twenty-¬rst century machines is still in a state of ¬‚ux, it is impossible to say whether the

spherical harmonics are doomed or whether ingenuity will continue to make them com-
petitive. Second, the new ¬nite order methods on the sphere have a shorter history than
spherical harmonics, and consequently many schemes are still nagged by stability prob-
lems and also by their low accuracy for smooth, large scale structures that are much better
resolved by spectral codes.
Most of the studies in Table 18.12 which have compared ¬nite difference and spectral
codes have found the latter were better. However, Held and Suarez(1994) have claimed
that a low order difference model is just as accurate and runs faster; Gustafsson and Mc-
Donald(1996) report a tie. The meteorological community is applying more sophisticated
benchmark tests including both smooth and discontinuous solutions, so these controver-
sies will eventually be resolved. In the meantime, spherical harmonics models soldier on.
The second assault is the replacement of spherical harmonics by an alternative high
order method. Spectral elements, which allow the work to be split among many processors
without sacri¬cing high accuracy, have great promise (Table 18.11. However, only two
groups, one each in meteorology and oceanography, are developing such models. Two-
dimensional Fourier series are very promising, too (Spotz, Taylor and Swarztrauber, 1998,
and Fornberg and Merrill, 1997, Shen, 1999, and Cheong, 2000). However, work to date
has been con¬ned to the shallow water wave equations ” two-dimensional, no clouds, no
chemistry. It seems likely that one or both of these alternatives will someday be used for
operational forecasting, but when is very murky.
The third assault is from conservative, shock-resolving van Leer/Roe schemes, which
have become very popular in aerodynamics. The Lin and Rood algorithm (1996) has gained
a wide following. However, experience with such low-order-but-conservative schemes is
still limited.
The only conclusion we can make at the turn of the twenty-¬rst century is: the siege

18.31 Elliptic and Elliptic Cylinder Coordinates
Two-dimensional elliptic coordinates (µ, ·) are de¬ned by

y = ’ sinh(µ) sin(·) (18.94)
x = cosh(µ) cos(·) &

where x and y are the usual Cartesian coordinates in the plane. The reason for the name
“elliptic” is that the surfaces of constant µ are all ellipses with foci at x = ±1 along the
x-axis as illustrated in Fig. 2-16 of Chapter 2. At large distances from these foci, the ellipses
degenerate into circles and elliptical coordinates become ordinary polar coordinates. In
this limit, µ ≈ log(r) where r is the polar radius and · ≈ θ, the polar angle.
Elliptic cylinder coordinates add a vertical coordinate z which is perpendicular to the
x-y plane. These three-dimensional coordinates generalize elliptic coordinates in the same
way that cylindrical coordinates extend polar coordinates to three dimensions. Since the
z-axis is straight, one should use either Chebyshev polynomials or rational Chebyshev
functions for z, just as for cylindrical coordinates.
Because the quasi-angular coordinate · is cyclic, all functions are periodic in ·, which
implies a Fourier basis in ·. However, the coordinate singularity implies that µ-dependent
Fourier coef¬cients have de¬nite parity, just as for spherical and polar coordinates. How-
ever, there are intriguing differences because the singular coordinate µ = 0 is a singular
line (the whole interval x ∈ [’1, 1]) rather than a single point.


If F (µ, ·) is an analytic scalar function which is expanded as

F (µ, ·) = fm (µ) cos(m·) + gm (µ) sin(m·)

then continuity of F (µ, ·) and all its derivatives at µ = 0 implies:

(i) fm (µ) is symmetric with respect to µ = 0 for all m

(ii) gm (µ) is antisymmetric for all m

Thus, one should expand the coef¬cients of cos(m·) as even Chebyshev polynomials in µ while
the sine terms should be represented as a sum of odd polynomials only.

PROOF: All functions in two dimensions can be approximated to arbitrary accuracy as
polynomials in the Cartesian coordinates x and y. This in turn implies that any well-
behaved function can be written as a sum of the powers

xj y k = cosj (·) sink (·) coshj (µ) sinhk (µ) (18.96)

The parity of xj y k with respect to both · and µ is controlled entirely by k. If k is even,
then sink (·) and sinhk (µ) are both symmetric functions; if k is odd, then (18.96) is also odd
in both elliptic coordinates. Now a function which is symmetric in · may be expanded as
a cosine series in ·. The fact that even parity in · ←’ even parity in µ implies that all
the cosine coef¬cients must be symmetric in µ, which is proposition (i). Similar reasoning
implies (ii).
One could alternatively prove the theorem by showing that it is necessary for continuity
of F (µ, ·) and all its derivative across the coordinate singularity at µ = 0. (By singularity,
we mean that differential equations in elliptical coordinates have coef¬cients that are sin-
gular at µ = 0; the solution is analytic on that line segment.)
Boyd and Ma(1990) and Mittal and Balachandar(1996) are ¬‚uid mechanics applications
of spectral methods in elliptic coordinates, Fourier in · and Chebyshev polynomials in µ.

18.32 Summary
The end-of-the-century state of the art is a mixture of the known and unknown. The known
includes: coping with coordinate singularities and implementing stable, fast algorithms in
spherical, cylindrical or polar coordinates and in some related coordinate systems; these
methods run well on coarse-grained parallel computers with sixteen or fewer processors.
One unknown is: Can the Legendre transform be accelerated to keep spherical harmon-
ics or One-Sided Jacobi polynomials competitive at very large N ? Another unknown is:
Can spherical harmonics schemes and related algorithms be ef¬ciently adapted to massively-
parallel machines? Will spherical harmonics be replaced by latitude-and-longitude Fourier
series or spectral elements? Will spectral methods of all kinds be routed by low order ¬nite
difference or ¬nite element methods for weather forecasting and climate modelling?
In the areas of this chapter, the Age of Exploration is still not over.
18.32. SUMMARY 441

Table 18.13: Reviews, Books, and Model Descriptions

References Comments
Bourke&McAvaney Comprehensive review, describing spherical harmonics code which is the
&Puri&Thurling(1977) basis for ECMWF forecast and NCAR climate models
Haltiner&Williams(1980) one chapter on spherical harmonic methods
Sela (1980, 1995) Reviews of National Meteorological Center (USA) forecast model
Daley(1981) Review concentrating on normal mode initialization
for weather prediction models on sphere
Gordon&Stern(1982) Description of GFDL global circulation model
Kanamitsu et al.(1983) Description of the operational model
of the Japanese Meteorological Agency
Boer et al.(1984) Review of Canadian Climate Centre model
Jarraud &Baede (1985) lengthy review of atmospheric codes
Daley(1991) Monograph about data analysis, initialization and
other issues for weather prediction
Hogan&Rosmond(1991) U. S. Navy global operational spectral code
Jakob-Chien & seven benchmark problems for shallow water eqs.;
&Hack&Williamson(1995) fronts, waves, Gibbs™ wiggles
Chapter 19

Special Tricks

“If it works once, it™s a trick. If it works twice, it™s a method. If it works three times, it™s a
” Source Unknown

“Card tricks are as important in numerical analysis as in stage magic.”
” J. P. Boyd

19.1 Introduction
In this chapter, we give brief, self-contained descriptions of a number of ideas that are use-
ful even if they are applicable only under special circumstances. These special tricks can be
safely skipped by the reader who is interested only in the major themes of spectral meth-
ods. Nonetheless, the assertion that “numerical analysis is as much an art as a science” has
become a proverb. Much of the reason is precisely because special technques are the best
way of solving certain exotic problems; card tricks, if we may so dub the ideas illustrated
in this chapter, are as important in computation as in conjuration.
Besides the special algorithms collected in this chapter, we have discussed a number of
useful tricks elsewhere including the following.
First, one can use nonlinear degrees-of-freedom, such as width or shape parameters,
in solving differential equations. Recall that when the numerical solution is approximated
by a series, uN (x) is a linear function of the unknown series coef¬cients, a0 , . . . , aN . But
it need not always be so. Nonlinear unknowns were used with great success in the H2
ion (Chapt. 3, Sec. 7) and the Korteweg-deVries soliton (Appendix G of the ¬rst edition
of this book). (An especially good candidate as a nonlinear degree-of-freedom is a width
parameter that matches the scale of the basis functions to the width of the solution.)
Second, the standard methods of Chapt. 17 for solving problems on an in¬nite interval
will fail unless the oscillations in u(x) decay exponentially as |x| ’ ∞. (This condition
may be met either because u(x) itself decays exponentially, with or without oscillations,
or because the oscillations rapidly decay even though u(x) asymptotes to a constant or
its mean value decreases algebraically with |x|.) If u(x) oscillates with slowly decaying
amplitude as |x| ’ ∞, then all standard methods fail because the rational Chebyshev


functions (or whatever) are being asked to faithfully mimic an in¬nite number of crests
and troughs. Chapt. 17, Sec. 13, shows that one may solve this problem, at least for the
Bessel functions Jn (x), by computing two simultaneous series: one for the amplitude and
one for the phase of the oscillation.
The Bessel double expansion is a particular case of a very general strategy: modifying
a standard basis set by using special functions that match the pathological behavior of the
solution. Other examples are described in Secs. 5, 6, and 7 below.
Third, when the solution is singular or nearly singular on the interior of an interval, one
can detour around the branch point by integrating on an arc in the complex plane. This de-
vice, very useful for computing weakly unstable eigenvalues in hydrodynamic instability
theory, is described in Chapt. 7, Sec. 11.
Fourth, when u(x) has regions of steep gradient, one can use a change-of-coordinate to
cluster grid points in the frontal region. Chapt. 13 describes a wide variety of real-valued
mappings; transformations for spatially periodic problems are in Sec. 7 and adaptive grids
are explained in Sec. 9.
These ideas and those explained below are only a sampler. No doubt the reader will
invent his own variations.

19.2 Sideband Truncation
Normally, we truncate a spectral expansion by including all basis functions up to a certain
cutoff, beginning with φ0 (x). For special problems, however, this may be inef¬cient. The
large eigenfunctions of Mathieu™s equation furnish a good example. These solve

uxx + [» ’ 2 q cos(2x)] u = 0 (19.1)

subject to the condition of periodicity with period 2π with » as the eigenvalue and q as a
constant parameter.
When q = 0, the eigenvalues and eigenfunctions are

» n = n2 un (x) = sin(nx) or (19.2)
; u0 = 1; cos(nx)

For non-zero q, the unperturbed solutions are still good approximations to those of the full
Mathieu problem, (19.1), whenever n2 q. If n = 15, for example, the spectral coef¬cient
with the largest amplitude will not be a0 or a1 , but a15 . If we use the most drastic possible
truncation ” one term! ” then sin(15 x) or cos(15 x) must be that term.
One can show by using the trigonometric identity

cos(2x) cos(nx) = 0.5 {cos[(n + 2)x] + cos[(n ’ 2)x]} (19.3)

(and its analogue for sin(nx)) that the perturbation, 2 q cos(2x), will only couple Fourier
terms whose degrees differ by two and are of the same type, either sines or cosines. The
eigenfunctions of Mathieu™s equation therefore fall into four separate classes. The eigen-
functions of each class can be computed independently of the eigenfunctions of the other
three classes.
These four groups of eigenfunctions illustrate the double-parity symmetry of Fourier
functions (Chapter 8, Sec. 2). Each eigenfunction has de¬nite parity both with respect to
x = 0 and with respect to x = π/2. The little table below lists the four classes in the
conventional notation. The middle column gives the double-parity symmetry of each class
with the parity with respect to x = 0 listed ¬rst. The third column lists the Fourier functions
that suf¬ce to represent all modes in a given class; each of these basis functions has the same
double symmetry as the Mathieu functions.


. 90
( 136 .)