Chebyshev coef¬cients for the function known variously as the “Lorentzian” or the “Witch

of Agnesi”, which is

∞

2

(2.107)

= an Tn (x)

(x2 + 2)

n=0

2.16. BOUNDARY LAYER RULE-OF-THUMB 57

where the odd coef¬cients are zero by symmetry and where the even coef¬cients are asymp-

totically

|an | ∼ 2 e’n (2.108)

n 1& 1

For simplicity, assume is small. Bounding the error by the sum of all the neglected coef-

¬cients as in (2.95) shows that the error in truncating after aN is approximately |aN |/ for

1. If we demand the truncation error = 10’d so that d is the number of decimal places

of accuracy, then we obtain the following.

Rule-of-Thumb 5 (WITCH-OF-AGNESI RESOLUTION ESTIMATE)

If a function has an internal boundary layer of width 2 due to poles at x = ± i , then the value

of N needed to reduce the truncation error for that function to 10’d is roughly

d

N≈ (2.109)

log(10) 1

What this rule-of-thumb asserts is that for more realistic functions ” with singularities

” the slope of the error will asymptote to a straight line on a log/linear plot as shown in

Fig. 2.20 rather than diving more steeply off the bottom of the graph, as true for the entire

function sin(x).

By “width of the internal boundary layer”, we mean that the spoor of the complex

singularities is that the witch-of-Agnesi function rises from 1/2 to 1 and then falls to 1/2

again on the small interval [’ , ] as shown on the inset in Fig. 2.20. Presumably, some

small-scale features in hydrodynamic ¬‚ows (or whatever) are similarly associated with

poles or branch points for complex x close to the center of the feature.

A couple of caveats are in order. First, singularities a distance from the real axis are

most damaging to convergence when near the center of the interval. As one moves√ to

x = ±1, the pole becomes progressively less harmful until N becomes a function of 1/

” not 1/ ” in the immediate neighborhood of the endpoints. This special endpoint effect

is described in the next section, but it signi¬cantly modi¬es (2.109) only when the branch

point or pole is quite close to x = ±1. (Observe the equiconvergence ellipses in Fig. 2.16;

the innermost curves (small ) are very eccentric and approximately parallel the real axis

except in the small, semi-circular caps around ±1.)

Second caveat: Eq.(2.108) should be multiplied by a factor of nk when the singularities

are other than simple poles where k = 1 for a second order pole (stronger singularity) and

k = ’1 for a logarithm (weaker singularity, the integral of a simple pole). In this case, one

can show that N (d, ) should be determined from the more complex rule

N + k log(N ) ≈ d log(10) (2.110)

Often, however, especially when d is large, the logarithm of N may be neglected. Usually

the simpler rule is best because one rarely knows the precise type or location of a complex

singularity of the solution to a differential equation; one merely estimates a scale, .

2.16 Boundary Layer Rule-of-Thumb

Viscous ¬‚ows have boundary layers next to solid surfaces where the tangential velocity

is reduced to zero. When the boundary layer is of width where 1, it is normally

necessary to increase the number of grid points as O(1/ ) as ’ 0. This is extremely

expensive for high Reynolds number ¬‚ows.

CHAPTER 2. CHEBYSHEV & FOURIER SERIES

58

0

10 1

0.8

0.6

0.4

-5

10 0.2

2µ

0

-0.5 0 0.5 1

-10

10

-15

10

10 15

0 5

N µ / log(10)

Figure 2.20: Witch-of-Agnesi Rule-of-Thumb. This predicts that if a function has a peak in

the interior of the interval of width 2 due to singularities at (x) = , then the error will be

10d where d = N /log(10). The errors in the Chebyshev approximation of f (x) ≡ 2 /(x2 +

2

) is plotted versus N / log(10), i. e., against the expected number of decimal digits in

the error for varying from 1 (circles) to 1/100. As predicted, the curves lie approximately

on top of one another (with some deviation for = 1. When N / log(10) ≈ 15, the error

is indeed O(10’15 ). The inset graph illustrates the “Witch of Agnesi” function itself along

with the meaning of the width parameter .

A much better approach is to use a mapping, that is, a change of variable of the form

(2.111)

y = f (x)

where y is the original, unmapped coordinate and x is the new coordinate. Unfortunately,

there are limits: If the mapping varies too rapidly with x, it will itself introduce sharp

gradients in the solution, and this will only aggravate the problem. However, the higher

the order of the method, the more rapidly the transformation may vary near the bound-

ary. Orszag and Israeli (1974) assert (without giving a proof) that with an “optimal map-

ping”, this change-of-coordinates trick can reduce the error (for a given total number of

grid points) by O( 2 ) for a second order method and O( 4 ) for a fourth order scheme.

One effective mapping for the most common case ” boundary layers near the end-

points ” was employed by Kalnay de Rivas(1972) who used

(2.112)

y = cos(x)

for a “physical” coordinate y ∈ [’1, 1]. This is precisely the change-of-variable that con-

verts a Chebyshev series in y into a Fourier series in x. It follows that her evenly spaced

grid in x is the same grid that we would use with an un-mapped Chebyshev series in y. To

put it another way, the uneven grid in y that is sometimes such a burden with Chebyshev

methods, particularly for time-integration schemes, is in fact just what we need to resolve

boundary layers.

2.16. BOUNDARY LAYER RULE-OF-THUMB 59

0

10 1

µ

0.8

0.6

0.4

-1 0.2

10

0.97 0.98 0.99 1

-2

10

-3

10

3

0 1 2 4 5

Νµ

Figure 2.21: Boundary Layer Rule-of-Thumb. The errors in the approximation of f (x) =

tanh([x ’ 1]/ ) are plotted versus Nscaled ≡ N 1/2 for = 1 (lower circles), = 1/10 (as-

terisks), = 1/100 (plus signs), = 1/1000 (x™s) and = 1/10000 (upper circles). Except for

= 1, which is offset a little below the others, the curves for different are superimposed

one atop another. Regardless of the thickness of the boundary layer, when Nscaled ≥ 3,

that is, when the number N of Chebyshev polynomials is greater than 3/ 1/2 , the approx-

imation has at least moderate (1%) accuracy. The inset shows f (x) = tanh((x ’ 1)/ ) for

= 1/100, together with the grid points (circles) for Nscaled = 3, i. e., N = 30. As long

as there are a couple of grid points in the boundary layer, the Chebyshev approximation is

okay.

Rule-of-Thumb 6 (BOUNDARY LAYER RESOLUTION REQUIREMENT)

To resolve a boundary layer of width at y = ±1 (as opposed to an internal boundary layer,

such as a front), one should use approximately

3

N=√ polynomials (2.113)

1

to obtain a moderate accuracy of O(1 %).

Source: Orszag & Israeli (1974).

This rule-of-thumb, like all such empirical principles, must always be applied with cau-

tion. (See Fig. 2.21.) It is interesting, however, that it agrees exactly with Kalnay™s (1972)

rule-of-thumb for her second order ¬nite difference method, even to the numerical factor

of 3, although hers was expressed in terms of grid points rather than polynomials.

The reason that N is proportional to the square root of 1/ , instead of 1/ itself, is that

the Chebyshev polynomials have a “built-in” variable stretching which is “quadratic” in

the sense that the interior grid points closest to the boundary are only O(1/N 2 ) away from

the endpoint [rather than O(1/N ) as for an evenly spaced grid]. This is easily proved by

recalling that the Chebyshev grid is the map of an evenly spaced grid in x via y = cos x,

and then using the Taylor expansion of the cosine near x = 0, π.

CHAPTER 2. CHEBYSHEV & FOURIER SERIES

60

The conclusion is that with Chebyshev polynomials and Fourier series, the best way

to deal with boundary layers may well be to do nothing! (When Fourier series are used,

the boundary condition is spatial periodicity; a periodic interval is equivalent to a ring-

shaped interval without boundaries or boundary layers.) For some problems, however, the

natural, quadratic boundary layer stretching which is, as it were, “built-in” to Chebyshev

methods may not be optimum. Then, an additional mapping may be helpful.

For boundary value and eigenvalue problems, it follows that boundary layers are not

a major problem. For time-dependent problems, alas, they are a serious complication be-

cause the narrow grid spacing near the wall implies that a very small time step is needed

with an explicit time-differencing scheme. This is not a problem peculiar to spectral meth-

ods, however; Kalnay™s ¬nite difference scheme was subject to exactly the same restrictions

as a Chebyshev method, and for the same reason.

Thus, boundary layers pose no special problems for spectral methods. Rather, the

Chebyshev grid is unfortunate only when boundary layers are not present. With Cheby-

shev polynomials, the high density of grid points near the walls requires either a very

small timestep or an implicit time integration even when there is no physical reason for

high resolution near the endpoints.

Table 2.2: Relationships Between Singularities and Asymptotic Spectral Coef¬cients.

Note: all rows apply to Fourier series as well as Chebyshev and Legen-

dre polynomials except where marked by (*), which exclude trigonometric series.

Form of singularity Type of singularity Asymptotic form

spectral coef¬cients

1/(x ’ a) [ ]pn

Simple pole

1/(x ’ a)2 [ ]n pn

Double pole

[ ]n’1 pn

log(x ’ a) Logarithm

√

[ ]n’1/2 pn

1/ x ’ a Reciprocal of square root

(x ’ a)1/3 [ ]pn /n4/3

Cube root

[ ] exp(’p n1/2 )

In¬nitely differentiable but

exp(’q/|x|)

singular at x = 0

xψ [ ] 1/nψ+1

Branch point within [-1,1]

Jump discontinuity

f (x) = sign(x) [ ] /n

[ ] /n2

Discontinuous ¬rst derivative

df /dx = sign(x)

(assuming f continuous)

(1 ’ x)ψ [ ] 1/n2 ψ+1 *

Branch point at endpoint

[ ] exp(’p n2/3 )

In¬nitely differentiable but

exp(’q/|x + 1|)

singular at endpoint

Chapter 3

Galerkin Method, the Mean

Weighted Residual Method &

Inner Products

“Six months in the lab can save you a day in the library.”

Albert Migliori, quoted by J. Maynard in Physics Today 49, 27 (1996)

3.1 Mean Weighted Residual Methods

One difference between algorithms is the “error distribution” principle, that is, the means

of making the residual function R(x; a0 , a1 , . . . , aN ) small where

R(x; a0 , a1 , . . . , aN ) ≡ HuN ’ f (3.1)

where H is the differential or integral operator and uN (x) is the approximation

N

uN ≡ (3.2)

an φn (x)

n=0

for some suitable basis functions φn (x).