. 69
( 136 .)


sides of the polygon. Mason is able to impose the boundary conditions by writing the
unknown f (u, v) as

f (u, v) = ¦(u, v) g(u, v)

where ¦(u, v) = 0 on all sides of the polygon. He then expands g(u, v) as a double Cheby-
shev series with no additional boundary constraints imposed upon it. He was able to
successfully calculate the lowest eigenvalue to ¬ve decimal places by computing the corre-
sponding eigenvalue of an 81—81 matrix.

Mason™s work is an excellent tour de force, but the use of a single mapping for the whole
domain is obviously complicated and demands both analytical and numerical ingenuity.
For this reason, Mason™s work has had few direct imitations. However, the Liverpool
school of “global element” enthusiasts has had great success in applying two-dimensional
mappings to parts of the domain. Separate Chebyshev series are then applied in each of
several subdomains (McKerrell, 1988, and Kermode, McKerrell, and Delves, 1985). This
work will be discussed in more detail in the chapter on domain decomposition methods.
An alternative to mapping is to augment the basis with singular functions chosen to
match the behavior of the corner singularities. Lee, Schultz, and Boyd (1988a) have applied
this to the ¬‚ow in a rectangular cavity; Boyd (1988d) gives one-dimensional illustrations.
Much earlier, this technique has been combined with ¬nite element methods as reviewed
by Strang and Fix (1973). Although a rigorous theory is lacking, the method of singular
basis functions seems to work well. The only caveat is that one must use only a small
number of singular basis functions. The reason is that very weak branch points such as
xk log(x) where k is large can be approximated to many decimal places of accuracy by the
N Chebyshev polynomials in the basis. Thus, the corresponding singular basis functions
are indistinguishable from a sum of Chebyshev polynomials on a ¬nite precision machine,
leading to poorly conditioned matrices. The moral is that one should use singular basis
functions only for mimicing those singularities which are so strong that the non-singular
(Chebyshev) part of the basis set cannot cope with them.
Both two-dimensional mappings and singular basis functions require precise knowl-
edge about the nature of the corner branch points. The exponential mapping is a brute
force technique which simply concentrates lots and lots of grid points near the boundaries.
It is thus the preferred method when the branch points are (i) too strong to be ignored and
(ii) the analytical form of the singularities is not known.

16.7 Periodic Problems with Concentrated Amplitude and
the Arctan/Tan Mapping
Many physical problems have solutions which are concentrated in a small portion of the
computational interval. The “cnoidal” waves of the Korteweg-deVries equation are a good
example. When the amplitude is small, this solution is indistinguishable from the cosine
function. As the amplitude increases, however, the crest of the wave becomes increasingly
tall and narrow. In the limit of in¬nite amplitude, the cnoidal wave has the shape of a delta
A change of coordinate which clusters the points of the pseudospectral grid in and
near the peak is a good strategy for such strongly “concentrated” solutions. Boyd (1978c)
describes a cubic-plus-linear stretching in conjunction with Chebyshev polynomials. Boyd
(1987c) describes similar changes-of-coordinate for use when the solution is periodic. The
latter article is worth discussing in greater detail because it illustrates how to choose a good
For functions that are periodic in y with a period of π, Boyd (1987c) suggests

x, y ∈ [0, π] (16.27)
y = arctan(L tan(x))

Here, y is the original coordinate, x is the new numerical variable, and L is a constant (“map
parameter”) which can be chosen to give as little or as much stretching as desired. This
transformation has several advantages. First, the mapping function can be expressed in
terms of elementary transcendentals, that is, functions built-in to most compilers. Second,







’π/4 π/4
’π/2 π/2
Figure 16.1: The “bicnoidal” wave of the ¬fth-degree Korteweg-deVries equation (Boyd,
1986b). Although the wave is spatially periodic, it is so sharply peaked about the center of
the interval and decays so rapidly (exponentially) towards the endpoints that it is wasteful
to use a Fourier spectral method with an evenly spaced grid.

the map has the equally trivial inverse

x = arctan

This explicit inverse is valuable in evaluating the numerical solution at arbitrary values
of the original coordinate y. Third, the derivative of y(x) is a trigonometric polynomial
in x. This, together with (16.28), implies that differential equations whose coef¬cients are
trigonometric polynomials in y will be transformed into differential equations whose co-
ef¬cients are still trigonometric polynomials in the numerical coordinate x ” the map is
“polynomial coef¬cient preserving”. Fourth, the mapping is periodic in the sense that

y(x) = x + P (x; L)

where P (x; L) is a periodic function of x whose Fourier series is given explicitly in Boyd
(1987c). Fifth, and perhaps most important, the mapping (16.27) is smooth so that it does
not create additional boundary layers or regions of rapid variation.
Figs. 16.1 and 16.2 illustrate the effectiveness of the mapping. Fig. 16.1 shows the graph
of the so-called “bicnoidal wave” of the Fifth-Degree Korteweg-deVries equation. Since
the function is graphically indistinguishable from 0 over half of the interval, it is obviously
wasteful to use an evenly spaced grid. Fig. 16.2 shows the logarithm of the relative error
(relative to the maximum of the wave) as a function of the number of terms in the Fourier
pseudospectral approximation for L = 1 (no change-of-coordinates) and L = 0.5. The
mapping allows one to achieve a given accuracy with roughly half the number of degrees
of freedom ” and 1/8 the work ” as without the mapping.
No Map


Map with



0 10 15 20 25 30

Figure 16.2: The logarithm of the relative error as a function of the number N of Fourier
cosine functions without the use of a mapping and with the application of the arctan / tan
map with L = 1/2. The change-of-coordinate allows one to achieve the same accuracy with
roughly half the number of basis functions.

16.8 Adaptive Methods
Adaptive algorithms, which automatically compress and stretch the grid so as to resolve
boundary layers, shocks, and other regions of rapid change, have become a “hot topic” in
every area of numerical analysis. Table 16.2 catalogues spectral examples.
As a prototype, consider Burger™s equation with periodic boundary conditions:

ut + u uy + ν uyy = 0

Because of the damping, the solutions are always single-valued and never break, but for
small ν, u(x, t) may develop a moving layer of very steep gradient, a quasi-shock wave. An
obvious strategy is apply a transformation so that an evenly spaced grid in the computa-
tional coordinate x is mapped into a grid which has many grid points clustered around the
shock. Because the front is propagating and the equation is nonlinear, it is almost impossi-
ble to choose a mapping non-adaptively. From scale analysis, we might be able to estimate
the width of the frontal zone, but we cannot predict its exact location without ¬rst solving
the problem.
We will discuss two different strategies for adaptation. The ¬rst is to choose arclength s
along the solution curve u(x) as the new coordinate (White, 1982). This is always de¬ned
even if we set ν = 0; u may become a multi-valued function of x, but it is always a single-
valued function of arclength.
The second strategy is to choose a particular functional form for the mapping and re-
strict adaptation to the more modest goal of choosing the parameters. For example, the
arctan/tan mapping is

y = yf + arctan[L tan(x)]

Table 16.2: A Selected Bibliography of Adaptive Pseudospectral Methods With and With-
out Mappings

References Comments
Bayliss&Matkowsky(1987) Arctan/tan mapping; combustion fronts
Guillard&Peyret(1988) Flame fronts and Burgers™ equations
Bayliss&Gottlieb Arctan/tan map
&Matkowsky&Minkoff(1989) reaction/diffusion problems
Augenbaum (1989) Adaption, based on Sobolev norm error indicators
for both Fourier and Chebyshev algorithms
Augenbaum(1990) Multidomain adaptation; acoustic waves in discontinuous media
Bayliss&Kuske&Matkowsky(1990) Parameters of arctan/tan map are allowed to vary
with the other spatial coordinate
Bayliss et al.(1994) thermo-viscoplastics
Bayliss&Matkowsky(1992) Arctan/tan map for nonlinear cellular ¬‚ames
Mavriplis(1994) Adaptive meshes for spectral elements
Gauthier et al.(1996) Moving interface with domain decomposition
Carcione(1996) Adaptive mapping for elastic wave equation
Renaud&Gauthier(1997) Both map and subdomain walls are moved
adaptively; compressible vortex roll-up
Mulholland,Huang&Sloan(1998) Adaptive ¬nite difference scheme de¬nes map
for pseudospectral solution for problems with frontal zones
Black(1998) spectral elements; map of semin¬nite domain into ¬nite interval

Table 16.3: A Few Mapping Functions
Note: y is the physical (unmapped) coordinate; x is the computational coordinate.
References Change-of-Coordinate Functions
y = β’1 + » arctan {±x ’ γ}
Bayliss&Matkowsky(1987) β+1
Bayliss&Matkowsky(1992) β = arctan(±(1’γ))
» = [(β + 1)/2] arctan (±(1 ’ γ))
4 parameters: ± [width] γ [location], β, »
1 + π arctan ± tan π βx’1 ’ 1
Bayliss&Gottlieb 4
& Matkowsky&Minkoff(1989) 2 parameters: ± [width] β [location]
± > 0, ’1 < β < 1
Guillard&Peyret(1988) Flame fronts and Burgers™ equations
y = (2/π) arctan ± π X
Augenbaum (1989, 1990) 2
X ≡ (γ + x)/(γ + 1)
2 parameters: ±[width] γ[location]

where yf and L are adaptively varied with time. The width parameter L(t) allows the
zone of high resolution to shrink or expand as needed. The shift parameter yf (t) allows
the center of the high resolution zone to track the front.
A simple tactic for determining yf and L is to demand that they minimize a “smooth-
ness” function such as
I(yf , L) ≈ |uxx |2 + |ux |2 dx (16.32)

It is easy to solve such minimization problems via Newton™s method, but a ¬rst guess is
needed. In a time-integration, we can use the known structure of the initial solution to
generate a ¬rst guess. At later times, we use the optimum parameters for time level j as
the ¬rst guess for minimizing the functional (16.32) at time level (j + 1). We then apply
(16.31) to generate the grid that will be used to march to time level (j + 2).

The extra costs are obvious: We must integrate over the interval (to generate the func-
tional) and then interpolate (with spectral accuracy!) from the old grid to the new. One
harsh reality which has discouraged pseudospectral applications is that the Fast Fourier
Transform cannot be applied to interpolate from the old grid to the slightly stretched-and-
shifted new grid. We must use the slower “off-grid” interpolation methods described in
Chapter 10, Sec. 7. (A useful Rule-of-Thumb is that these schemes, whether Boyd™s or
multipole-based, cost roughly four times as much as an FFT of the same N .)
There are some good coping strategies. For example, one might shift into a coordinate
system moving at the approximate phase velocity of the front. One would then need to
update the mapping only every ten time steps, instead of every step, because the front will
not move very far relative to the center of the high-resolution zone in that short length of
time. Even so, pseudospectral adaptation is a rather costly business.
The philosophy is the same, however: to minimize the residual with respect to width-
and-shift parameters so as to optimize the pseudospectral solution. Adaptive grids are
very promising, but because of the extra overhead, it is important that the minimization be
carried out only occasionally as we march in time.

16.9 Chebyshev Basis with an Almost-Equispaced Grid:
Kosloff/Tal-Ezer Arcsine Map
Kosloff and Tal-Ezer(1993) and Tal-Ezer (1994) have reduced the severe time-step require-
ments of standard Chebyshev methods by using a mapping which distributes the grid


. 69
( 136 .)