. 95
( 136 .)


REDUCE2, required only 90 seconds on a CDC-7600 ... and corrected six errors in the
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


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)

Galerkin or Methods of Moments
Chebyshev polynomials
Legendre or Gegenbauer or Powers
(i)Fourier and rational Chebyshev =’ Chebyshev
(ii)Transcendental data and forcing =’ Polynomial
(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-
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

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
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;
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)
= 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:
x2j R4 (x; a0 , a1 , a2 , a3 )dx = 0, (20.7)
j = 0, 1, 2, 3

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,


. 95
( 136 .)