unknown f (u, v) as

(16.26)

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.

CHAPTER 16. COORDINATE TRANSFORMATIONS

330

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

function.

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

mapping.

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,

16.7. PERIODIC PROBLEMS & THE ARCTAN/TAN MAP 331

275000

225000

175000

125000

75000

25000

-25000

’π/4 π/4

’π/2 π/2

0

y

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

1

(16.28)

tan(y)

x = arctan

L

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

(16.29)

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.

CHAPTER 16. COORDINATE TRANSFORMATIONS

332

0

-2

No Map

-4

-6

-8

Map with

L=1/2

-10

-12

-14

-16

-18

5

0 10 15 20 25 30

N

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:

(16.30)

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

(16.31)

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

16.8. ADAPTIVE METHODS 333

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 ’ γ}

1

Bayliss&Matkowsky(1987) β+1

arctan(±(1+γ))

Bayliss&Matkowsky(1992) β = arctan(±(1’γ))

» = [(β + 1)/2] arctan (±(1 ’ γ))

4 parameters: ± [width] γ [location], β, »

β’x

1 + π arctan ± tan π βx’1 ’ 1

4

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)

0

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).

CHAPTER 16. COORDINATE TRANSFORMATIONS

334

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