. 78
( 136 .)


Fourier or Chebyshev coefficients
0 20 40 60 80 100

Figure 17.13: Absolute values of spectral coef¬cients for the rational Chebyshev series
(solid) and the Cloot-Weideman algorithm (disks) for the function f (y) = sech(y). The
dashed line is a ¬t to the leading coef¬cients of the T B series.

jump discontinuity at t = π whose magnitude is the domain truncation error. The ratio-
nal Chebyshev series ¬‚attens out, as shown by comparison with the dashed dividing line,
which was chosen to match the slope of the coef¬cients aj for j < 40. However, the error
continues to decrease exponentially fast, albeit at a subgeometric rate. The curves provide
some useful empirical guidance for optimizing L as expressed by the following:

Plot the coef¬cients aj versus degree on a log-linear plot. If the graph abruptly ¬‚attens at some N ,
then this implies that L is TOO SMALL for the given N , and one should increase L until the ¬‚atten-
ing is postponed to j = N .

Fig. 17.14, which is identical to Fig. 17.12 except for a different u(y) and N , illustrates
how dif¬cult it is make rigid rules about the superiority of one method over another. For
the function u(y) = sech(y) cos(3y), the smallest errors attainable for a given N are roughly
the same for each choice of basis. One would probably be more comfortable using the
rational Chebyshev basis because it is less sensitive to the choice of the map parameter.
For solving model time evolution equations, Weideman and Cloot (1990) concluded
that the sinh-map gave the smallest errors and was most ef¬cient when the width of the
solution changes little with time. However, the rational Chebyshev basis was superior
when the width of the solution varied greatly with time. Both authors have continued to
employ both schemes in their own later work.
Their own inability to choose one over the other reiterates a key theme: on the in¬nite
and semi-in¬nite intervals, there are a lot of good choices, but no perfect choice.

17.15 Summary
When the target solution u(y) decays exponentially as the coordinate tends to in¬nity, there
are many good spectrally-accurate algorithms. Sinc functions are the simplest, Hermite
and Laguerre functions have the closest connection with the physics for some classes of
problems, rational Chebyshev functions and the Weideman-Cloot sinh-mapped Fourier

-3 -3
10 10

-4 -4
10 10

-5 -5
10 10

-6 -6
10 10

-7 -7
10 10
0 5 10 15 0 1 2
L L/Lopt (WC)

Figure 17.14: Errors in the rational Chebyshev series (left panel) and the Weideman-Cloot
algorithm for the function u(y) = sech(y) cos(3y) for various map parameters L or domain
parameters L/Lopt where Lopt ≡ (1/π) log(N ) is the Cloot and Weideman™s own estimate
of the best L. Twenty-two rational Chebyshev functions were used for the left panel; the
same number of cosine functions for the right graph.

scheme exploit the Fast Fourier Transform. We have not made a general recommenda-
tion about the “best” basis because there is no compelling reason to choose one basis over
another that applies to all cases.
When the function u(y) decays as an inverse power of y, then rational Chebyshev func-
tions are the clear favorite. This basis is equivalent to making a change of coordinate, such
as y = Lcot(t) for the in¬nite interval, and then applying a Fourier basis in the new coordi-
nate t. However, one must analyze the asymptotics of the solution carefully to ensure that
the function u(y[t]) is well-behaved in t. If u(y) decays slowly but as some fractional power
of y, then one can still obtain spectral accuracy. However, one must modify the cotangent
mapping so that u(y) is in¬nitely differentiable everywhere in the new coordinate t.
Another class of nasty problems is when a function de¬ned on a semi-in¬nite interval
is singular at one endpoint or on or near the interior of the interval. The singularities can
be resolved by combining the ideas of this chapter with those of the previous chapter. For
an endpoint singularity, for example, one can map the semi-in¬nite interval to t ∈ [0, π]
using a change-of-coordinate that clusters the grid points with a spacing that decreases
exponentially near the singular endpoint y = 0 as illustrated by the series for K1 Bessel
function above. Similarly, one can choose the map to detour around a singularity on the
interior of the interval whether the interval is ¬nite or unbounded.
When u(y) decays to an oscillation rather than to zero as |y| ’ ∞, the best strategy is
to use special basis functions, added to a standard in¬nite interval basis, which match the
oscillations, leaving the standard basis to approximate the part of the function which de-
cays without oscillation. We have illustrated this strategy using the amplitude-and-phase
approximations to the J0 Bessel function. However, the same idea works well in quantum
scattering and computations of weakly nonlocal solitary waves.
For all basis sets on an unbounded interval, there is a scale factor or map parameter that
(if chosen properly) matches the width of the basis functions to the width of the solution.
One must use much narrower basis functions to compute the structure of a molecule than
to compute a solitary wave approximation to the Great Red Spot of Jupiter, which is 40,000
17.15. SUMMARY 379

km across.
For Whittaker cardinal functions, the scaling factor is the grid spacing. When Cheby-
shev polynomials are mapped from [’1, 1] to [’∞, ∞] to give the rational basis functions
denoted by T Bn (y), the scaling factor is the map parameter L. The Hermite functions do
not explicitly contain a free parameter, but we always have the option of performing the
expansion in terms of ψn (±y) instead of ψn (y) where ± is an arbitrary scaling factor. The
proper choice of ± is just as important to maximizing the ef¬ciency of the Hermite series
as is a good choice of L for mapped Chebyshev polynomials. When the domain is un-
bounded, there is always, either explicitly or implicitly, a scale parameter to be chosen.
Chapter 18

Spherical, Cylindrical, Toroidal &
Elliptical Geometry

“The book of nature is written in the characters of geometry”
” Galileo Galilei

18.1 Introduction
Cylindrical and spherical geometry require special theory, special grids, and special basis
functions. There is no easy road to compute ¬‚ows on the surface of a sphere or solve a
partial differential equation on the interior of a disk.
The heart of the problem is that at the origin of a polar coordinate system, or at the north
and south poles of the surface of a sphere, the lines of constant polar angle or constant
longitude converge in a single point. This convergence-of-meridians has two unfortunate
consequences. First, there is a severe time-stepping limit, the so-called “pole problem”.
Second, the differential equation is almost invariably singular at this convergence point
even when the solution is everywhere smooth and in¬nitely differentiable.
In this chapter, we hope to initiate the reader into the mysteries of disk and sphere. In
the ¬rst part, we discuss polar coordinates. In the second, we concentrate on latitude and
longitude, the coordinates of the surface of a sphere. As explained in the next section, polar,
cylindrical, spherical, and toroidal coordinates are very closely related: there is really only
one species of “pole problem”, manifesting itself as a disease in four different coordinate
Spectral algorithms for all these coordinates have been developed. The numerical tech-
nology has become so mature that spectral methods for weather forecasting and climate
modelling have conquered a global empire embracing more countries than that of Caesar
or Alexander. Nevertheless, as the century ends and massively-parallel computing gathers
momentum, spherical harmonic models are under siege. As we describe the techniques,
we shall narrate both the virtues and faults of each spectral and non-spectral option and
let the reader be his own prophet.





Figure 18.1: Left: De¬nitions of radius r and angle θ for a typical point (black disk) in polar
coordinates. Right: In polar coordinates, the contours of constant radius r are concentric
circles. The contours of constant angle θ are semi-in¬nite rays through the origin (thick

18.2 Polar Coordinates and Their Relationship to Cylindri-
cal, Spherical and Toroidal Coordinates
In polar coordinates, a point is speci¬ed by the radius r, which is the distance from the
origin to point, and the angle θ, as illustrated in the left side of Fig. 18.1. The Cartesian
coordinates x and y are connected to polar coordinates through the relation

x = r cos(θ), y = r sin(θ)

The contours of constant radius r are concentric circles centered on the origin. The isolines
of polar angle θ are semi-in¬nite rays all converging at the origin as shown on the right of
Fig. 18.1.
Cylindrical coordinates are obtained by adding a third coordinate z which measures
distance along a vertical axis perpendicular to the plane de¬ned by polar coordinates. The
isosurfaces of constant r are cylinders. The rays of constant θ now converge at every point
on the vertical axis. Each plane of constant z is a two-dimensional surface where points are
located by polar coordinates.
Toroidal coordinates are obtained by wrapping the cylinder around into the shape of
a bagel or doughnut. The surfaces of constant r are now tori, and the coordinate z, now
renamed the “toroidal angle” » (Fig. 18.3), is always cyclic.
In spherical coordinates, the radial coordinate is now the distance from the origin to a
point in three dimensions rather than two, and the surfaces of constant radius are spheres.
Latitude and longitude, as on any geographical map, locate points on the surface of a
given sphere. (Actually, in scienti¬c applications, it is customary to replace latitude by
“colatitude”, which is latitude minus π/2, but this is merely a convention that does not
alter the pole problem.)
A “polar cap” is a small region centered on either the north or the south pole, bounded
by colatitude θ = θ0 where θ0 is small. The polar cap is a slightly puckered disk. As
θ0 ’ 0, that is, as the polar cap becomes smaller and smaller, the polar cap becomes ¬‚atter
and ¬‚atter. In the limit, the convergence of meridians (lines of constant longitude) at the

Table 18.1: Spectral Methods in a Disk or Cylinder
Note: Exterior ¬‚ows and solutions in an unbounded domain are listed in the next table.
References Comments
Mercier & Raugel (1982) Theoretical justi¬cation for pole conditions
in mixed spectral/¬nite element [in French]
Orszag&Patera (1983) Robert-type basis in polar coordinates
Tan(1985) 3D Poisson equation in cylinder
Randriamampianina& ¬‚ows with buoyancy and rotation
Bontoux&Roux(1987) in a cylinder (in French)
Lewis&Bellan(1990) Symmetry conditions for scalars & vectors
across r = 0 in cylindrical coordinates
Eisen&Heinrichs mapped unit disk onto rectangle
&Witsch(1991) and imposed pole conditions
Bouaoudia&Marcus (1991) Robert basis (now obsolete)
Pulicani&Ouazzani(1991) 3D ¬‚ow and diffusion in cylinder
Huang &Sloan(1993b) singularity analysis; good discussion of
¬nite difference preconditioning
Bessel eigenproblem & Poisson eq. in a disk
Launaye et al.(1994) axisymmetric 2D crystal growth in a cylinder
Van Kemenade&Deville(1994b) spectral elements; non-Newtonian ¬‚ow
Godon&Shaviv(1993,1995) Chebyshev in r; Fourier in θ
Godon(1995,1996abc,1997a,b) astrophysics problems
Godonet al.(1995) accretion disk boundary layers
Priymak(1995),Priymak&Miyazaki(1998) Turbulence
Fornberg(1995) Chebyshev/Fourier basis with strong polar ¬ltering
Matsushima and Marcus(1995) One-Sided Jacobi basis
Raspo& Ouazzani & Multidomain scheme for axisymmetric ¬‚ow within
Peyret(1994,1996) a cylinder
Shen(1997) fast algorithms for Poisson eq. and related
operators in disk or cylinder
Verkley(1997a,1997b) One-sided Jacobi basis; hydrodynamics
Matsushima&Marcus(1997) unbounded vortical ¬‚ows, polar coordinates
Lopez&Shen(1998) Navier-Stokes equations with fast algorithms
for the semi-implicit time marching

north pole is indistinguishable from the convergence of the rays of constant θ at the origin
in polar coordinates.
Thus, the special dif¬culties caused by the convergence of isolines of an angular coor-
dinate at a single point are identical in all four coordinate systems.

18.3 Polar Coordinates: Apparent Singularity at the Origin


. 78
( 136 .)