. 80
( 136 .)


φm (r, θ) ≡ Wn (r)
n sin(mθ)

where in Verkley™s notation
Wn (r) ≡ rm P(n’m)/2 (2r2 ’ 1),
m = 0, 1, 2, . . . ; n = m, m + 2, m + 4, . . .
where Pk (s) is the Jacobi polynomial of degree k of order (0, m) in its argument s. (Note
that because s ≡ 2r2 ’ 1, the Jacobi part of the basis function is of degree (n ’ m) in r.)
This basis function explicitly enforces all the constraints of the parity theorem; each basis
function is symmetric in r for even m and odd in r for odd m. The factor of rm ensures that
each basis function has an m-th order zero at the origin. Verkley(1997a, Table 1) provides
the explicit form of these basis functions for small n and m.
The Jacobi polynomials are orthogonal so that
±,β ±,β
Pk (s)Pk (s)(1 ’ s)± (1 + s)β ds = 0, (18.17)

Changing variables from s to r and setting ± = 0, β = m gives
0,m 0,m
Pk (2r2 ’ 1)Pk (2r2 ’ 1)r2m rdr = 0, (18.18)

or in other words,
m m
Wn (r)Wn (r) rdr = 0, n=n

These Jacobi polynomials of order (0, m) are “one-sided” in the sense that the weight
function in Eq.(18.18) diminishes rapidly (as r2m ) as the origin is approached, but varies
only slightly near the outer edge of the disk. If the Jacobi polynomials oscillated uniformly

in r like the Robert functions, they would be nearly-dependent and ill-conditioned like the
Robert functions. The orthogonality constraint, however, forces the polynomials of differ-
ent degree to be as independent as possible. This requires that the polynomials oscillate
mostly near r = 1: the roots of the basis functions move closer and closer to the outer
boundary for ¬xed degree as the order m increases.
The one-sideness is the magic. The basis functions of high wavenumber m are the
villians that describe the fast tangential advection or diffusion that limits the timestep.
However, the radial parts of these basis functions, if one-sided Jacobi polynomials, have
only negligible amplitude close to the origin where the grid points are close. The rapid
advection or diffusion is suppressed, and a much longer timestep is possible.
As with spherical harmonics, discussed later, there are two options for truncating the
basis: “rectangular” and “triangular”. A rectangular truncation employs the same num-
ber of radial basis functions for each angular wavenumber m. It is the easiest to program,
but has fallen from grace in spherical coordinates. The triangular truncation decreases the
number of radial basis functions by m with each increase in m until the highest wavenum-
ber has but a single radial basis function. It has been shown for spherical harmonics that
this gives the most uniform resolution (and allows the longest timestep) of any truncation
that includes a given maximum zonal wavenumber M . This property of “equiareal resolu-
tion”, de¬ned formally below in Sec. 18.13, has not been proved for the One-Sided Jacobi
basis, but it seems plausible that a triangular truncation is preferable here, too.
The usual strategy of choosing the roots of the polynomials as the interpolation points
are different for each polar wavenumber
does not work here because the zeros of the Pk
m. To take transforms in θ, we need the unknowns de¬ned on a single set of radial grid
points. Therefore, one is forced to abandon a pseudospectral strategy and instead apply
a Galerkin method. The integration points are the Legendre-Radau grid which includes
r = 1 but not r = 0. To avoid quadrature errors, it is necessary to use more grid points than
unknowns, but the transform from radial grid point values to coef¬cients can be still be
expressed as the multiplication of a vector of grid point values or coef¬cients by a matrix,
albeit now a rectangular matrix.
These basis functions, unlike the spherical harmonics they mimic, are not eigenfunc-
tions of Laplace operator. However, both Matsushima and Marcus(1995) and Verkley(1997a)
derive recurrence relations which allow the Laplacian to be inverted by solving a penta-
diagonal matrix problem. Similarly, derivatives can be evaluated at an operation count
directly proportional to the number of spectral coef¬cients in the truncation. Overall, the
ef¬ciency of this basis is roughly the same as for spherical harmonics. Indeed, the algorith-
mic connection with spherical harmonics is so close that Verkley™s shallow-water solver
for the disk is merely a rewrite of a previous code that solved the same hydrodynamic
equations on the surface of a sphere via spherical harmonics.
Unfortunately, the Jacobi basis shares the vices of spherical harmonics. First, there is no
Fast Fourier Transform for the Jacobi polynomials, so all radial transforms from grid point
values to spectral coef¬cients must be done by Matrix Multiplication Transform (MMT).
This costs O(N 2 ) for each wavenumber m versus only O(N log2 (N )) for an FFT. (Trans-
forms in θ are still done by FFT, thank goodness.)
Second, the messy recurrence relations and basis functions require some investment of
both learning time and programming time. (We may hope that someday the necessary
subroutines will become widely available in software libraries, and then this objection will
There is one profound difference from spherical harmonics, however. The One-sided
Jacobi basis has the usual high density of grid points near r = 1 where grid points are sepa-
rated in r by O(1/N 2 ). With a triangular truncation so that the maximum polar wavenum-
ber M is equal to N , the tangential distance between grid points on the outermost ring is

Chebyshev Expansions of J0


T (2 r - 1)

Tj(2 r - 1),

T (r)
-10 2j

0 2 4 6 8 10
degree j

Figure 18.2: Chebyshev expansions of the Bessel function J0 (r). The slowly-convergent
series (dashed) is the expansion in Shifted-Chebyshev polynomials: Tj— (r) ≡ Tj (2r ’ 1).
The two good techniques are (i) an expansion in Shifted-Chebyshev polynomials with a
QUADRATIC argument Tj (r 2 ) = Tj (2r 2 ’ 1) and (ii) an expansion in the EVEN Chebyshev
polynomials without a shift: T2j (r). These two options generate only a single curve be-
cause (i) and (ii) are in fact IDENTICAL because of the identity Tj (2r2 ’ 1) = T2j (r).

also O(1/N 2 ). This implies that the One-sided Jacobi basis has not completely solved the
dif¬culty of a restrictive timestep; indeed, if the time scales for tangential advection or dif-
fusion near the origin are similar to the time scales for radial advection or diffusion near
the boundary, then the switch from a Chebyshev basis to the Jacobi basis will not increase
the time step at all.
Nevertheless, this basis has worked well in Matsushima and Marcus(1995) and Verkley

18.5.2 Boundary Value & Eigenvalue Problems on a Disk

When there is no timestep, the pressure to use One-Sided Jacobi basis functions in radius is
greatly reduced. Orszag(1974) and Boyd(1978c) have shown that satisfying all the pole con-
ditions is not terribly important to numerical ef¬ciency; a regular Chebyshev basis seems
to work just ¬ne. (Their papers actually solve problems on the sphere, but because the
polar caps are well-approximated by small disks, their conclusions should apply to polar
coordinates, too.) Matsushima and Marcus(1995) show by solving the eigenproblem whose
exact solution is the Jm Bessel function that the Jacobi basis can reach a given error toler-
ance with only half as many coef¬cients as the Parity-Restricted Chebyshev polynomials
when m = 50. However, in a Fourier-Chebyshev series in two dimensions, the Fourier
series in m usually converges rapidly, so the Chebyshev basis is converging slowly only
for wavenumbers whose amplitude is negligible anyway.
Our recommendation for boundary and eigenvalue problems is to try a Chebyshev
basis ¬rst. If the speed of the resulting (rather easy-to-program) code is unsatisfactory,
the One-Sided Jacobi basis can be used as a fallback with the Chebyshev program used to
check the newer and much more intricate Jacobi code.

Table 18.2: Cylindrical or Polar Coordinates: Unbounded Domain or External to a Cylinder

References Comments
Zebib(1987a) Stability of ¬‚ow past a cylinder
Don&Gottlieb(1990) unsteady ¬‚ow past a cylinder
Deane & Kevrekidis spectral elements; low order basis of
&Karniadakis & Orszag(1991) empirical orthogonal eigenfunctions
Mayer&Powell(1992) Eigenvalue calculation of stability of trailing
vortex through domain truncation
Mittal & Balachandar(1996) ¬‚ow past elliptic cylinder in elliptic coordinates
good discussion of boundary conditions,
blending in¬‚ow & out¬‚ow
Blackburn&Henderson(1996) ¬‚ow past a vibrating cyclinder
Matsushima&Marcus(1997) Unbounded domain including the origin;
2D and 3D vortex ¬‚ows

18.5.3 Unbounded Domains Including the Origin in Cylindrical Coor-

Matsushima&Marcus(1997) have extended their earlier work by applying an algebraic
change-of-coordinate in r to map an unbounded domain into a disk. They then apply the
One-Sided Jacobi basis described above. Like the rational basis T Bj for the in¬nite interval
described in Chapter 17, the Matsushima-Marcus basis functions are rational functions of
r ∈ [0, ∞].
The alternative is the even T B functions for even polar wavenumber m and the T B2j’1 (r)
for odd m: a Parity-Restricted Rational Chebyshev basis. The simpler functions are prefer-
able for boundary and eigenvalue problems, but the Matsushima-Marcus basis may give
faster convergence and allow a longer timestep.
Flows exterior to a cylinder in an unbounded domain have been the subject of much
study, both analytical and numerical. Exterior ¬‚ows are free of the pole problem, but im-
posing boundary conditions on a ¬‚ow at large radius can be tricky. We have therefore
collected some illustrations in Table 18.2.

18.6 Annular Domains
Flow in an annulus that does not include the origin is free from the pole problem. We
recommend a Fourier basis in the angle θ and a standard Chebyshev basis in the radial
coordinate r. In three dimensions, one should use a Chebyshev basis in the axial coordinate
z unless periodic boundary conditions in z are imposed, in which case a Fourier basis is
greatly preferable.
There is little difference between annular ¬‚ow and channel ¬‚ow except that cylindri-
cal coordinates introduce an r-dependent “metric” factor into the Laplace operator, etc.,
whereas the Laplace operator has constant coef¬cients in Cartesian coordinates. The usual
methods for separable partial differential equations, discussed in Chapter 15, Sec. 11, still
For most practical purposes, however, an annular channel is simply a straight channel
which is periodic in the downstream direction. A standard Fourier/Chebyshev basis is
best. Some representative works are listed in Table 18.3.

Table 18.3: Annular Flows
References Comments
Marcus(1984a,1990) Taylor-Couette ¬‚ow
Le Qu´ r´ &Pecheux(1989,1990)
ee axisymmetric convection; bifurcations
Chaouche(1990a) axisymmetric ¬‚ows;
Chaouche et al.(1990) in¬‚uence matrix method
Randriamampianina(1994) 3D ¬‚ow; vorticity-vector-potential

18.7 Spherical Coordinates: An Overview
The pattern of latitude and longitude lines on the surface of the sphere is, near the poles,
the same as that of a polar coordinate system in a disk. So, it is hardly surprising that there
are many similarities between cylindrical coordinates and spherical coordinates. However,
there are important differences.
One is that a sphere has two poles rather than one. A second, more important difference
is that the surface of a sphere has no boundaries. This has the profound consequence that
spherical basis sets automatically and individually satisfy the behavioral boundary condi-
tions on the sphere. In contrast, it is necessary to impose numerical boundary conditions at
the boundary circle of the disk.
A third difference is that for the sphere, there is an “obvious” basis set, the spherical har-
monics, which are nearly ideal. Spherical harmonics give equiareal resolution, exponential
convergence, and trivial inversion of the Laplace operator, which is the eigenoperator for
these functions. In contrast, the obvious basis for the disk, the cylindrical harmonics, are
Bessel functions of radius with only a poor, algebraic rate of convergence. One can match
most of the good properties of spherical harmonics by using the One-Sided Jacobi basis for
the disk, but these are not eigenfunctions of the Laplace operator.
The fourth difference is that the number of simulations in a cylinder is fairly modest
compared to the large body of work that has been done in spherical coordinates for weather
forecasting, climate modelling, mantle convection and stellar ¬‚ows. For this reason, we
shall discuss spectral methods on the sphere in great detail.

18.8 The Parity Factor for Scalars: Sphere versus Torus
Topologically, the sphere is a two-dimensional manifold of genus zero while the torus
(“doughnut”) is of genus one. Ironically, however, it is easier to compute on the surface
of a torus than the surface of a sphere.
Science ¬ction writers have described toroidal planets (Boyd, 1981c, 1984c) and plasma
physicists solve ¬‚ows in tori to model the ring-shaped fusion generators known as toka-
maks (Schnack et al., 1984). In this section, we will discuss a two-dimensional basis for the
surface of a torus. However, one may de¬ne a three-dimensional orthogonal coordinate
system (“toroidal coordinates”, Morse and Feshbach, 1953) in which the third coordinate
is radial distance from the surface of a torus to its centerline.


. 80
( 136 .)