published manual derivation.”

“ R. A. d™Inverno (1976)

20.1 Introduction

Computer algebra has had a strange and checkered history. The notion of manipulating

symbols instead of merely numbers is as old as the concept of the computer itself as shown

by the quote from the Countess of Lovelace. Algebraic manipulation (AM) languages such

as REDUCE and MACSYMA have been available for more than twenty years. And yet the

spread of applications of algebraic manipulation has been, to use the punning words of a

1979 reviewer, “a case of creeping ¬‚ow”. Two decade later, the ¬‚ow is still creeping, but

with a noticeably smaller viscosity.

Several problems have retarded the growth of “AM” languages. First, symbolic cal-

culations are best done interactively with a live human watching the stream of formulas

materialize on the screen, and intervening as needed with substitution and simpli¬cation

rules. Symbolic manipulation was therefore a rather and slow painful business on screen-

less, teletype terminals. Second, AM programs are terrible memory hogs; many research

461

CHAPTER 20. SYMBOLIC CALCULATIONS

462

projects were prematurely truncated through running out of storage. Third, symbolic ma-

nipulation is inherently an order-of-magnitude more expensive than numerical manipula-

tion: a number is a number is a number (with apologies to Gertrude Stein), but a symbolic

variable must be tagged with all sorts of baggage to indicate precisely what the symbol

symbolizes. Any bold soul who wants to push the state of the art in calculating perturba-

tion series will sadly turn away from Maple and Mathematica to slog through the agonies

of writing a FORTRAN program to do the same job “ the numerical program will be ten

times as long, but it may run 100 times faster and require only a tenth the storage to reach

a given perturbation order.

Paradoxically, however, it is as certain as tomorrow™s sunrise that algebraic manipula-

tion will triumph in the end. The scientist of the twenty-¬rst century will use AM every-

day, and will no more be able to imagine life without Maple or Mathematica or their ilk

than his counterpart of the 30™s could imagine designing without a slide rule. It seems un-

likely, after twenty years of waiting, that algebraic manipulation will ever become the major

computational tool for more than a tiny fraction of scienti¬c projects. As more and more

engineers have AM on their desks, however, algebraic manipulation will likely become an

indispensable tool for small, “scratchpad” calculations of the kind we do every day.

It follows that symbolic spectral codes will always remain few compared to spectral

models that manipulate only numbers. Nevertheless, there are several reasons why spec-

tral algorithms and algebraic manipulation languages are an especially favorable pairing.

First, interpolatory and non-interpolatory spectral methods both manipulate polyno-

mials and trigonometric functions, which are precisely what REDUCE, Maple, and Mathe-

matica multiply, divide, differentiate, and integrate most easily. Second, spectral methods

have always been used at low order to compute analytical solutions to simple problems (Fin-

layson, 1973). Computer algebra greatly extends the “analytic range” of spectral methods

to more complex problems and higher resolution. Third, a major use of symbolic manip-

ulation has been to compute perturbation series. In many applications, the perturbation

series is a spectral series (Rand, 1984; Rand and Armbruster, 1987; Mills, 1987; Haupt and

Boyd, 1988, 1991).

In the remainder of this chapter, we discuss the striking differences between numeri-

cal and symbolic spectral methods. When calculating by hand or by REDUCE, one must

apply very different rules to choose the most ef¬cient method, the optimum strategy. In-

deed, many of the precepts for numerical calculations must be reversed or inverted when

the number of degrees of freedom is small and the goal is not a table or a graph but an

algebraic expression. The next two sections discuss such “inversions”. The fourth section

is a catalogue of the little tricks and devices that are important in writing ef¬cient pro-

grams. Sec. 5 describes several examples. The ¬nal section is a summary of the current

state-of-the-art and lists several open problems.

20.2 Strategy

Table 20.1 summarizes four strategies for symbolic calculations. All spectral methods sub-

stitute a series with unknown coef¬cients aj into the differential equation to obtain the

so-called residual function, R(x; a0 , a1 , . . . , aN ’1 ) where N is the number of basis func-

tions in the approximation. If the trial solution uN (x) were exact, then the residual would

be zero for all x. With a truncated spectral approximation, the best we can do (barring

special cases) is to make the residual as small as possible in some sense.

The most widely used discretization strategy for numerical applications is known var-

iously as “collocation”, “orthogonal collocation”, “discrete ordinates” or the “pseudospec-

tral” method. One obtains N equations in N unknowns by demanding that the residual

20.2. STRATEGY 463

Table 20.1: Precepts for spectral methods: Numerical (N >> 1) versus Symbolic (N small)

SYMBOLIC [N ∼ O(3)] NUMERICAL [N >> 1]

Collocation

Galerkin or Methods of Moments

Chebyshev polynomials

Legendre or Gegenbauer or Powers

Unnecessary

Polynomialization

(i)Fourier and rational Chebyshev =’ Chebyshev

(ii)Transcendental data and forcing =’ Polynomial

Rationalization

(i)Collocation at rational points (i) Collocation at roots

of TN , which are irrational

(ii)Replace irrational numbers (ii) Unnecessary

&¬‚oating point numbers by rational numbers

function should be zero at each of N interpolation points. In contrast, the Galerkin method

demands that the residual should be orthogonal to each of the ¬rst N basis functions. The

accuracy difference between the two discretization strategies is negligible for large N , so

collocation is usually preferred because of it requires fewer computations and less pro-

gramming (Chapter 4).

For small N , however, Galerkin™s method may be two or three times more accurate than

collocation. To put it another way, the rule of thumb quoted in Chapter 3, Sec. 9, is that

the Galerkin solution with N = m is usually as accurate as the pseudospectral solution

with N = m + 1 or m + 2, provided that the solution is suf¬ciently smooth that a small

N is suf¬cient to resolve it. Since the complexity of a symbolic answer usually grows

exponentially with N , we must reverse the usual precept, and suggest that for calculations

in Maple or REDUCE, etc., one should use Galerkin™s method wherever possible.

Unfortunately, when a problem is cursed with an exponential nonlinearity as for the

Bratu equation, uxx + exp(u) = 0, Galerkin™s method may be completely unfeasible. How-

ever, it may be still be possible to get a simple but accurate answer by using collocation, as

in Boyd (1986c), Finlayson (1973) and our Example Two.

In ¬‚oating point computations, it is almost always necessary, unless N is very small,

to represent both the basis functions and the weight functions as orthogonal polynomials.

Because algebraic manipulation languages compute in exact, rational arithmetic wherever

possible and also because N is very small, one may use the equivalent powers of x as

both basis functions and as the Galerkin test functions, in which case Galerkin™s method is

usually called the “method of moments”. Numerical ill-conditioning, which is a perennial

worry in high resolution number-crunching, is almost never a worry in symbolic manipu-

lation.

A second issue is the choice of the basis functions. As explained in Chapter 2, the

Chebyshev and Legendre polynomials are part of a large family (“Gegenbauer polynomi-

als”) which are orthogonal on the interval [-1, 1] with weight function w(x) = (1 ’ x2 )± .

The Chebyshev polynomials are the special case ± = ’1/2; the Legendre are ± = 0, i. e., no

weighting at all. The Chebyshev weighting, heavily biased towards the endpoints, is op-

timum for expanding general functions; the Legendre error is smaller than the Chebyshev

over most of the interval, but is much larger (roughly by O(N 1/2 ) ) near the endpoints as

illustrated in Sec. 13 of Chapter 2.

The solution to a boundary value problem is not arbitrary, however; its boundary be-

havior is strongly constrained by the boundary conditions. Therefore “ especially for an

approximation with a very small number of degrees of freedom “ it is preferable to replace

Chebyshev polynomials by Legendre polynomials or perhaps even Gegenbauer polyno-

mials of larger ± (Finlayson, 1973). If ± is a positive integer, there is the further bene¬t that

CHAPTER 20. SYMBOLIC CALCULATIONS

464

the Galerkin integrals are simpler because of the absence of the square root of the Cheby-

shev weight, increasing the probability that the integration routines, built-in to all major

algebraic languages, will be able to perform the integrations analytically.

The third strategic suggestion has no counterpart in numerical calculations: It is to

convert all transcendental functions to polynomials. A Fourier cosine series, for example,

can be converted into a Chebyshev series by making the change of variable t = arccos(x)

and using the identity Tn (x) = cos(nt) for all n. The rational Chebyshev functions, which

are a basis for in¬nite and semi-in¬nite domains (Chapter 17), are also just Chebyshev

polynomials after a change of variable. Similarly, powers of sech(x), which often arise in

soliton theory, can be “polynomialized” by de¬ning z = tanh(x) and applying the identity

sech2 (x) = 1 ’ tanh2 (x) = 1 ’ z 2 . Lastly, transcendental coef¬cients of a differential

equation can be approximated by polynomials. The error in

cos([π/2]x) ≈ 1 ’ 0.4967(πx/2)2 + 0.03705(πx/2)4 , (20.1)

for example, is less 0.0009 on x ∈ [’1, 1].

It is often convenient to delay “polynomialization” until an intermediate stage. Maple,

for example, can differentiate and multiply trigonometric polynomials and then convert

them to a standard Fourier series through the “combine(expression,trig)” command. How-

ever, it has no command to pick off the coef¬cients of a Fourier series, so polynomializa-

tion through a command like “subs(z=cos(x),z**2=cos(2*x), z**3=cos(3*x), expression)” is

essential to exploiting the “coeff” command, which does return a particular coef¬cient of a

polynomial.

Algebraic manipulation languages are more ef¬cient at calculations with integers and

rational numbers than ¬‚oating point numbers, so the fourth strategy is “rationalization”,

that is, converting ¬‚oating point numbers to rational numbers. The usual collocation

points, for example, are the irrational roots or extrema of the (N + 1)-st basis function.

When collocation, as opposed to the preferred Galerkin method, is necessary, it is usually

better to approximate an irrational collocation point such as (1/5)1/2 = 0.44721 . . . by a

rational number such as 4/9 = 0.444 . . . . For small N , the error is small. Other ¬‚oating

point numbers in the boundary conditions, differential equation, etc., should usually be

“rationalized” too.

Arnold (1983) shows that for any irrational number µ, there exists a sequence of rational

approximations, p/q, which converge to µ. The error of a given member of the sequence is

less than the reciprocal of the square of the denominator, i. e.,

|µ ’ p/q| < 1/q 2 (20.2)

The well-known approximation π ≈ 355/113, which has a relative error of less than 10’6 ,

was long thought to be very special. What (20.2) shows, in words, is that for any µ, there ex-

ists an approximation with a three-digit denominator and six-digit accuracy. The sequence

of optimum rational approximations can be computed by an algorithm known variously as

the “Euclidean algorithm”, “algorithm of continued fractions”, or “the algorithm of stretch-

ing the noses” (Arnold, 1983). The ¬‚oating point coef¬cients in (20.1) can be approximated

by 75/151 and 1/27, for example, with errors of only 0.00011 and 0.000014, respectively.

There are two major exceptions to “rationalization”. First, when the answer depends

only a single irrational such as 21/2 , it may be more ef¬cient to carry a symbolic variable

ROOT2 through the remainder of the algebra than to replace it by a rational number. The

second is that rationalization may be unnecessary if the later stages of the computation

have to be done numerically, as often happens. It is not, alas, unusual to see hundred

digit numbers when a long calculation is performed in exact arithmetic (“integer swell”) as

shown in Table 20.5 below. The precepts in Table 20.1 should be interpreted as suggestions

and guidelines rather than commandments. Let common sense be the highest law!

20.3. EXAMPLES 465

Table 20.2: Maple program for Example One: Linear BVP with Boundary Layer

u:= (1-x*x) * (a0 + a1*x*x + a2*x**4 + a3*x**6);

R4:=tau*tau*diff(u,x,x) - u + 1;

eq0:=integrate(R4,x=-1..1); eq1:=integrate(R4*x*x, x=-1..1);

eq2:= integrate(R4*x**4, x=-1..1); eq3:= integrate(R4*x**6, x=-1..1);

solve( eq0, eq1, eq2, eq3, a0, a1, a2, a3 );

A couple of other important pieces of advice have been omitted from the table because

they are not inversions of usual practice. First, exploiting parity and other symmetries

(Chapter 8) is always highly advantageous in numerical calculations. The sort of simple,

idealized problems that are amenable to symbolic spectral methods are precisely the kind

of problem that are most likely to have symmetries. Therefore, the standard precept: Look

for Symmetry! is redoubled for symbolic computations.

Second, one may impose boundary conditions through either “boundary-bordering”

or “basis recombination” as explained in Chapter 6. In symbolic calculations, the latter is

strongly preferred because incorporating the boundary conditions into the basis functions

may halve the number of degrees of freedom, and enormously reduce the complexity of the

algebraic formula which is the answer. Furthermore, it is particularly easy to incorporate

homogeneous boundary conditions into a symbolic answer: for homogeneous Dirichlet

conditions: Multiply the trial solution by (1 ’ x2 ).

20.3 Examples

All examples were programmed in REDUCE or Maple. The choice of languages has no

particular signi¬cance except that both are widely used and the author knows them.

Example One: Linear ODE BVP (Carrier & Pearson, 1968)

uxx ’ u = ’1;

2

(20.3)

u(’1) = u(1) = 0

where is a constant. The exact solution is

u(x) = 1 ’ cosh(x/ )/ cosh(1/ ) (20.4)

Because u(x) is symmetric, we shall assume an expansion in even powers of x. Fur-

ther, we incorporate the homogeneous boundary condition into the series. Arbitrarily, we

choose N = 4:

u4 (x) = (1 ’ x2 ) a0 + a2 x2 + a4 x4 + a3 x6 (20.5)

Substituting (20.5) into (20.3) gives the residual function

≡ u4,xx ’ u4 + 1 (20.6)

R4

= 1 ’ [1 + 2 2 ]a0 ’ 2 2 a2 + a0 ’ [1 + 12 2 ]a1 + 12 2 a2 x2

+ a1 ’ [1 + 30 2 ]a2 + 30 2 a3 + a2 ’ [1 + 56 2 ] 2 a3 x6 + a3 x8

Because u(x) has a boundary layer of width O( ) “ all the variation of u(x) is close to the

endpoints “ we shall not include a factor of (1 ’ x2 ) in the Galerkin conditions, but instead

demand that the unweighted integrals of R4 with each of the ¬rst four even powers of x

should be zero:

1

x2j R4 (x; a0 , a1 , a2 , a3 )dx = 0, (20.7)

j = 0, 1, 2, 3

’1

CHAPTER 20. SYMBOLIC CALCULATIONS

466

Table 20.3: The maximum pointwise errors in the four-term Legendre approximation as a

function of the parameter for Example One: Linear BVP with Boundary Layer

L∞ error

1/20 0.081

3/40 0.030

1/10 0.012

3/20 0.0046

1/5 0.00040

1/4 0.000099

3/10 0.000028

4/10 0.0000033

1/2 0.00000056

3/4 1.7E-8

1 1.3E-9

This gives four equations in four unknowns. This can be cast as a matrix equation of the

form H A = F and solved in REDUCE by the single statement A:=(1/H)*F. A REDUCE

program for general 2d order BVPs is Table 17-1 of the ¬rst edition of this book (Boyd,