. 21
( 136 .)


N +1
TN +1 (x) ≡ (x ’ xi ) (4.12)
2N i=1

To minimize the error in the Cauchy Remainder Theorem, the polynomial part of the re-
mainder should be proportional to TN +1 (x). This implies that the OPTIMAL INTERPO-

Figure 4.4: Graphical construction of the unevenly spaced Chebyshev interpolation grid. If
a semicircle of unit radius is cut into evenly spaced segments, and then vertical lines are
drawn from these “Fourier” gridpoints to the line segment [’1, 1], which is the base of the
semicircle, the vertical lines will intersect the horizontal line at the Chebyshev grid points.
The polar coordinates of the grid points on the semicircle are unit radius and angle θ =
π(2i ’ 1)/(2N ) where i = 1, . . . , N . The Chebyshev grid points are xi = cos(θi ).

Since the Chebyshev polynomials are just cosine functions in disguise, these roots are

(2i ’ 1)π
xi ≡ ’ cos (4.13)
i = 1, 2, . . . , N + 1
2(N + 1)

By Taylor expansion of the cosine function, we verify the assertion made above that the
grid spacing is O(1/N 2 ) near the endpoints:

π2 9π
x1 ≈ ’1 + x2 ≈ ’1 + (4.14)
; [N 1]
8N 2 8N 2
and similarly near x = 1.
The grid is illustrated in Fig. 4.4.
A number of important questions still remain. First, how accurate is Chebyshev in-
terpolation? Does it converge over as wide a region as the usual Chebyshev expansion
in which we compute a polynomial approximation by integration instead of interpolation?
Second, many problems (such as those in spherical coordinates or an in¬nite domain) re-
quire something exotic like spherical harmonics or Hermite functions. What interpolation
scheme is best for those geometries? Third, it is obvious that (4.13) de¬nes a choice of in-
terpolation points that we can use to solve differential equations via the pseudospectral
method, but what is the relationship, if any, between such an algorithm and the Galerkin
In the next section, we generalize the choice of a best grid to other basis sets appropri-
ate for other geometries. In later sections, we will return to the Chebyshev and Fourier
methods and show that the price we pay for interpolation instead of integration is modest

4.3 Gaussian Integration & Pseudospectral Grids
The reason that “collocation” methods are alternatively labelled “pseudospectral” is that
the optimum choice of the interpolation points, such as (4.13) for Chebyshev polynomials,

makes collocation methods identical with the Galerkin method if the inner products are
evaluated by a type of numerical quadrature known as “Gaussian integration”.
Numerical integration and Lagrangian interpolation are very closely related because
one obvious method of integration is to ¬t a polynomial to the integrand f (x) and then
integrate PN (x). Since the interpolant can be integrated exactly, the error comes entirely
from the difference between f (x) and PN (x). The standard textbook formulas are of the
f (x) dx ≈ (4.15)
wi f (xi )
a i=0

The weight functions wi are given by
wi ≡ [“quadrature weight”] (4.16)
Ci (x) dx

where the Ci (x) are the cardinal functions on the set of points {xi } as de¬ned by (4.8) above.
As if there were not already enough jargon, the interpolation points are “abscissas” and
the integrals of the cardinal functions are “quadrature weights” in the context of numerical
The order of a numerical method is directly related to the highest polynomial for which
it is exact. For example, the usual three-point and ¬ve-point centered difference formulas
(Chapter 1) are exact when applied to general second-degree and fourth-degree polynomi-
als, respectively, so the errors are O(h2 ) and O(h4 ), respectively. Similarly, a quadrature
formula with (N + 1) points will be exact if the integrand is a polynomial of degree N .
Gauss made the observation that an evenly spaced grid is nothing sacred. In particular,
if we allow the interpolation points {xi } as well as the weights {wi } to be unknowns, we
have twice as many parameters to choose to maximize the accuracy of the method, and we
can make it exact for polynomials of degree (2N + 1).

If the (N + 1) “interpolation points” or “abscissas” {xi } are chosen to be the zeros of PN +1 (x)
where PN +1 (x) is the polynomial of degree (N + 1) of the set of polynomials which are orthogonal
on x ∈ [a, b] with respect to the weight function ρ(x), then the quadrature formula
f (x) ρ(x) dx = wi f (xi )
a i=0

is exact for all f (x) which are polynomials of at most degree (2N + 1).

PROOF: Davis (1975, pg. 343)
As noted earlier, we can take almost any smooth function ρ(x), insert it into the integral
inner product, and then apply Gram-Schmidt orthogonalization to 1, x, x2 , . . . to create an
in¬nite set of polynomials which are orthogonal with respect to the inner product
(f, g) = f (x) g(x) ρ(x) dx

The Chebyshev, Legendre, Gegenbauer, Hermite, and Laguerre polynomials merely corre-
spond to different choices of weight functions and of the interval [a, b].
The general theory of orthogonal polynomials shows that the (N + 2)-th member of a
family of polynomials is always of degree (N + 1) and always has exactly (N + 1) real zeros

within the interval [a, b]. We might worry about the possibility that, for exotic basis sets like
Hermite functions, the zeros might be complex, or otherwise outside the physical domain,
but this fear is groundless. The only dif¬culty is that the roots are not known in closed form
for general N except for Chebyshev polynomials or trigonometric interpolation. However,
the roots and weights wi for various N and various weight functions ρ(x) can be found in
most mathematical handbooks such as Abramowitz and Stegun (1965) and also are easily
calculated by the FORTRAN routines in the appendix of Canuto et al.(1988). The Legendre-
Lobatto roots up to degree N are known in analytical form and are given in Appendix F,
Sec. F.8.

The “Composite Trapezoidal Rule” and “Composite Midpoint Rule”, which both use abscissas
which are uniformly spaced, are Gaussian quadratures in the sense that these formulas are exact
with N points for TRIGONOMETRIC polynomials of degree 2N ’ 2, that is, for polynomials
composed of the constant plus the ¬rst (N ’ 2) cosines plus the ¬rst (N ’ 2) sines.

PROOF: Implicit in the discrete orthogonality relationships between the sines and cosines
given in the next section.
The Trapezoidal Rule and Midpoint Rule are described in introductory numerical anal-
ysis courses as simple but crude quadrature schemes with an accuracy of O(1/N 2 ). This is a
libel, true only when integrating a NON-PERIODIC f (x). For PERIODIC functions, these nu-
merical integration schemes are GAUSSIAN QUADRATURES with an error which decreases
GEOMETRICALLY FAST with N for f (x) which are analytic for real x.
What does all this have to do with interpolation and the pseudospectral method? The
answer is that since quadrature formulas are obtained by analytical integration of the in-
terpolating polynomial, it follows that the best choice of quadrature points is also the best
choice of interpolation points and vice-versa. This implies the following principle that ap-
plies to Fourier series, Chebyshev, Legendre and Gegenbauer polynomials, Hermite and
Laguerre functions and so on.

The grid points should be the abscissas of a Gaussian quadrature associated with the basis set.

The vague phrase “a Gaussian” quadrature is necessary because there are actually two
useful Gaussian quadratures associated with each basis. In the Fourier case, one can use
either the Trapezoidal Rule or the Midpoint Rule; which are equally accurate. The Trape-
zoidal Rule, as used in the next section, is the default choice, but the Midpoint Rule (alias
Rectangle Rule), which does not include either x = 0 or x = 2π as grid points, is convenient
when solving differential equations which have singularities at either of these endpoints.
Lobatto showed that one could obtain a numerical integration formula with all the good
properties of Gaussian quadrature by using the zeros of the ¬rst derivative of a Chebyshev
or Legendre polynomial plus the endpoints x = ±1 as the quadrature points. This “Lo-
batto” or “endpoints-plus-extrema” grid is now more popular than the “Gauss” or “roots”
grid for solving boundary value problems because the boundary conditions instantly fur-
nish two grid point values for the unknown u(x). When the differential equation is singular
at the endpoints ” for example, in cylindrical coordinates, differential equations usually
have singularities at r = 0 even when the solution is perfectly smooth and well-behaved
at the origin ” then the “roots” grid is preferable . There is little difference between the
Trapezoidal and Midpoint Rule, or between the Gauss and Lobatto grids, in terms of accu-
racy or other theoretical properties.
The one mild exception to this Rule-of-Thumb is that weather forecasting models, which
use a spherical harmonics basis, usually employ a single latitudinal grid for all zonal

wavenumbers. The spherical harmonics for a given zonal wavenumber m are Gegenbauer
polynomials whose order is m + 1/2. Technically, the optimum grid for a given m is a func-
tion of m. However, using different grids for different m would entail enormous amounts
of interpolation from one grid to another to calculate the nonlinear terms where different
wavenumbers m interact. So, a single grid is employed for all m. However, this grid is the
Gauss grid for the Legendre polynomials, so interpreted broadly, the rule applies to this
case, too.

4.4 Pseudospectral: Galerkin Method via Gaussian Quadra-
The integrals that de¬ne the matrix elements in the Galerkin method usually must be eval-
uated by numerical quadrature. What scheme is most ef¬cient? The answer is Gaussian
quadrature, using the roots of the very same orthogonal polynomials (or trigonometric
functions) that are the basis functions. This choice connects the collocation and Galerkin
Theorem 16 (GALERKIN™s with QUADRATURE)
If the matrix elements are evaluated by (N +1)-point Gaussian quadrature, then for linear prob-
lems the spectral coef¬cients {an } in
uN, G (x) = an φn (x)

as calculated by the Galerkin method will be identical with those computed by the collocation method
using the same number of terms in the series and the same set of interpolation points.
It is because of this close relationship between the collocation and Galerkin methods that collocation-
at-the-Gaussian-quadrature-abscissas is also known as the “pseudospectral” method. (Historically,
before the development of collocation algorithms, “spectral” was a synonym for Galerkin.)

PROOF: Let the pseudospectral matrix equation be L a = f . The matrix elements are
Lij ≡ L φj (xi ) fi ≡ f (xi ) (4.20)
where L (without vector symbol) is the operator of the differential or integral equation and
f (x) is the forcing function. A single row of the matrix equation is then
Lij aj = fi

or written out in full,
[L φj (xi )] aj = f (xi )

Multiply each row of (4.22) by wi φn (xi ) for some particular choice of n and then add.
(For simplicity, assume the basis functions are orthonormal.) The constants wi are the Gaus-
sian quadrature weights. This gives
Hnj aj = gn

Hnj ≡ (4.24)
wi φn (xi ) [L φj (xi )]

But this is simply the Galerkin matrix element

Hnj ≈ (φn , L φj ) (4.25)

if the inner product integral is evaluated using Gaussian quadrature with (N + 1) points.
gn ≡ (4.26)
wi φn (xi ) f (xi )

gn ≈ (φn , f ) (4.27)

which is the n-th row of the column vector on the R. H. S. of the Galerkin matrix equation,
H a = g.
If some rows of the matrix are used to impose boundary conditions, then the proof is
only slightly modi¬ed. The boundary rows are the same for either the Galerkin or pseu-
dospectral methods. With two boundary rows, for example, we impose only (N ’ 1) collo-
cation conditions. We then add those collocation rows using the weights wi appropriate to
Gaussian integration with (N ’ 1) points, and again reproduce the corresponding (N ’ 1)
rows of the Galerkin matrix. Q. E. D.
Theorem 16 implies that collocation ” with the right set of points ” must inherit the
aura of invincibility of the Galerkin method.

If a set of (N + 1) basis functions are orthogonal under the integral inner product

(φi , φj ) = δij

where δij is the Kronecker δ-function, then they are still orthogonal with respect to the discrete inner


. 21
( 136 .)