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(1974)

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

18.31. ELLIPTIC AND ELLIPTIC CYLINDER COORDINATES 439

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

continues.

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.

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

440

Theorem 37 (ELLIPTICAL COORDINATES: PARITY in QUASI-RADIAL COORD.)

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

∞

(18.95)

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

m=0

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

law.”

” 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

442

19.2. SIDEBAND TRUNCATION 443

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.

CHAPTER 19. SPECIAL TRICKS

444