Figure 13.1: Splitting by dimension: a two-dimension process such as diffusion, which

causes a drop of food coloring to widen in all directions simultaneously (top), is treated as

two one-dimensional processes, happening one-after-the-other instead of simultaneously

(bottom).

One strategy is to use sophisticated matrix methods. For for constant coef¬cient and

separable BVPs, there are special fast spectral algorithms. For non-separable BVPs, pre-

conditioned matrix iterations can be very effective.

The alternative is to attack and simplify before the spatial discretization is applied: The

“splitting” or “fractional step” family of algorithms.

De¬nition 30 (SPLITTING/FRACTIONAL STEPS) An implicit time-marching algorithm which

replaces simultaneous processes by sequential steps to increase ef¬ciency is called a SPLITTING

or FRACTIONAL STEPS method. The split may be by dimensions: three-dimensional diffusion,

for example, may be split into three one-dimensional substeps of diffusion in the x-direction only fol-

lowed by diffusion in y and then z. The splitting may also be by physics: advection on one fractional

step, pressure adjustment (waves) on another, and diffusion/viscosity on a third.

The key idea of splitting, simultaneous-into-sequential, can be illustrated by the two-

dimensional diffusion equation, ut = uxx + uyy . If the matrices which are the spatial dis-

cretizations of the x and y derivatives are denoted by Λ1 and Λ2 , then the exact solution

CHAPTER 13. SPLITTING & ITS COUSINS

254

is

(13.1)

u(„ ) = exp([Λ1 + Λ2 ] „ ) u(0)

Suppose, however, that we solve the pair of equations that results from assuming that the

solution diffuses in the x-direction ¬rst and then in the y-direction (Fig. 13.1):

vt = vxx ; v(t = 0) = u(0)

(13.2)

wt = wyy ; w(t = 0) = v(„ )

The discretized solution is

(13.3)

w(„ ) = exp(Λ1 „ ) exp(Λ2 „ ) u(0)

In general, it is not automatic that these two solutions are the same because for matrices

(or linear operators), the usual exponential identity must be generalized to, for arbitrary

matrices A, B,

exp(A + B) = exp(A) exp(B) exp( ’ [A, B]/2) Weyl formula (13.4)

where

[A, B] ≡ A B ’ B A Commutator of A, B (13.5)

However, if both A and B are proportional to the time step „ , then the commutator must

be O(„ 2 ). It follows that in the limit „ ’ 0, the splitting system (13.2) provides at least a

¬rst-order-in-time approximation to the diffusion equation, regardless of how the time and

space derivatives are discretized. Furthermore, we can replace the matrices for diffusion

by any other pair of physical processes and the argument still applies.

The great advantage of splitting for the diffusion problem is that it requires only one-

dimensional boundary value problems. Even though we solve many such problems, one

BVP in x for each interior grid point in y, for example, this is still much cheaper than

solving a single two-dimensional BVP.

For the hydrodynamic equations, the advantage of splitting-by-process is that advec-

tion can be treated by a different algorithm than pressure adjustment, which in turn can be

different from diffusion. We need not even use the same spatial discretization; we can use

¬nite differences for diffusion and spectral methods for other parts (as in Lambiotte et al.,

1982, for example).

Because of these virtues, the splitting method, also known as the method of fractional

steps, was extensively developed in the Former Soviet Union as described in books by

Yanenko (1971) and Marchuk (1974). Strang (1968) independently invented this family of

algorithms, which are often called “Strang splitting” in the U. S. where these ideas caught

on more slowly. Although especially popular with spectral methods, splitting can be used

with ¬nite difference or ¬nite element spatial discretizations just as well, and indeed was

developed by Yanenko and Marchuk only for non-spectral methods.

Splitting sounds almost too wonderful to be true. Not quite, but there are snags. The

biggest is boundary conditions. Viscosity furnishes the highest derivatives for the Navier-

Stokes equations, so we cannot impose the usual rigid boundary conditions on the non-

viscous fractional steps. This can and does lead to growth of error at the boundaries. (This

problem does not arise in Fourier or spherical harmonics applications where there are no

boundaries.) However, at the price of sacri¬cing the simplicity of straightforward splitting,

some good remedies have been developed. Splitting is still widely used.

13.2. FRACTIONAL STEPS FOR DIFFUSION 255

13.2 The Method of Fractional Steps for the Diffusion Equa-

tion

The great vice of the Crank-Nicholson method is that it is cheap only in one-dimension.

In two or three dimensions, even with a ¬nite difference method, we must invert a large

matrix at each time step. With direct matrix-solving methods (LU factorization, also known

as Gaussian elimination), the cost per time step in two dimensions is O(N 4 ) and O(N 6 )

operations for ¬nite difference and spectral methods, respectively, where N is the number

of grid points in either x or y. Ouch! These costs can be greatly reduced by clever iterative

matrix-solvers as described later, but are still formidable.

The three-dimensional diffusion equation is

(13.6)

ut = uxx + uyy + uzz + g(x, y, z)

where g(x, y, z) is a known forcing function.

The method of fractional steps avoids inverting full matrices by advancing the solution

in one coordinate at a time. First, write

(13.7)

Λ = Λ 1 + Λ2 + Λ 3

where Λ1 = ‚xx , Λ2 = ‚yy , and Λ3 = ‚zz . Then the algorithm is

un+1/4 ’ un 1

Λ1 un+1/4 + un (13.8)

=

„ 2

’ un+1/4

n+1/2

u 1

Λ2 un+1/2 + un+1/4 (13.9)

=

„ 2

’ un+1/2

n+3/4

u 1

Λ3 un+3/4 + un+1/2 (13.10)

=

„ 2

’ un+3/4

n+1

u

(13.11)

=g

„

Each fractional step is the Crank-Nicholson method for a one-dimensional diffusion equa-

tion. Consequently, if Λ1 is generated by applying the usual three-point difference formula

to ‚xx , it is only necessary to invert tridiagonal matrices at each time step. As shown in Ap-

pendix B, a tridiagonal matrix may be inverted in O(N ) operations. The operation count

for the fractional steps method is directly proportional to the total number of grid points

even in two or three dimensions; it lacks the O(N 2 ) penalty of the unsplit Crank-Nicholson

method (versus an explicit scheme).

The pseudospectral method, alas, generates non-sparse matrices ” but it is far better to

invert 3N 2 matrices each of dimension N than a single matrix of dimension N 3 since the

expense of inverting a matrix scales as the cube of its dimension. Later, we shall discuss

cost-reducing iterations. An alternative direct method, faster but restricted to constant co-

ef¬cient differential operators, is to use Galerkin™s method. For this special case of constant

coef¬cients, the Galerkin matrix is banded and cheap to invert as explained in the Chapter

15.

It is obvious that (13.8)“(13.11) is cheap; it is not so obvious that it is accurate. However,

eliminating the intermediate steps shows that the method of fractional steps is equivalent

to the scheme

un+1 = M (’„ )’1 M („ ) un + „ M (’„ )’1 g (13.12)

CHAPTER 13. SPLITTING & ITS COUSINS

256

in whole time steps where the operator M is given by

„2 „3

„

M („ ) ≡ 1 + Λ + {Λ1 Λ2 + Λ1 Λ3 + Λ2 Λ3 } + Λ1 Λ2 Λ3 (13.13)

2 4 8

Eq. (13.12) is identical with the corresponding formula for the Crank-Nicholson scheme

except for the terms in „ 2 and „ 3 in the de¬nition of M („ ). Because of cancellations, the

overall error in both the splitting method and Crank-Nicholson is O(„ 3 ) for u(x, y, z, t).

At ¬rst glance, the entire premise of the method of fractional steps seems ridiculous.

Physically, diffusion is not a sequential process. However, if we consider only a small time

interval of length „ , then it is a good mathematical approximation to split the time interval

into three, and pretend that there is only diffusion in x on the ¬rst fractional step, only

diffusion in y on the second, and only diffusion in z on the third. (Note that there are no

factors of (1/3) in (13.8)“(13.10): we allow each second derivative to diffuse over the whole

time interval even with the splitting-up.) Our ¬ction introduces extra error; the accuracy

for a given time step „ is poorer than that of the Crank-Nicholson method. However, it

is the same order-of-magnitude ” O(„ 3 ). The cost-per-time-step of splitting is orders-of-

magnitude smaller than for the Crank-Nicholson scheme.

13.3 Pitfalls in Splitting, I: Boundary Conditions

Fractional step methods introduce some new dif¬culties peculiar to themselves. Variables

at fractional time levels are not solutions of the full set of original equations. What bound-

ary conditions should we impose on these intermediates?

For behavioral boundary conditions, this dif¬culty disappears because behavioral con-

ditions are never explicitly imposed. Instead, one merely expands all the unknowns, in-

cluding intermediates at fractional time levels, in Fourier series (for spatial periodicity) or

rational Chebyshev functions (on an in¬nite interval).

Numerical boundary conditions require careful analysis, however. If the domain is rect-

angular and the boundary values are independent of time, then the treatment of boundary

conditions is straightforward. For simplicity, consider the two-dimensional heat equation.

With splitting, one must solve a two-point boundary value problem in x at each value

of y on the interior of the grid (Fig. 13.2). The values of u(0, y) and u(1, y) give us just

what we need, however. Similarly, the known values of u at the top and bottom of the

rectangle give the boundary conditions for the second splitting step: solving a set of two-

point boundary value problems in y, one for each value of x. The matrix representation

of the x-derivatives is a large matrix, but it decouples into independent subproblems as

illustrated in Fig. 13.3.

When the boundary conditions vary with time, we have troubles. The intermediates

in splitting ” un+1/2 , for example ” do not correspond to the physical solution at t =

(n + 1/2)„ , but rather are merely intermediate stages in advancing the solution by a whole

time step.

When the domain is irregular, however, the proper boundary conditions require real

art. Spectral elements, which subdivide a complicated domain into many subdomains

only slightly deformed from rectangles, greatly alleviates this problem. Therefore, we re-

fer the curious reader to Yanenko (1971), who explains how to apply splitting with ¬nite

differences to a problem with curvy boundaries.

13.3. PITFALLS IN SPLITTING, I: BOUNDARY CONDITIONS 257

Figure 13.2: In a split scheme, it is necessary to solve one-dimensional boundary value

problems on vertical and horizontal lines. One such (horizontal) line is marked by the four

circles which are the unknowns. Instead of a single 20x20 matrix, splitting requires (on this

6x7 grid) the solution of ¬ve 4x4 problems and four 5x5 problems.

Figure 13.3: Matrix of the second x-derivative: Dots denote nonzero matrix elements. There

is one block for each interior grid point in y. Each block is the matrix of a one-dimensional

boundary value problem in x.

CHAPTER 13. SPLITTING & ITS COUSINS

258

13.4 Pitfalls in Splitting, II: Consistency

The consistency problem is that stability does not necessarily guarantee accuracy. General

truth: absolutely stable methods gain their stability from being inaccurate for the rapidly-

varying components of the solution. For example, the Crank-Nicholson method achieves

stability by arti¬cially slowing down the extremely fast decay of the high wavenumber

terms in the solution. This is acceptable if the only goal is the asymptotic, steady-state

solution, but we must always ask the engineering question: do we need to correctly model

the transient, too?

Even when only the steady state solution matters, we must still be careful because it is

possible that the inaccuracies of the absolutely stable method will poison the asymptotic

solution as well. Mathematicians describe this as a “lack of consistency”. Like stability,

consistency with the original differential equations may be either conditional ” we get the

correct answer only if „ is suf¬ciently small ” or absolute. The Crank-Nicholoson scheme

is “absolutely consistent” because it will correctly reproduce the asymptotic solution as

t ’ ∞ even if „ is very large.

The three-dimensional splitting method, alas, is only conditionally consistent with the

heat equation. To show this, write

(13.14)

u = T (x, y, z, t) + v(x, y, z)

where T (x, y, z, t) is the “transient” part of the solution and v(x, y, z) is the asymptotic

solution. It follows that

lim T (x, y, z, t) ≡ 0 (13.15)

t’0

2

(13.16)

v+g =0

Similarly, the steady solution of the splitting equations will satisfy (from (13.12))

[M (’„ ) ’ M („ )] v = „ g (13.17)

Recalling that

„2 „3

„

M („ ) ≡ 1 + {Λ1 Λ2 + Λ1 Λ3 + Λ2 Λ3 } + Λ1 Λ2 Λ3 (13.18)

Λ+

2 4 8

we see that only the odd powers of „ survive the subtraction in (13.17) to give (after divid-

ing out a common factor of „ )

„2

’Λ + (13.19)

Λ 1 Λ2 Λ3 v = g

4

Oops! When „ is suf¬ciently small, we can neglect the second term in the {} in (13.19)