. 97
( 136 .)


# Pick off the coef¬cient of A1xx, A3xx, and call them div1, div3;
div1:=coeff(eq1,A1xx,1); div3:=coeff(eq3,A3xx,1);
eq1:=expand(simplify(eq1/div1)); eq3:=expand(simplify(eq3/div3));
eq1f:=evalf(eq1); eq3f:=evalf(eq3);
20.3. EXAMPLES 471

and similarly at higher order. By applying the usual Galerkin method, we obtain (for N =
2) the set of two coupled nonlinear equations in two unknowns:

1 1 1
’(3 + 1 )A1 ’ 1.09A2 ’ 1.45A1 A3 ’ 19.8A2 = 0
A1,xx 1 3
c c c c
1 1 1
’(7 + 1 )A3 ’ 5.93A2 ’ 1.65A1 A3 ’ 0.0302A2 = 0 (20.30)
A3,xx 3 1
c c c c
The Maple listing is Table 20.6.
This example demonstrates one of the strengths of algebraic manipulation: Correction
of errors in hand calculations. The set (20.30) was originally given in Boyd (1989e), but
with the coef¬cient of A1 A3 in the second equation equal to 0.82 instead of the correct
value shown above. The error was found only when the derivation was repeated using
computer algebra.
When (20.30) is further approximated by setting one or the other of the unknowns to
zero, the result is the travelling-wave form of the Korteweg-deVries equation. The general
solution of this N = 1 approximation is an elliptic function; limiting cases are in¬nitesimal
cosine waves (which recovers (20.28)) and solitary waves.
When both modes are retained, the result is the system (20.30) which cannot be analyti-
cally solved. However, “maximally simpli¬ed” models like the two-mode system are often
very useful in illuminating complicated behavior. E. N. Lorenz employed this philosophy
in modelling two-dimensional convection by a system of three ordinary differential equa-
tions which had been previously derived by B. A. Saltzman. Lorenz™ discovery of chaos
in this system, which now bears his name, was a great milestone in the theory of chaotic
dynamics, strange attractors, and fractals. Similarly, Boyd (1989e) shows that a two-mode
model is the smallest which displays the “weakly nonlocal” solitary waves of the full AEW
The single-mode approximation, which is just the KdV equation, is usually obtained
by a very different route which employs the singular perturbation technique known as
the “method of multiple scales” (Bender and Orszag, 1978, Boyd, 1980c). The end result,
though, is identical with applying the spectral method with a single basis function. The
perturbative analysis is really just a systematic way to answer the question: When is a
one-mode approximation justi¬ed? The answer is that when the solitary wave is very
wide with a length scale which is O(1/ ) and an amplitude of O( 2 ) where << 1 is a
small parameter, then the nonlinear terms coupling A1 (x) with A3 (x) will be small. The
self-interaction term cannot be neglected, however, if the corresponding linear term is also
small, i. e., of O( 4 ). This happens when the phase speed differs from the k = 0 (“long
wave”) limit of the linear phase speed by an amount of O( 2 ). It is also necessary to use a
spectral basis whose functions have the structure of the in¬nitesimal waves, i. e., Hermite
functions for this example.
The formal justi¬cation for small is heartening, but experience with the method of
multiple scales (and few-mode models) is that both are often qualitatively accurate even
when ∼ O(1), that is, far outside their region of formal accuracy (Bender and Orszag,
1978, Finlayson, 1973). The symbolic spectral method is a very quick and general alterna-
tive to the perturbation theory.
The only caveat is one illustrated by all the examples above: To get the most out of
the spectral method (or any other solution algorithm), it is important to understand the
physics. Does u(x) have a boundary layer? Is it smooth so that small N will be okay, or does
it have a complicated structure that will require very large N (and probably make symbolic
calculations unfeasible)? Is some choice of basis functions, such as Hermite functions for
the AEW equation, particularly well-matched to the solutions? Computer algebra and
numerical calculations are not alternatives to analysis but merely its extension.

20.4 Summary and Open Problems
There has been considerable progress in marrying spectral methods with algebraic ma-
nipulation languages as summarized here and in the review, (Boyd, 1993). None of the
guidelines should be interpreted too rigidly; a couple of examples deliberately broke some
of the precepts of earlier sections to emphasize the need for common sense.
The major open problem is coping with nonlinearity. The examples above show some
success, but only for fairly easy problems. Perturbation expansions or a single Newton
iteration may be essential in obtaining an answer which is simple enough to be useful.
Fox and Parker (1968) summarize a large body of hand calculations which exploited
recurrence relations of the Chebyshev polynomials instead of Galerkin™s method. These
recurrences signi¬cantly improve accuracy for small N . This line of attack has been contin-
ued by Ortiz and his students as reviewed in Chapter 21. The Galerkin method is simpler
and more general, so we have omitted these alternative Chebyshev algorithms from this
discussion. However, these recurrences lend themselves well to symbolic computations.
Algebraic manipulation is still dominated by a “perturbation-and-power series” men-
tality. However, there is much more to symbolic mathematics than perturbation series or
the derivation of exact algebraic relationships. The main algorithms of this work, spec-
tral methods and Newton™s iteration, are still badly underutilized in symbolic calculations.
Perhaps some readers will change this.
Chapter 21

The Tau-Method

“I never, never want to be a pioneer . . . It™s always best to come in second, when you can
look at all the mistakes the pioneers made ” and then take advantage of them.”
” Seymour Cray

21.1 Introduction
The tau-method is both an algorithm and a philosophy. It was invented by Cornelius Lanc-
zos in the same (1938) paper that gave the world the pseudospectral method.
As an algorithm, the tau-method is a synonym for expanding the residual function as
a series of Chebyshev polynomials and then applying the boundary conditions as side
constraints. Thus, it is indistinguishable from what we have earlier described as Galerkin™s
method with “boundary bordering”. The formal de¬nition is the following.

De¬nition 43 (Tau-Method) A mean-weighted residual method is a TAU-METHOD if the “test
functions” are the Chebyshev polynomials and the inner product is the usual Chebyshev integral
inner product.

For example, to solve a second order ODE with u(’1) = u(1) = 0, the „ -method would
impose the constraints

m = 0, . . . , N ’ 2 (21.1)
(Tm , R(x; a0 , . . . , aN ) = 0
an Tn (±1) = 0

where the {an } are the Chebyshev coef¬cients of u(x) and R(x; a0 , . . . , aN ) is, as usual,
the residual function obtained by substituting the truncated Chebyshev series into the dif-
ferential equation.
In contrast, Canuto et al. (1988) prefer to apply the label “Galerkin” only to that mean
weighted residual method that chooses the “test” functions to be basis functions that satisfy
the homogeneous boundary conditions. For example,

(φm , R(x; b2 , . . . bN ) = 0 m = 2, . . . , N


where the basis functions are

φ2n (x) ≡ T2n (x) ’ 1 φ2n+1 (x) ≡ T2n+1 (x) ’ x (21.4)

These two methods do give slightly different answers because the two highest-degree con-
straints in (21.3) involve the inner product of the residual with TN ’1 (x) and TN (x) whereas
the highest Chebyshev polynomial multiplying the residual in the tau-method inner prod-
uct (21.1) is only TN ’2 (x). For this reason, this terminological distinction is popular in
the literature (Gottlieb and Orszag, 1977, for example). However, the accuracy differences
between the “tau” and “Galerkin” methods are likely to be negligible.
At the end of the chapter, we will return to issues of nomenclature.
First, however, we shall discuss the philosophy of the „ -method. Lanczos recognized
that there are two ways of attacking a differential equation. The obvious approach is to
compute an approximate solution to the exact, unmodi¬ed differential equation. The sec-
ond is to compute the exact solution to a modi¬cation of the original differential equation.
If the “modi¬cation” is small, then the solution to the modi¬ed problem will be a good
approximation to that of the original problem.
This second strategy ” to solve approximate equations exactly ” is the philosophy of the
„ -method. In the next section, we shall apply this tactic to approximate a rational function.
In Sec. 3, we shall extend this philosophy to differential equations. Finally, in Sec. 4, we
shall discuss Lanczos™ “canonical polynomials”, which have been the jumping-off point
for a long series of papers by E. L. Ortiz and his collaborators.

21.2 „ -Approximation for a Rational Function
Suppose the goal is to approximate a rational function f (x) by a polynomial fN of degree
N . Let

f (x) = Pp (x)/Qq (x)

where P (x) and Q(x) are polynomials of degree p and q, respectively. Substituting f ’ fN
and multiplying through by Q(x) gives

fN (x) Q(x) = P (x)

which has only polynomials on either side of the equal sign. At ¬rst glance, it appears
that we could compute the (N + 1) coef¬cients of fN (x) merely by matching powers of
x. However, the degree of fN (x)Q(x) is (N + q); thus, we have more conditions on the
coef¬cients of fN (x) than we have unknowns to satisfy them.
Lanczos observed that there is an exact, polynomial solution to the perturbed problem

fN (x) Q(x) = P (x) + (x)

where (x) is a polynomial of degree (N + q) with q undetermined coef¬cients. Eq. (21.7)
has polynomials of degree (N + q) on both sides with a total of (N + q + 1) unknowns.
Matching the powers of x in (21.7) gives a set of (N + q + 1) linear equations to determine
the unknown coef¬cients of fN (x) and (x).
The secret to success is to choose (x) ” more accurately, to choose N + 1 of the (N +
q + 1) coef¬cients of (x) ” such that the perturbation is small in some appropriate sense.
The obvious choice is to set all j = 0 for j ¤ N . In this case, fN (x) is the usual power
series approximation to f (x).The trouble with this choice is that (x) is highly non-uniform;
because it is O(xN +1 ) for small |x|, this “power series perturbation” is very, very small near

the origin and then grows rapidly for larger x. The result is that (i) the accuracy of fN (x)
is similarly non-uniform (good for small |x| and increasingly bad as |x| increases) and (ii)
the error may be huge. Indeed, if |x| is larger than the absolute value of any of the roots
of the denominator, Q(x), then the approximation (which is just a truncated power series)
will diverge as N ’ ∞ even if f (x) is bounded and smooth for all real x.
Lanczos™ second key observation is that the Chebyshev polynomials oscillate as uni-
formly as possible on their canonical interval, [-1, 1]. It follows that if we de¬ne
N +q
[“Lanczos perturbation”] (21.8)
(x) = „j Tj (x)
j=N +1

the perturbation will spatially uniform in magnitude on x ∈ [’1, 1] and so will the error
|f (x) ’ fN (x)|.
Of course, the error may still be large if the coef¬cients „j are large. However, observe
that the coef¬cients of (x) are simply those of the smooth function fN (x) Q(x) ’ P (x). We
saw in Chapter 2 that the coef¬cients of any well-behaved function fall off exponentially
fast with N (for suf¬ciently large N ), so that it follows that the „j will be exponentially
small, too, at least for N 1.
If we write
fN (x) ≡ (21.9)
fn Tn (x)

then the fn are determined by solving the linear equations
(Ti , Tj Q) fj = (Ti , P ) i = 0, . . . , N


„j = (Tj , P ’ fN Q) (21.11)
j = N + 1, . . . , (N + q)

Several points are important. First, it is not necessary to explicitly compute the „ -
coef¬cients in order to determine the approximation fN (x); (21.11) is a second, separate
step that can be applied only after fN (x) has already been computed.
Second, the „ -approximation is not simply the Chebyshev expansion of f (x); the solu-
tion of the linear system (21.10), whose matrix elements are inner products of polynomials,
is usually different from the coef¬cient integrals of the Chebyshev series of f (x), which
are the inner products of the Chebyshev polynomials with the rational function f (x). As
N ’ ∞, of course, the differences between the „ -approximation, the Chebyshev series,
and the pseudospectral interpolant decrease exponentially fast with N .
Third, (21.10) is not the way Lanczos himself performed the calculation. In his time B. C.
(Before Computers), it was more convenient to represent fN (x) as an ordinary polynomial,
expand the Chebyshev polynomials in (x) as powers of x, and then match the powers
of x. This saved work, but it required simultaneously determining both fN (x) and (x).
Lanczos™ variant obscures the fact that it is possible to determine fN (x) independently of the
perturbation (x) if one represents all quantities as Chebyshev polynomials.
Fourth, the „ -coef¬cients are useful for a posterior error analysis because
N +q
(x) 1
f (x) ’ fN (x) = ¤ |„ | (21.12)
min |Q(x)|
j=N +1

if we can ¬nd a lower bound on Q(x) on x ∈ [-1, 1]. In practice, this error analysis, although
highly praised in Fox & Parker (1968), is not useful in the age of microcomputers: it is easier
to simply evaluate the difference between f (x) and fN (x) for a large number of points and
take the maximum.
Furthermore, the „ -method is not normally the method of choice even for generating a
Chebyshev approximation to a given f (x). Numerically evaluating the inner products of
f (x) with the polynomials, i. e. the usual Chebyshev series, is easy and avoids the costly
inversion of a matrix. The „ -method is only useful in conjuction with algebraic manipula-
tion methods for small N where one wants an approximation with rational coef¬cients or
where f (x) may contain symbolic parameters, making (21.11) more attractive.
Nevertheless, Lanczos™ philosophy of exactly solving a perturbed problem is useful in
other areas of applied mathematics. The technique that evolved from that philosophy is
still useful even today for solving differential equations.

21.3 Differential Equations
It is not possible to solve

ux + u = 0 & u(’1) = 1


. 97
( 136 .)