Fourier or Chebyshev coefficients

-4

10

-6

10

-8

10

TB

-10

10

-12

10

Weideman-Cloot

-14

10

-16

10

0 20 40 60 80 100

j

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:

Rule-of-Thumb 16 (OPTIMIZING INFINITE INTERVAL MAP PARAMETER)

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

CHAPTER 17. METHODS FOR UNBOUNDED INTERVALS

378

-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

systems.

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.

380

18.2. POLAR, CYLINDRICAL, TOROIDAL, SPHERICAL 381

θ=π/2

θ=π/4

θ=3π/4

r

θ=0

θ=π

θ

θ=7π/4

θ=5π/4

θ=3π/2

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

lines).

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

(18.1)

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

CHAPTER 18. SPHERICAL & CYLINDRICAL GEOMETRY

382

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