N +1

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-

LATION POINTS are the ROOTS of the CHEBYSHEV POLYNOMIAL of DEGREE (N + 1).

CHAPTER 4. INTERPOLATION, COLLOCATION & ALL THAT

86

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

method?

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

indeed.

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,

4.3. GAUSSIAN INTEGRATION & PSEUDOSPECTRAL GRIDS 87

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

form

N

b

f (x) dx ≈ (4.15)

wi f (xi )

a i=0

The weight functions wi are given by

b

wi ≡ [“quadrature weight”] (4.16)

Ci (x) dx

a

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

integration.

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

Theorem 14 (GAUSS-JACOBI INTEGRATION)

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

N

b

(4.17)

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

b

(4.18)

(f, g) = f (x) g(x) ρ(x) dx

a

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

CHAPTER 4. INTERPOLATION, COLLOCATION & ALL THAT

88

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.

Theorem 15 (PERIODIC GAUSSIAN QUADRATURE)

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.

Rule-of-Thumb 7 (CHOICE OF PSEUDOSPECTRAL GRIDS)

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

4.4. PSEUDOSPECTRAL IS GALERKIN METHOD VIA QUADRATURE 89

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-

ture

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

methods.

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

N

(4.19)

uN, G (x) = an φn (x)

n=0

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

N

(4.21)

Lij aj = fi

j=0

or written out in full,

N

(4.22)

[L φj (xi )] aj = f (xi )

j=0

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

N

(4.23)

Hnj aj = gn

j=0

CHAPTER 4. INTERPOLATION, COLLOCATION & ALL THAT

90

where

N

Hnj ≡ (4.24)

wi φn (xi ) [L φj (xi )]

i=0

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.

Similarly,

N

gn ≡ (4.26)

wi φn (xi ) f (xi )

i=0

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.

Theorem 17 (ORTHOGONALITY under the DISCRETE INNER PRODUCT)

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

(4.28)

(φi , φj ) = δij

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

product