. 15
( 136 .)


cally the “Witch of Angesi” Rule-of-Thumb™ because it is derived by examing the analytical
Chebyshev coef¬cients for the function known variously as the “Lorentzian” or the “Witch
of Agnesi”, which is

= an Tn (x)
(x2 + 2)

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

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

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

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.

10 1



10 0.2

-0.5 0 0.5 1


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

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

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.

10 1

-1 0.2
0.97 0.98 0.99 1


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

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

N=√ polynomials (2.113)

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, π.

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
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
uN ≡ (3.2)
an φn (x)

for some suitable basis functions φn (x).


. 15
( 136 .)