. 54
( 136 .)


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

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


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)
wt = wyy ; w(t = 0) = v(„ )

The discretized solution is

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)


[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 The Method of Fractional Steps for the Diffusion Equa-

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

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

Λ = Λ 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
u 1
Λ2 un+1/2 + un+1/4 (13.9)
„ 2
’ un+1/2
u 1
Λ3 un+3/4 + un+1/2 (13.10)
„ 2
’ un+3/4

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

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.

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.

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

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)

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

’Λ + (13.19)
Λ 1 Λ2 Λ3 v = g

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


. 54
( 136 .)