where ck = 2 if k = 0 and ck = 1 for k > 0.

The one ¬‚aw is that there are still two expensive procedures that must be executed at

each step: (i) the calculation of the coef¬cients {ak } from the grid point values {u(xj )}

(interpolation) and (ii) the evaluation of the derivative from the series (9.15) at each of

the N grid points (summation). If we execute these two steps the easy way ” by matrix

multiplication ” both cost O(N 2 ) operations, and we are worse off than before.

What makes the hybrid pseudospectral method ef¬cient is that we can replace the

costly, rate-determing matrix multiplications by the Fast Fourier Transform (FFT) algo-

rithm. This requires only O(5N log2 N ) real operations to compute either grid point values

from coef¬cients or vice versa. The bad news is that the FFT is applicable only to Fourier

series and to its offspring-by-change-of-variable, Chebyshev polynomials; this is a major

reason for the primacy of Fourier and Chebyshev expansions over all their competitors.

The FFT has been the foundation for the modern rise of spectral methods and is there-

fore the theme of the next chapter.

Fig. 9.1 is a schematic of the pseudospectral algorithm. The ¬gure shows clearly that we

jump back and forth between the grid point representations of u(x) [to multiply by c(x)]

and the spectral coef¬cients [to differentiate u(x)].

9.3 The Shamrock Principle

A millenium and a half ago, a Welsh bishop named Patrick was chie¬‚y responsible for

converting the Irish people from paganism to Christianity. According to tradition, the doc-

trine of the Trinity was an especially hard one to understand: One God, yet three Persons:

Father, Son and Holy Spirit.

CHAPTER 9. EXPLICIT TIME-INTEGRATION METHODS

178

St. Patrick devised an ingenious audiovisual aid: a little plant called the shamrock,

rather like a clover but with squarish leaves, which grows wild all over the country. Pluck-

ing a shamrock, he said, “Look! It is one plant, but it has three leaves. Three can be One.”

Pseudospectral time-marching require a similar conceptual awakening:

Principle 1 (SHAMROCK PRINCIPLE) In the limit N ’ ∞, the set of spectral coef¬cients

{an } and the set of gridpoint values {u(xj )} are fully equivalent to each other and to the function

u(x) from whence they come. These three different, seemingly very different representations, are

really just a single function.

Fig. 9.2 shows a shamrock with each leaf appropriately labeled. The grid point values

and the spectral coef¬cients of a function u(x) and the function itself are all equivalent to

one another. We must freely use whichever form is most convenient.

u(x)

Analytical

form

{a n} {u(xj)}

n=0,...,N j=0,...,N

Spectral Grid Points

Coefficients Values

Summation

Interpolation

Figure 9.2: A visualization of the Shamrock Principle, or the trinity of spectral methods.

The ¬rst N Chebyshev coef¬cients, {an }, and the set of N values on the interpolating grid,

{u(xj )}, are equivalent to the analytical de¬nition of u(x), represented by the top leaf, in the

sense that either of these two sets allows one to compute u(x) to arbitrarily high precision

in the limit N ’ ∞.

Via interpolation, one may pass from the “grid point representation”, {u(xj )}, to the “spec-

tral” representation, {an }. Summation of the Chebyshev series transforms {an } ’ {u(xj )}.

9.4 Linear and Nonlinear

Reviews of spectral methods in hydrodynamics (e. g., Orszag and Israeli, 1974) often justify

the pseudospectral algorithm by showing the inef¬ciency of Galerkin™s method for nonlin-

9.5. EXAMPLE: KDV EQUATION 179

Table 9.1: MATLAB code for Right-Hand Side of the System of ODEs in Time

function F=KdVRHS(t,u); global k kcub

a=ifft(u); % Compute Fourier coef¬cients an from grid point values u(xj )

ax=k .* a; axxx = kcub .* a; % Compute coef¬cients of 1st and 3rd derivative

ux=real(fft(ax)); uxxx=real(fft(axxx) ); % Reverse FFT to get grid point values

% of ¬rst and third derivatives

F= - u .* ux - uxxx; % RHS of KdV ODE system. Nonlinear term evaluated by

%pointwise multiplication of u by ux

% In a preprocessing step, either in the main program or in a subroutine called once,

% one must execute the following to initialize the vectors k and kcub with the product

% of i with the wavenumber k and with the (negative) of its cube, respectively

for j=1:(n/2), k(j)=-i*(j-1); end; for j=(n/2+1):n, k(j)=-i*(j-1 - n); end

for j=1:n, kcub(j)=k(j)*k(j)*k(j), end

ear equations. This is a little misleading because the real problem is not the nonlinearity

per se, but rather the fact that coef¬cients of the differential equation vary with x.

Consequently, Sec. 2 deliberately discussed a linear equation. The methodology is ex-

actly the same for a nonlinear problem.

9.5 The Time-Dependent KdV Equation: An Example

To illustrate explicit time-marching with spectral discretization in space, examine the Korteweg-

deVries (KdV) equation:

(9.19)

ut + u ux + uxxx = 0

subject to the boundary condition of spatial periodicity with period 2π:

(9.20)

u(x + 2π) = u(x)

for all x. If the grid point values of u are used for time-advancement, as in a ¬nite difference

code, the system of ordinary differential equations in time is

u(xj , t)

= Fj (t) ≡ ’u(xj , t) ux (xj , t) ’ uxxx (xj , t) (9.21)

dt

A complete code must create the initial conditions, declare parameters, loop forward

in time and call the Runge-Kutta subroutine. However, all these steps are necessary with

¬nite difference or ¬nite element discretizations of x, too. The only part of the overall code

that needs any knowledge of spectral methods is the subroutine that computes the right-

hand side of the system of ODEs in time, that is, the vector F . A ¬nite difference code can be

upgraded to spectral accuracy merely by replacing the subroutine which calculates F .

Table 9.1 shows a Fourier pseudospectral subroutine (in the language MATLAB) for

evaluating this vector of grid point values for the KdV equation. It calls “fft” and “ifft”,

which are built-in MATLAB Fast Fourier Transform routines. It employs MATLAB™s oper-

ation for elementwise multiplication of one vector by another vector, which is denoted by

“.*”. A FORTRAN 77 subroutine would be a little longer because the elementwise multi-

plications would have to be replaced by DO loops and so on. Nevertheless, the brevity of

the subroutine is startling.

CHAPTER 9. EXPLICIT TIME-INTEGRATION METHODS

180

The ¬rst line computes the Fourier coef¬cients, a vector a, by taking a Fourier Trans-

form. The second line computes the coef¬cients of the ¬rst and third derivatives by multi-

plying those of u(x) by (ik) and (’ik 3 ), respectively. (It is assumed that in a preprocessing

step, the vectors “k” and “kcub” have been initialized with the appropriate wavenumbers).

The third step is to take inverse Fourier Transforms to convert these coef¬cients for the

derivatives into the corresponding grid point values. The ¬nal step is to add the grid point

values together to form the vector F . Note that the nonlinear term is evaluated by point-

by-point multiplication of the grid point values of u with those of ux .

Tables 9.2 and 9.3 are short bibliographies of the simplest class of time-marching prob-

lems: PDEs with periodic boundary conditions, solved by Fourier series.

Table 9.2: Time-Marching with a One-Dimensional Fourier Basis

References Comments

Abe&Inoue (1980) Korteweg-deVries (KdV) nonlinear wave equation

Chan&Kerkhoven (1985) Explicit & semi-implicit schemes for KdV eqn.

™

Fla (1992) Energy-conserving scheme for Derivative-NLS equation

Fornberg(1975) Hyperbolic equations

Fornberg&Whitham(1978) Nonlinear waves.: KdV, MKdV, Benjamin-Ono & others

Fornberg (1987) Elastic waves; comparisons with ¬nite differences

Frutos et al. (1990) “Good” Boussinesq;

Hamiltonian (symplectic) time-marching

Frutos&Sanz-Serna (1992) 4-th order time integration for KdV & other wave eqns.

Garc´

ia-Archilla(1996) ˜Equal Width˜ equation

Guo&Manoranjan (1985) Regularized Long Wave (RLW) equation

Herbst&Ablowitz (1992) Sine-Gordon eqn.; numerical instabilities;

integrable-to-chaos transition because of numerical errors

Herbst&Ablowitz (1993) Symplectic time-marching, numerical chaos,

exponentially small splitting of separatrices

¨

If&Berg&Christiansen Split-step spectral for Nonlinear Schrodinger Eq.

& Skovgaard(1987) with absorbing (damping) boundary conditions

Li&Qin (1988) Symplectic time-marching

Ma&Guo (1986) Korteweg-deVries equation; uses “restrain operator”

Mulholland&Sloan(1991) Filtering and its effects on solitons for RLW, KdV, etc.

Mulholland&Sloan(1992) Implicit & semi-implicit with preconditioning

for wave equations

Nouri&Sloan(1989) Comparisons between pseudospectral algorithms for KdV

Qin&Zhang (1990) Multi-stage symplectic schemes for Hamiltonian systems

Salupere&Maugin KdV; many-soliton ¬‚ows

&Engelbrecht&Kalda(1996)

Sanders&Katopodes KdV; RLW, Boussinesq eqs.

&Boyd (1998)

Sanugi&Evans (1988) Nonlinear advection equation

Sloan (1991) Regularized Long Wave (RLW) eqn.

Wang (1991) Hamiltonian systems and conservation laws

Weideman&James (1992) Benjamin-Ono equation

Zheng&Zhang&Guo(1989) SLRW equation

9.6. IMPLICITLY-IMPLICIT: RLW & QG 181

Table 9.3: A Selected Bibliography of Time-Marching with Multi-Dimensional Fourier

Bases

References Comments

Three-dimensional turbulence: 5123 resolution

Chen et al.(1993)

Fox&Orszag(1973) Two-dimensional turbulence

Fornberg (1977) Two-dimensional turbulence in doubly periodic box

Ghosh&Hossain 2D MHD turbulence (non-quadratic nonlinearity)

&Matthaeus(1993)

Haidvogel(1977, 1983) Review of quasi-geostrophic models

Hald(1981) Navier-Stokes equations

Hua&Haidvogel(1986) Quasi-geostrophic turbulence; details of algorithms

Hua(1987) Review of quasi-geostrophic ocean models

Kosloff&Reshef Elastic waves (seismology)

&Loewenthal(1984)

Ma&Guo(1987a) Two-dimensional vorticity equation

Tan&Boyd(1997) Two-dimensional generalization of quasi-geostrophic eq.

Boyd&Tan(1998) Solitary vortices, topographic deformations

Rashid&Cao Low Mach number compressible ¬‚uid ¬‚ow

&Guo(1993, 1994a,b)

Vallis(1985) Doubly-periodic quasi-geostrophic ¬‚ow

9.6 Implicitly-Implicit: the Regularized Long Wave equa-

tion and the Quasi-Geostrophic equation

Some partial differential equations have time derivatives multiplied by a differential oper-

ator L, i. e.,

(9.22)

L ut = G(u, x, t)

To restate this in the canonical form

(9.23)

ut = F (u, x, t),

we must solve a boundary value problem for F :

(9.24)

LF =G

For such an equation, it is as laborious to apply an explicit method as an implicit method

because one must solve a boundary value problem at every time step. Hence, we shall

dub such equations “implicitly-implicit” because the computational labor of an implicit

time-marching algorithm is inevitable.

Two examples are the Regularized Long Wave (RLW) equation of water wave theory,

{1 ’ ‚xx }ut = ’ux ’ uux (9.25)

and the quasi-geostrophic equation of meteorology and physical oceanography,

{‚xx + ‚yy }ψt = ψy (ψxxx + ψxyy ) ’ ψx (ψxxy + ψyyy ) (9.26)

When the boundary conditions are spatial periodicity, the inversion of the operator L is

trivial for both examples.

For the RLW equation, for example, the ¬rst step is to evaluate

G(xj , t) = ’ (1 + u(xj ) ) ux (xj ) (9.27)

CHAPTER 9. EXPLICIT TIME-INTEGRATION METHODS

182

The second step is compute the Fourier coef¬cients gn by a Fast Fourier Transform (FFT).

In a Galerkin representation, the boundary value problem Eq.(9.24) is

gn

’ fn =

(1 + n2 )fn = gn (9.28)

1 + n2

where the {fn } are the spectral coef¬cients of F (u, x, t) and where we have the used the

fact that the second derivative of exp(inx) is ’n2 exp(inx). An inverse FFT then gives the

grid point values of F .

Fourth order Runge-Kutta requires four evaluations of F per time step. One might

think that this need for four boundary value solutions per time step would make “implicitly-

implicit” equations rather expensive to solve. The irony is that both our examples were

originally proposed as cheaper alternatives to other PDEs that could be solved by fully ex-

plicit methods!