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

equation.

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.

CHAPTER 20. SYMBOLIC CALCULATIONS

472

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

N

(21.2)

an Tn (±1) = 0

n=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,

(21.3)

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

473

CHAPTER 21. THE TAU-METHOD

474

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

(21.5)

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

(21.6)

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

(21.7)

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

21.2. „ -APPROXIMATION FOR A RATIONAL FUNCTION 475

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

N

fN (x) ≡ (21.9)

fn Tn (x)

n=0

then the fn are determined by solving the linear equations

N

(21.10)

(Ti , Tj Q) fj = (Ti , P ) i = 0, . . . , N

i=0

where

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

Q(x)

j=N +1

CHAPTER 21. THE TAU-METHOD

476

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

(21.13)

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