. 39
( 136 .)



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.

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.


{a n} {u(xj)}
n=0,...,N j=0,...,N
Spectral Grid Points
Coefficients Values


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-

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:

ut + u ux + uxxx = 0

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

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

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

Table 9.3: A Selected Bibliography of Time-Marching with Multi-Dimensional Fourier
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)
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)
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.,

L ut = G(u, x, t)

To restate this in the canonical form

ut = F (u, x, t),

we must solve a boundary value problem for F :


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)

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
’ 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!


. 39
( 136 .)