COMPUTATIONAL FLUID DYNAMICS - Anderson, Dale A. Computational Fluid Mechanics and Heat Transfer100-150
COMPUTATIONAL FLUID DYNAMICS - Anderson, Dale A. Computational Fluid Mechanics and Heat Transfer100-150
COMPUTATIONAL FLUID DYNAMICS - Anderson, Dale A. Computational Fluid Mechanics and Heat Transfer100-150
3.18 Consider a steady-state conduction problem governed by Laplace's equation with convective
boundary conditions (see Fig. P3.2). The formal statement of the boundary condition is
-k a T / aY)bdy =h(T - Too),which can be readily cast into finite-difference form as -k[ (T
o
-
Too)j ,l y) +O(~y) =h(T
o
- Too). Use the control-volume approach to develop an expression for the
boundary condition at point O. Evaluate the T.E. in this expression assuming that Laplace's equation
applies at the boundary point.
3.19 Consider a heat conduction problem governed by a T / a t =a (a
2
T / a x
2
). Develop a finite-
difference representation for this equation by the control-volume approach. Do not assume that the
grid is uniform.
3.20 For 2-D steady-state conduction in a solid, apply the control-volume method to derive an
appropriate difference expression for the boundary temperature in control volume Bin Fig. 3.7for
a di a ba t i c wa l l boundary conditions.
3.21 Solve the 1-D heat equation using forward-time centered-space differences with a (,l t j ,l x
2
)
=t. Let the grid consist of fivepoints, including three interior and two boundary points. Assume a
constant unity wall temperature and a zero initial temperature on the interior. Complete this
calculation for 10 integration steps. Compare your results with those obtained in the example of
Section 3.3.4.
3.22 Refer to Fig. 3.10. Following the methodology illustrated in the text material associated with
Fig. 3.10, develop an appropriate finite-volume expression for the heat flux across the boundary from
c to d.
Too' h
2 0
x
Figure P3.2
BASICS OF DISCRETIZATION METHODS 99
3.23 Refer to Fig. 3.10and the associated text material. In an example in Section 3.5, an expression
was developed for the heat flux across boundary a -b for control volume B (Eq. 3.1(0). Following a
similar methodology, develop an appropriate finite-volume expression for the heat flowinto volume
B across the boundary from d to a . If the difference scheme is to be conservative, this heat flow
should be equal in magnitude but opposite in direction to the inflowcomputed for volume A across
boundary d-a . Check to see if this is true.
3.24 Show that the amplification factor derived for the [mite-difference solution of the heat
equation, Eq. (3.101), could be obtained by direct substitution of a solution of the form
In this form em represent the Fourier coefficients of the initial error distribution and g m is the
amplification factor. Identify g m with Eq. (3.106). Discuss the convergence of the solution and relate
your conclusions to the Lax equivalence theorem.
3.25 Use avon Neumann stability analysis to showfor the waveequation that asimple explicit Euler
predictor using central differencing in space is unstable. The difference equation is
"+1 " f1t(U7 +
I
-u7 -
1
)
u, = u: -c-
I I f1x 2
Nowshow that the same difference method is stable when written as the implicit formula
f1t ( u" +I - u~ +I )
u~+1=u" - c- 1+1 I-I
I I f1x 2
3.26 The DuFort-Frankel method for solving the heat equation requires solution of the difference
equation
-r = : ' a
= --(u~ - u"+1 - u"-I +u" )
2f1t (f1X)2 I +I I I I-I
Develop the stability requirements necessary for the solution of this equation.
3.27 Prove that the CFL condition is the stability requirement when the Lax-Wendroff method is
applied to solvethe simple 1-Dwave equation. The difference equation is of the form
c M c2(~d
urI =u[ - --(u7+1 - u7-1) +---2 (U[+I - 2uJ +UJ -I)
2 f1x 2(~x)
3.28 An implicit scheme for solving the heat equation is given by
a M [1 2 ]
U"+I =u" +-- _(u"+1 - 2u~+1+U"+I) +_(un - Zu" +u" )
I I (f1X)2 3 1+1 I I-I 3 1+1 I I-I
Apply the Fourier stability analysis to this scheme and determine the stability restrictions, if any.
3.29 An implicit scheme for solvingthe first-order wave equation is given by
c f1t
u"+1 =u~ - __ (U"+I - U~+I)
I I f1x 1+I I
Apply the Fourier stability analysis to this scheme and determine the stability restrictions, if any.
3.30 The leap frog method for solvingthe 1-Dwaveequation is given by
urI-urI UJ +I-uJ -I
+c =0
2f1t 2f1x
Apply the Fourier stability analysis to this method, and determine the stability restrictions, if any.
100 FUNDAMENTALS
3.31 Determine the stability requirement necessary to solve the 1-D heat equation with a source
term
a u a
2
u
- =a -- + ku
a t a x
2
Use the central-space, forward-time difference method. Does the von Neumann necessary condition,
Eq. (3.121), make physical sense for this type of computational problem?
3.32 Use the matrix method to determine the stability of the Lax method used to solve the
first-order wave equation on a mesh with two interior points and two boundary points. Assume the
boundaries are held at constant values u
1eft
=1, uright =O.
3.33 Use the matrix method and evaluate the stability of the numerical method used in Prob. 3.21
for the heat equation using a five-point mesh. Howmany frequencies must one be concerned with in
this case?
3.34 In attempting to solve a simple PDE, a system of finite-difference equations of the form
ur I =[A)uJ has evolved, where
[ A)= [ l ~vV l i v l~J
Investigate the stability of this scheme.
3.35 The application of a finite-difference scheme to the heat equation on a three-point grid results
in the following systemof equations:
uri =[A)uJ
where
[A) =[~
o
1- 2r
o
and r =a ,l t /(,~x)2. Determine the stability of this scheme.
3.36 The upstream scheme
U)' ! + 1 - u 7 u ' : - u r : 1
-'----'- +c) )- =0
,It ,l x
is used to solvethe wave equation on a four-point grid for the boundary conditions
and the initial conditions (n =1)
ui =1 u~ =u~ =ul =0
Use the matrix method to determine the stability restrictions for this method.
CHAPTER
FOUR
APPLICATION OF NUMERICAL METHODS
TO SELECTED MODEL EQUATIONS
In this chapter weexamine in detail various numerical schemes that can be used
to solve simple model partial differential equations (PDEs). These model
equations include the first-order wave equation, the heat equation, Laplace's
equation, and Burgers' equation. These equations are called model equations
because they can be used to "model" the behavior of more complicated PDEs.
For example, the heat equation can serve as a model equation for other
parabolic PDEs such as the boundary-layer equations. All of the present model
equations have exact solutions for certain boundary and initial conditions. We
can use this knowledge to quickly evaluate and compare numerical methods that
we might wish to apply to more complicated PDEs. The various methods
discussed inthis chapter were selected because they illustrate the basic properties
of numerical algorithms. Each of the methods exhibits certain distinctive features
that are characteristic of a class of methods. Some of these features may not be
desirable, but the method is included anyway for pedagogical reasons. Other
very useful methods have been omitted because they are similar to those that
are included. Space does not permit a discussion of all possible methods that
could be used.
1 01
102 FUNDAMENTALS
4.1 WAVE EQUATI ON
The one-dimensional (I-D) wave equation is a second-order hyperbolic PDE
given by
a
2
u a
2
u
-- = c
2
--
a t
2
a x
2
This equation governs the propagation of sound waves traveling at awavespeed
c in a uniform medium. A first-order equation that has properties similar to
those of Eq. (4.1) is given by
(4.1)
au au
- + c- = 0
a t a x
Note that Eq. (4.1) can be obtained from Eq. (4.2). We will use Eq. (4.2) as our
model equation and refer to it as the first-order I-D wave equation, or simply
the "wave equation." This linear hyperbolic equation describes a wave
propagating inthe x direction with avelocity c, and it can be used to "model" in
arudimentary fashion the nonlinear equations governing inviscidflow. Although
we will refer to Eq. (4.2) as the wave equation, the reader is cautioned to be
aware of the fact that Eq. (4.1) isthe classical waveequation. More appropriately,
Eq. (4.2) is often called the I-D linear convection equation.
The exact solution of the wave equation [Eq. (4.2)] for the pure initial value
problem with initial data
c > 0
(4.2)
u(x,O) =F(x) -oo<x<oo (4.3)
is given by
u(x, t) =F(x - ct)
(4.4)
Let us now examine some schemes that could be used to solve the wave
equation.
4.1.1 Euler Explicit Methods
The following simple explicit one-step methods,
u
n
+
1
- u" u" - u,n
J J +C J +1 =
A.! dx
c> o
(4.5)
Uj+1 - Uj-1
+c-'-----'-- =
2dx
(4.6)
have truncation errors (T.E.s) of O[A.!, dx] and O[A.!, (dX)2], respectively. We
refer to these schemes as being first-order accurate, since the lowest-order term
in the T.E. isfirst order, i.e., A.! and dx for Eq. (4.5) and A.! for Eq. (4.6). These
schemes are explicit, since only one unknown uj+ 1 appears in each equation.
Unfortunately, when the von Neumann stability analysis is applied to these
schemes, we find that they are unconditionally unstable. These simple schemes,
APPLICATION OF NUMERICAL METHODS TO SELECTED MODEL EQUATIONS 103
therefore, prove to be worthless in solving the wave equation. Let us now
proceed to look at methods that have more utility.
4.1.2 Upstream (First-Order Upwind. or Windward) Differencing Method
The simple Euler method, Eq. (4.5), can be made stable by replacing the forward
space difference by abackward space difference, provided that the wavespeed c
is positive. If the wave speed is negative, a forward difference must be used to
assure stability. This point is discussed further at the end of the present section.
For a positive wave speed, the following algorithm results:
Uj -Uj _l
+ c dx
=0 c> o
(4.7)
This is afirst-order accurate method with T.E. of O[M, dx] . The von Neumann
stability analysis shows that this method is stable, provided that
(4.8)
where v =c dt /dx.
Let us substitute Taylor-series expansions into Eq. (4.7) for uj + l and Uj -l .
The following equation results:
1([ (dd (dd ])
M uj +Mu, +-2-Ut t +-6-Ut t t +... - uj
+ ~(un - [ un - dxu +(dX)2 U - (dX)3 U +... ]) =0 (4.9)
dx ] ] x 2 xx 6 xxx
Equation (4.9) simplifies to
Note that the left-hand side of this equation corresponds to the wave equation
and the right-hand side is the T.E., which is generally not zero. The significance
of terms in the T.E. can be more easily interpreted if the time-derivative terms
are replaced by spatial derivatives. In order to replace u., by aspatial-derivative
term, we take the partial derivative of Eq. (4.10) with respect to time, to obtain
(4.11)
104 FUNDAMENTALS
and take the partial derivative of Eq. (4.10) with respect to x and multiply by
-c:
c t xt c
2
AX c(At)2 C
2
(AX)2
-cu
t x
- c
2
u
xx
=-2-
Ut t x
- -2-
uxXX
+--6-
Ut t t x
+ 6 u
xxxx
+...
(4.12)
Adding Eqs. (4.11) and (4.12) gives
In a similar manner, we can obtain the following expressions for Ut t l ' U
t t X
' and
U
xxt
:
U
t t t
=-c
3
u
xxx
+O[ At , AX]
U
t t x
=C2~xx +O[ At , AX]
u
xxt
=-cu
xxx
+O[ At , Ax]
(4.14)
Combining Eqs. (4.10), (4.13), and (4.14) leaves
2
C Ax c(AX) 2
u , +CU
x
=-2-(1 - v)u
xx
- 6 (2v - 3v +l )u
xxx
+ O[ (AX)3,(Ax)2 At ,AX(Ad,(Ad] (4.15)
An equation, such as Eq. (4.15), is called a modi f i ed equa t i on (Warming and
Hyett, 1974). It can be thought of as the PDE that is actually solved (if sufficient
boundary conditions were available) when a finite-difference method is applied
to a PDE. It is important to emphasize that the equation obtained after
substitution of the Taylor-series expansions, i.e., Eq. (4.10), must be used to
eliminate the higher-order time derivatives rather than the original PDE, Eq.
(4.2). This is due to the fact that a solution of the original PDE does not in
general satisfy the difference equation, and since the modified equation
represents the difference equation, it is obvious that the original PDE should
not be used to eliminate the time derivatives.
The process of eliminating time derivatives can be greatly simplified if a
table is constructed (Table 4.1). The coefficients of each term in Eq. (4.10) are
placed in the first row of the table. Note that all terms have been moved to the
left-hand side of the equation. The u., term is then eliminated by multiplying
Eq. (4.10) by the operator
At a
2 a t
o o o o
o o o o o
o o o
o o
o o
o o o
o
o
<l
~ : a
~ ~
w w
,......1('1) _ 1 " ' 1 "
I +
~ 0-
<l ~
_N, ) ...,~IC'" \ ~
---:--- I ...-..
<l ~ <l
<l ~
: a - I ' " ! t
~I~ +
'-" '
o
o
o
o
o
o
o
o
o
1 05
106 FUNDAMENTALS
and adding the result to the first row, i.e., Eq. (4.10). This introduces the term
-(c At /2)u
t x
' which is eliminated by multiplying Eq. (4.10) by the operator
cAt a
2 a x
and adding the result to the first two rows of the table. This procedure is
continued until the desired time derivatives are eliminated. Each coefficient in
the modified equation is then obtained by simply adding the coefficients in the
corresponding column of the table. The algebra required to derive the modified
equation can be programmed on a digital computer using an algebraic
manipulation code.
The right-hand side of the modified equation [Eq. (4.15)] isthe T.E., since it
represents the difference between the original PDE and the finite-difference
approximation to it. Consequently, the lowest order term on the right-hand side
of the modified equation gives the order of the method. In the present case, the
method is first-order accurate, since the lowest order term is O[ At , AX] . If
v =1, the right-hand side of the modified equation becomes zero, and the wave
equation is -solved exactly. In this case, the upstream differencing scheme
reduces to
which is equivalent to solving the wave equation exactly using the method of
characteristics. Finite-difference algorithms that exhibit this behavior are said to
satisfy the shi f t condi t i on (K utler and Lomax, 1971).
The lowest order term of the T.E. in the present case contains the partial
derivative U
w
which makes this term similar to the viscous term in 1-D fluid
flowequations. For example, the viscous term in the 1-DNavier-Stokes equation
(see Chapter 5) may be written as
a 4
a x (Tx) =" 3 J Luxx (4.16)
if aconstant coefficient of viscosity is assumed. Thus, when v = 1= 1, the upstream
differencing scheme introduces an a r t i f i ci a l vi scosi t y into the solution. This is
often called implicit artificial viscosity, as opposed to explicit artificial viscosity,
which is purposely added to a difference scheme. Artificial viscosity tends to
reduce all gradients in the solution whether physically correct or numerically
induced. This effect, which is the direct result of even derivative terms in the
T.E., is called di ssi p a t i on.
Another quasi-physical effect of numerical schemes iscalled di sp er si on. This
isthe direct result of the odd derivative terms that appear in the T.E. As aresult
of dispersion, phase relations between various waves are distorted. The combined
effect of dissipation and dispersion is sometimes referred to as di f f usi on.
Diffusion tends to spread out sharp dividing lines that may appear in the
computational region. Figure 4.1 illustrates the effects of dissipation and
dispersion on the computation of a discontinuity. In general, if the lowest order
APPLICATION OF NUMERICAL METIIODS TO SELECTED MODEL EQUATIONS 107
(a)
(b)
(c)
Figure 4.1 Effects of dissipation and dispersion. (a) Exact solution, (b) Numerical solution distorted
primarily by dissipation errors (typical of first-order methods), (c) Numerical solution distorted
primarily by dispersion errors (typical of second-order methods),
term in the T.E. contains an even derivative, the resulting solution will
predominantly exhibit dissipative errors. On the other hand, if the leading term
is an odd derivative, the resulting solution will predominantly exhibit dispersive
errors.
In Chapter 3wediscussed atechnique for finding the relative errors in both
amplitude (dissipation) and phase (dispersion) from the amplification factor. At
this point it seems natural to ask if the amplification factor is related to the
modified equation. The answer isdefinitely yes!Warming and Hyett (1974) have
developed a "heuristic" stability theory based on the even derivative terms in
the modified equation and have determined the phase shift error by examining
the odd derivative terms. However, the analysis of Warming and Hyett has been
shownby Chang (1987) to be restricted to schemes involvingonly two time levels
(n, n +1). Before showing the correspondence between the modified equation
and the amplification factor, let us first examine the amplification factor of the
present upstream differencing scheme:
G =(1 - v +v cos f3) - i( v sin f3)
The modulus of this amplification factor,
IGI =[(1- v +v cos f3)2 +(- v sin f3)2f/2
isplotted in Fig. 4.2for several values of v, It isclear fromthis plot that v must
be less than or equal to 1if the von Neumann stability condition IGI ~ 1isto be
met.
The amplification factor, Eq. (4.17), can also be expressed in the exponential
formfor a complex number:
(4.17)
G =IGle
i t l
)
where cp is the phase angle given by
cp=tan-1[Im(G)]=tan-1( -vsinf3 )
Re( G) 1- v +v cos {3
108 FUNDAMENTALS
UNIT CIRCLE
1 . 50 1 . 00
0. 50
0. 00 0. 50
1 . 00
IGI
Figure 4.2 Amplification factor modulus for upstream differencing scheme.
The phase angle for the exact solution of the wave equation (< /1e) is determined
in asimilar manner once the amplification factor of the exact solution is known.
In order to find the exact amplification factor we substitute the elemental
solution
into the wave equation and find that a =-i kmc, which gives
The exact amplification factor is then
u(t +Lit)
G= ---
e u(t)
eikmrx -c(t +dl)l
eikm(x-cI)
which reduces to
where
and
Thus the total dissipation (amplitude) error that accrues from applying the
upstream differencing method to the wave equation for N steps is given by
APPLICATION OF NUMERICAL METIIODS TO SELECTED MODEL EQUATIONS 109
where Ao is the initial amplitude of the wave. Likewise, the total dispersion
(phase) error can be expressed as N(< P e - < p ). The relative phase shift error
after one time step is given by
tan-
1
[( - v sin f3)/(l - v +v cos /3)]
- /3v
(4.18)
and is plotted in Fig. 4.3for several values of v, For small wave numbers (i.e.,
small /3) the relative phase error reduces to
< P 1 2 ) 2
- :::;1 - -(2v - 3v +1 {3
< P e 6
(4.19)
If the relative phase error exceeds 1for a given value of /3, the corresponding
Fourier component of the numerical solution has awave speed greater than the
exact solution, and this is a l ea di ng p ha se er r or . If the relative phase error is less
than 1, the wave speed of the numerical solution is less than the exact wave
speed, and this is a l a g g i ng p ha se er r or . The upstream differencing scheme has a
leading phase error for 0.5 <v <1and a lagging phase error for v <0.5.
Exa mp l e 4.1 Suppose the upstream differencing scheme is used to solve the
wave equation (c =0.75) with the initial condition
U(X,O) =sin(67 TX)
O~ x~ l
and periodic boundary conditions. Determine the amplitude and phase errors
after 10steps if At =0.02and Ax =0.02.
Sol ut i on In this problem a unique value of /3can be determined because the
exact solution of the wave equation (for the present initial condition) is
UNIT CIRCLE
1 . 50 1 . 00
0. 50
0. 00 0. 50
1 . 00
F" J gUre 4.3 Relative phase error of upstream differencing scheme.
110 FUNDAMENTALS
represented by asingleterm of aFourier series. Sincethe amplification factor is
also determined using a single term of a Fourier series that satisfies the wave
equation, the frequency of the exact solution is identical to the frequency
associated with the amplification factor, i.e., t ; =k
m
/27r. Thus the wave
number for the present problem is given by
2m7r 67r
k =-- =- =67r
m 2L 1
and f3 can be calculated as
f3 =; Ax =(67r)(0.02) =0.127r
Using the Courant number,
v =-- =----- =0.75
Ax (0.02)
cI1t (0.75)(0.02)
the modulus of the amplification factor becomes
[
2 2]1/2
IGI = (1 - v +v cos (3) +(- v sin (3) =0.986745
and the resulting amplitude error after 10steps is
(1 - IGIN)A
o
=(1 - IGl
lO
)(1) =1- 0.8751 =0.1249
The phase angle (cf after one step,
cf>=tan-1[ -vsinf3 ]=-0.28359
1- v +v cos f3
can be compared with the exact phase angle (cf > e) after one step,
cf > e = - f3v = - 0.28274
to givethe phase error after 10steps:
10(cf >e - cf =0.0084465
Let us now compare the exact and numerical solutions after 10steps where the
time is
t =10I1t =0.2
The exact solution is given by
u(x,0.2) =sin[67r(x - 0.15)]
and the numerical solution that results after applying the upstream differencing
scheme for 10steps is
u(x,0.2) =(0.8751) sin[67r(x - 0.15) - 0.0084465]
I n order to show the correspondence between the amplification factor and
the modified equation, we write the modified equation [Eq. (4.15)] in the
APPLICATION OF NUMERICAL METIIODS TO SELECTED MODEL EQUATIONS 111
following form:
(4.20)
where C
Zn
and C
2n
+ 1 represent the coefficients of the even and odd spatial-
derivative terms, respectively. Warming and Hyett (974) have shown that a
necessary condition for stability is
(4.21)
where C
ZI
represents the coefficient of the lowest order even derivative term.
This is analogous to the requirement that the coefficient of viscosity in viscous
flow equations be greater than zero. In Eq. (4.15) the coefficient of the lowest
order even derivative term is
cAx
C
z
=-2-0 - v)
(4.22)
and therefore the stability condition becomes
cAx
--O-vO
2
(4.23)
or v <1, which was obtained earlier from the amplification factor. It should be
remembered that the "heuristic" stability analysis, i.e., Eq. (4.21), can only
provide a necessary condition for stability. Thus, for some finite-difference
algorithms, only partial information about the complete stability bound is
obtained, and for others (such as algorithms for the heat equation) a more
complete theory must be employed.
Warming and Hyett have also shown that the relative phase error for
difference schemes applied to the wave equation is given by
< P 1 ;.. n Zn
-: i : =1- - 1... (-1)(k
m
) C
2n
+
1
o/e C n= l
(4.24)
where k
m
=f31 Ax is the wave number. For small wave numbers, we need only
retain the lowest order term. For the upstream differencing scheme, wefind that
< P 1 ( f3 )Z 1 Z Z
- ~ 1- - (-1) - C =1- - (2v - 3 v +1)f3
< P e c Ax 3 6
(4.25a )
whichisidentical to Eq. (4.19). Thus wehave demonstrated that the amplification
factor and the modified equation are directly related.
The upstream method given by Eq. (4.7) may be written in a more general
form to account for either positive or negative wave speeds. The method is
112 FUNDAMENTALS
normally written separately for these two cases as
At
U,!+l =un - C-(U
n
- un )
J J Ax J J-l
At
U
n
+ 1=U'! - C-(U
n
- Un)
J J Ax J +1 J
However, if we make use of the following definitions,
c' =t(c +[cl)
c-= t(c -leD
the upstream scheme may be written as the single expression
c> o
c< o
At
u
n
+
1
=u" - -[c+(u
n
- u" ) +c-(u
n
- un)]
J J Ax J i : 1 J + 1 J
Upon substituting for the values of c" and c", the final form becomes
At [c]At
n+ 1_ n ( n n) + (n 2 n + n ) (425b)
u
j
-u
j
-c
2AX
uj+1-U
j
-
1
2Ax u
j
+
1
- u
j
u
j
_
1
.
It is interesting to note that this form of the upstream scheme gives the
impression that it is a centered method. We recognize the first difference term
as a central-difference approximation and interpret the last term as an artificial
viscosity term. The function of this last term isto add the appropriate dissipation
to produce the upstream scheme when c is either positive or negative.
4.1.3 Lax Method
The Euler method, Eq. (4.6), can be made stable by replacing u' / with the
averaged term (U,/+l +u'/_1)/2. The resulting algorithm is the well-known Lax
method (Lax, 1954), which was presented earlier:
n+ 1 (n + n )/2 n n
u
j
- u
j
+
1
u
j
_
1
u
j
+
1
- u
j
-
1
......:....---'-----'--- +c =0
At 2Ax
(4.26)
This explicit one-step scheme is first-order accurate with T.E. of O[At, (AX)2/
At] and is stable if Ivi ~ 1. The modified equation is given by
2
CAX(1 ) c(Ax)
U +cu =-- - - v U + (1- v
2
)u +...
t x 2 v xx 3 xxx
Note that this method is not uniformly consistent, since (Ax)2 / At may not
approach zero in the limit as At and Ax go to zero. However, if v is held
constant as At and Ax approach zero, the method is consistent. The Lax
method isknown for its large dissipation error when v -=/ = 1.This large dissipation
is readily apparent when we compare the coefficient of the u
xx
term in Eq.
(4.27) with the same coefficient in the modified equation of the upstream
(4.27)
APPLICATION OF NUMERICAL METIIODS TO SELECTED MODEL EQUATIONS 113
differencing scheme for various values of II. The large dissipation can also be
observed in the amplification factor
G =cos {3- i II sin (3
(4.28)
which is described in Section 3.6.1. The modulus of the amplification factor is
plotted in Fig. 4.4(a). The relative phase error is given by
cf > tan- 1( - II tan (3)
cf > e - (311
which produces a leading phase error, as seen in Fig. 4.4(b).
4.1.4Euler Implicit Method
The algorithms discussed previously for the waveequation have all been explicit.
The following implicit scheme,
u
n
+
1
- u" C
} } +__ (U,!+l - u
n
+
1
) =0
At 2 Ax }+l }-l
is first-order accurate with T.E. of O[M,(AX)2] and, according to a Fourier
stability analysis, is unconditionally stable for all time steps. However, a system
of algebraic equations must be solved at each new time level. To illustrate this,
let us rewrite Eq. (4.29) so that the unknowns at time level (n +1) appear on
the left-hand side of the equation and the known quantity uj appears on the
right-hand side. This gives
(4.29)
(4.30)
v =1.0
0.75
0.5
o . 5
1.00 0.00 1.00
1& 1
2.00 1.00
0.00
1.00
(a) (b)
Figure 4.4 Lax method. (a) Amplification factor modulus. (b) Relative phase error.
114 FUNDAMENTALS
or
(4.31)
where a
j
=v/2, d
j
=1, b
j
=- v/2, and C, =u' j . Consider the computational
mesh shown in Fig. 4.5, which contains M +2grid points in the x direction and
known initial conditions at n =O. Along the left boundary, uZ+1 has a fixed
value of uo. Along the right boundary, u~/+\ can be computed as part of the
solution using characteristic theory. For example, if v =1, then u~: \ =u~.
Applying Eq. (4.31) to the grid shown in Fig. 4.5, we find that the following
system of M linear algebraic equations must be solved at each (n +1) time
level:
b d a U
n+ 1
M-l M-l M-l M-l
[A]
o
I n Eq. (4.32), C
1
and C
M
are given by
C
1
=u~ - buZ+
1
C
M
=u~ - a u~++\
o
o
[u]
(4.32)
(4.33)
where uZ+ 1 and u~++\ are the known boundary conditions.
Matrix [A] in Eq. (4.32) is a tridiagonal matrix. A technique for rapidly
solving a tridiagonal system of linear algebraic equations is due to Thomas
(1949) and is called the Thomas algorithm. I n this algorithm, the system of
equations is first put into upper triangular form by replacing the diagonal
t~ :::1111111111111111
j'" 0
2
Figure 4.5 Computational mesh.
M M+ 1
APPLICATION OF NUMERICAL MElHODS TO SELECTED MODEL EQUATIONS 115
elements d, with
i =2,3, ... ,M
and the C; with
b;
C; - dC; -1 i =2,3, ... ,M
;-1
The unknowns are then computed using back substitution starting with
n+ 1 C
M
u
M
=-
d
M
and continuing with
j=M - 1, M - 2, ... ,1
Further details of the Thomas algorithm are given in Section 4.3.3.
In general, implicit schemes require more computation time per time step
but, of course, permit a larger time step, since they are usually unconditionally
stable. However, the solution may become meaningless if too large atime step is
taken. This is due to the fact that a large time step produces large T.E.s. The
modified equation for the Euler implicit scheme is
u
t
+cU
x
=(i c
2
M)u
xx
- [ i c(Ll x)2 +t c3(Mi ] uxxx +... (4.34)
which does not satisfy the shift condition. The amplification factor
1-i"sinf3
G = (4.35)
1+,,2 sin? f3
and the relative phase error
cP
C P e
tan- 1 ( - " sin f3)
- f3"
(4.36)
are plotted in Fig. 4.6. The Euler implicit scheme is very dissipative for
intermediate wave numbers and has a large lagging phase error for high wave
numbers.
\ I: LO b
o.~
I
1.00 0.00 1.00
cp/cpe
1.00
0.00
IGI
1.00
(a)
(b)
Figure 4.6 Euler implicit method. (a) Amplification factor modulus. (b) Relative phase error.
116 FUNDAMENTALS
4.1.5 Leap Fr og Method
The numerical schemes presented so far in this chapter for solving the linear
waveequation are all first-order accurate. I n most cases, first-order schemes are
not used to solve PDEs because of their inherent inaccuracy. The leap frog
method is the simplest second-order accurate method. When applied to the
first-order waveequation, this explicit one-step three-time-Ievel scheme becomes
U,!+I - un-I
] ]
2 a t
n n
u
j
+
1
- u
j
_
1
+c =0
2 a x
(4.37)
The leap frog method is referred to as a three-time-Ievel scheme, since U must
be known at time levels nand n - 1in order to find u at time level n +1. This
method has a T.E. of O[(a t)2,(a x)2] and is stable whenever 1,,1 ~ 1. The
modified equation is given by
4
c(a x) 4 2
120 (9" - 10" +l )uxxxxx +...
(4.38)
The leading term in the T.E. contains the odd derivative u
xxx
' and hence the
solution will predominantly exhibit dispersive errors. This is typical of second-
order accurate methods. I n this case, however, there are no even derivative
terms in the modified equation, so that the solution will not contain any
dissipation error. As a consequence, the leap frog algorithm is neutrally stable,
and errors caused by improper boundary conditions or computer round-off will
not be damped (assuming periodic boundary conditions and I " I ~ 1). The
amplification factor
G =(1 - ,,2sirr' (3)1/2 - i v sin f 3
and the relative phase error
tan-I [ - " sin f 3/ (1 - ,,2 sin
2
f 3 )1/2]
- f 3"
(4.39)
~e
(4.40)
are plotted in Fig. 4.7.
The leap frog method, whilebeing second-order accurate with no dissipation
error, does have its disadvantages. First, initial conditions must be specified at
two-time levels. This difficulty can be circumvented by using a two-time-Ievel
scheme for the first time step. A second disadvantage is due to the "leap frog"
nature of the differencing (i.e., u' r 1 does not depend on u'[), so that two
independent solutions develop as the calculation proceeds. And finally, the leap
frog method may require additional computer storage because it is a three-
time-level scheme. The required computer storage is reduced considerably if a
simple overwriting procedure isemployed, whereby -: 1 isoverwritten by u' r I.
APPLICATION OF NUMERICAL MElHODS TO SELECTED MODEL EQUATIONS 117
m
0. 5
1.00
0.00
IGI
1.00
1.00 0.00 1.00
~/~e
(a)
(b)
Figure 4.7 Leap frog method. (a) Amplification factor modulus. (b) Relative phase error.
4.1.6 Lax- Wendr of T Method
The Lax-Wendroff finite-difference scheme (Lax and Wendroff, 1960) can be
derived from a Taylor-series expansion in the following manner:
(4.41)
Using the wave equations
U
t
=-cu
x
Ut t =c
2
u
xx
(4.42)
Equation (4.41) may be written as
ur 1 =u'j - ca t u; +~c2(a t)2 U
xx
+o[(a d]
(4.43)
And finally, if u; and u
xx
are replaced by second-order accurate central-
difference expressions, the well-known Lax-Wendroff scheme is obtained:
This explicit one-step scheme is second-order accurate with a T.E. of
O[(a x)2,(a t)2] and is stable whenever Ivl ~ 1. The modified equation for this
method is
2
(a x) 2
u,+CU
x
=-c--(1 - v )u
6 xxx
3
c(a x) 2
8 v(1 - v )u
xxxx
+... (4.45)
The amplification factor
G =1 - v
2
(1 - cos (3) - i v sin f3
(4.46)
118FUNDAMENTALS
1.00
0. 00
IGI
1.00
(a)
1.00
0. 00
, I ' e
(b)
1.00
Figure 4.8 Lax-Wendroff method. (a) Amplification factor modulus. (b) Relative phase error.
and the relative phase error
tan-I {- " sin /3/[ 1 - ,,2(1 - cos /3)]}
- /3"
(4.47)
are plotted in Fig. 4.8. The Lax-Wendroff scheme has a predominantly lagging
phase error except for large wave numbers with .f03 <" <1.
4.1.7 Two- Step Lax- Wendr of T Method
For nonlinear equations such as the inviscidflowequations, atwo-step variation
of the original Lax-Wendroff method can be used. When applied to the wave
equation, this explicit two-step three-time-level method becomes
Step 1:
n+ 1/2 _ (n + n)/2 n _ n
Uj+
I
/2 Uj+1 U
j
+C Uj+1 Uj =0
tJ .112 tJ .X
(4.48)
Step 2:
U
n
+ 1/2 _ ul " 1/2
j+I/2 j-I/2
+c tJ .X
=0
(4.49)
This scheme is second-order accurate with a T.E. of O[(tJ .X)2,(tJ .t)2] and is
stable whenever I" I ~ 1. Step 1isthe Lax method applied at the midpoint j+t
for a half time step, and step 2is the leap frog method for the remaining half
time step. When applied to the linear waveequation, the two-step Lax-Wendroff
scheme is equivalent to the original Lax-Wendroff scheme. This can be readily
shown by substituting Eq. (4.48) into Eq. (4.49). Since the two schemes are
equivalent, it follows that the modified equation and the amplification factor are
the same for the two methods.
APPLICATION OF NUMERICAL MElHODS TO SELECTED MODEL EQUATIONS 119
4.1.8 MacCor mack Method
The MacCormack method (MacCormack, 1969) is a widely used scheme for
solving fluid flow equations. It is a variation of the two-step Lax-Wendroff
scheme that removes the necessity of computing unknowns at the grid points
j+~and j- ~.Because of this feature, the MacCormack method isparticularly
useful when solvingnonlinear PDEs, as is shown in Section 4.4.3. When applied
to the linear wave equation, this explicit, predictor-corrector method becomes
Corrector:
(4.50) Predictor:
(4.51)
The term U 'j +1 is a temporary "predicted" value of U at the time level
n +1. The corrector equation provides the final value of U at the time level
n +1. Note that in the predictor equation a forward difference is used for
a Uf a x, while in the corrector equation a backward difference is used. This
differencing can be reversed, and in some problems it is advantageous to do so.
This is particularly true for problems involving moving discontinuities. For the
present linear wave equation, the MacCormack scheme is equivalent to the
original Lax-Wendroff scheme. Hence the truncation error, stability limit,
modified equation, and amplification factor are identical with those of the
Lax-Wendroff scheme.
4.1.9 Second- Or der Upwind Method
The second-order upwind method (Warming and Beam, 1975) is avariation of
the MacCormack method, which uses backward (upwind) differences in both the
predictor and corrector steps for c > 0:
Predictor: (4.52)
Corrector:
1[ _ ca t (_ _)
u
n
+1 =_ un +un+
1
- __ u
n
+1 _ un+1 _
J 2 J J a x J J -I
ca t ]
a x (u'j - 2u'j_1 +U'j-2)
(4.53)
The addition of the second backward difference in Eq. (4.53) makes this scheme
second-order accurate with T.E. of O[(a t)2, (a tXa x), (a x)2]. If Eq. (4.52) is
substituted into Eq. (4.53), the following one-step algorithm is obtained:
n+1 _ n (n n) +1 ( 1)( n 2 n + n )
u
j
- u
j
- " u
j
- U
j
_
1
2"" - U
j
- U
j
_
1
U
j
-
2
(4.54)
120 FUNDAMENTALS
The modified equation for this scheme is
2 4
c(a x) (a x) 2
U
t
+CU
x
= (1 - v)(2 - v)u - --v(1 - v) (2 - v)u +..,
6 xxx 8a t XXXX
(4.55)
The second-order upwind method satisfies the shift condition for both v =1
and v =2. The amplification factor is
G =1- 2v ( v +2(1 - v) sirr' ~) sin? ~ - iv sin J3(1+2(1 - v) sirr' ~)
(4.56)
and the resulting stability condition becomes 0 ~ v ~ 2. The modulus of the
amplification factor and the relative phase error are plotted in Fig. 4.9. The
second-order upwind method has a predominantly leading phase error for
o <v <1and a predominantly lagging phase error for 1 <v <2. We observe
that the second-order upwind method and the Lax-Wendroff method have
opposite phase errors for 0 <v <1. This suggests that aconsiderable reduction
in dispersive error would occur if alinear combination of the two methods were
used. Fromm's method of zero-average phase error (Fromm, 1968) is based on
this observation.
4.1.10 Time- Center ed I mplicit Method (Tr apezoidal
Dif f er encing Method)
A second-order accurate implicit scheme can be obtained if the twoTaylor-series
expansions
n (a t)2 n (a!)3 n
ur i =u' j +a t (ut )j +-2-(U
t t
)j +-6-(U
t t t
)j +...
161
v = 1.25 and 0.75
2. 00 and 1. 00
1.50 and 0.50
1. 75 and 0.25
1.00
0,00
1.00
(a)
(b)
Figure 4.9 Second-order upwind method. (a) Amplification factor modulus. (b) Relative phase error.
APPLICATION OF NUMERICAL MElHODS TO SELECTED MODEL EQUATIONS 121
are subtracted and (utt)r I is replaced with
n+Inn
(Utt)j =(Utt)j +a t(Uttt)j +...
The resulting expression becomes
a t
urI =u'j +T[(ut)n +(ut)n+t +O[(ad]
(4.58)
The time differencing in this equation is known as trapezoidal differencing or
Crank-Nicolson differencing. Upon substituting the linear wave equation u, =
- cux' we obtain
ca t [ ]
urI =u'j - -2- (u)n +(u)n+1 j+O[(ad] (4.59)
And finally, if the U
x
terms are replaced by second-order central differences, the
time-centered implicit method results:
(4.60)
This method has second-order accuracy with T.E. of O[(a x)2,(a t)2] and is
unconditionally stable for all time steps. However, a tridiagonal matrix must be
solved at each new time level. The modified equation for this scheme is
_ [c
3
(a t)2 +c(a x)2]u
ut +CUx = 12 6 x x x
_ [c(a x)4 +c
3
(a t)2(a x)2
120 24
c
4
(a t )4]
+ 80 Uxxxxx +... (4.61)
Note that the modified equation contains no even derivative terms, so that the
scheme has no implicit artificial viscosity. When this scheme is applied to the
nonlinear fluid dynamic equations, it often becomes necessary to add some
explicit artificial viscosity to prevent the solution from "blowing up." The
addition of explicit artificial viscosity (i.e., "smoothing" term) to this scheme will
be discussed in Section 4.4.7. The modulus of the amplification factor,
1 - (;v/2) sin f3
G = (4.62)
1+(;v/2) sin f3
and the relative phase error are plotted in Fig. 4.10.
The time-centered implicit method can be made fourth-order accurate in
space if the difference approximation given by Eq. (3.31) is used for u
x
:
1 a x 4
(U)j =2a x 1 +8; /6
Uj
+o[(a x)]
(4.63)
The modified equation and phase error diagram for the resulting scheme can be
found in the work by Beam and Warming (1976).
122 FUNDAMENTALS
1 . 00
0. 00
IGI
1 . 00 1 . 00
0. 00
, I ' e
(a)
(b)
Figure 4.10 Time-centered implicit method. (a) Amplification factor modulus. (b) Relative phase
error.
4.1.11 Rusanov (Bnrstein-Mirlnl Method
The methods presented thus far for solvingthe wave equation have either been
first-order or second-order accurate. Only asmall number of third-order methods
have appeared in the literature. Rusanov (1970) and Burstein and Mirin (1970)
simultaneously developed the following explicit three-step method:
St 1
(I) - 1.( n + n) _ 1. (n _ n)
ep . U
j
+
I
/
2
- 2 U
j
+
1
U
j
3" U
j
+
1
U
j
Step 2:
(2) _ n _ 1(I) _ (I) )
U
j
- U
j
3" U
j
+
I
/2 U
j
-I/2
1
n+ 1 _ n (2 n +7 n 7 n +2 n )
U
j
- U
j
- 24" - U
j
+
2
U
j
+
1
- U
j
_
1
U
j
_
2
(4.64)
Step 3:
3 (2 2)
- -" U{ ) - U{ )
8 }+I }-I
w
- 24(U'j+2 - 4u'j+1 +6u'j - 4u'j_1 +U'j_2)
Step 3contains the fourth-order difference term
~4 n n 4 n +6 n 4 n + n
U
x
u
j
=u
j
+
2
- u
j
+
1
u
j
- u
j
_
1
U
j
-2
which is multiplied by a free parameter w. This term has been added to make
the scheme stable. The need for this term is apparent when we examine the
stability requirements for the scheme:
I " I ~ 1
4,,2- ,,4 ~ W ~ 3 (4.65)
If the fourth-order difference term were not present (i.e., w =0), we could not
satisfy Eq. (4.65)for 0 <" ~ 1. The modified equation for this method is
c(a x)3 ( w )
u
t
+cU
x
=- 24 --;; - 4" +,,3 U
xxxx
c(a x)4
+ (-5w +4+15,,2- 4,,4)U +... (4.66)
120 xxxxx
APPLICATION OF NUMERICAL ME1HODS TO SELECTED MODEL EQUATIONS 123
In order to reduce the dissipation of this scheme, wecan make the coefficient of
the fourth derivative equal to zero by letting
(4.67)
In alike manner, wecan reduce the dispersive error by setting the coefficient of
the fifth derivative to zero, which gives
w=
(4.68)
The amplification factor for this method is
G =1 - ~ sirr' f3 - 2w sin" f3 - i u sin f3 [1 +~(1 - ,,2) sin? f3] (4.69)
2 3 2 3 2
The modulus of the amplification factor and the relative phase error are plotted
in Fig. 4.11. This figure shows that the Rusanov method has a leading or a
lagging phase error, depending on the value of the free parameter w.
4.1.12 War ming- Kutler - Lomax Method
Warming et al. (1973) developed a third-order method that uses MacCormack's
method for the first two steps and has the same third step as the Rusanov
method. This so-called WK L method is given by
Step 1:
Step 2:
v. 0. 50. W Z1 . 5
v = 1 . 00. W ' "3. 0
V ' " 0. 75. w" 2.5
1 . 00
0. 00
IGI
1 . 00 1 . 00
0. 00
, I ' e
1 . 00
(a)
(b)
Figure 4.11 Rusanov method. (a) Amplification factor modulus. (b) Relative phase error.
124 FUNDAMENTALS
Step 3:
1
n+ 1 - n (2 n +7 n 7 n +2 n )
u
j
- u
j
- 24 v - u
j
+
2
u
j
+
1
- u
j
_
1
u
j
-
2
(4.70a )
3
- - v(U< 2) - U(2) )
8 J +1 J -1
W
- 24(Uj+2 - 4Uj+1 +6uj - 4Uj_1 +Uj_2)
This method has the same stability bounds as the Rusanov method. In addition,
the modified equation is identical to Eq. (4.66) for the present linear wave
equation. The WK L method has the same advantage over the Rusanov method
that the MacCormack method has over the two-step Lax-Wendroff method.
4.1.13 Runge- Kutta Methods
Runge-K utta methods are frequently employed to solve ordinary differential
equations (ODEs). They can also be applied to solvePDEs (Lomax et aI., 1970;
and Jameson et aI., 1981, 1983). In fact, several of the methods described
previously in this section can be derived using Runge-K utta methodology. The
first step in this process is to convert the PDE into a "pseudo-ODE." This is
accomplished by separating out a partial derivative with respect to a single
independent variable in the marching direction and placing the remaining
partial derivatives into a term that is a function of the dependent variable. For
example, the linear wave equation can be written as
u, =R(u) (4.70b)
where R(u) =-cu
x
. This pseudo-ODE is a time-continuous equation, and any
integration scheme applicable to ODEs, including Runge-K utta methods, may
be used. Once the time differencing is completed, the partial derivatives
contained in R(u) can be differenced using appropriate spatial differences. To
illustrate this approach, let us apply the second-order Runge-K utta method, also
referred to as the improved Euler's method (Carnahan et aI., 1969), to Eq.
(4.70b), which gives
Step 1:
u(l) =u" +D.t R"
Step 2:
D.t
u
n
+
1
=u" +-(Rn +R(l ))
2
where
R" ==R(u
n
) =-cu~
APPLICATION OF NUMERICAL METHODS TO SELECTED MODEL EQUATIONS 125
The term R(l) in step 2 can be evaluated by making use of step 1 in the
following manner:
R(1) =- CU~l)
-c(u~ +LltR~)
-cu~ +c
2
Lltu~x
Substituting this expression for R(l) into step 2yields
Llt
u
n
+
1
=u" +-( + Lcu" . +c
2
Lltu
n
)
2 x xx
If second-order accurate central differences are then used to approximate the
spatial derivatives, the resulting scheme becomes
c Llt c
2
(Llt)2
n+ 1 _ n ( n n) + (n 2 n + n )
uj - uj - 2Llx uj+1 - uj-1 2(Llx)2 uj+1 - uj uj_1
which is the second-order accurate Lax-Wendroff scheme, Eq. (4.44).
Procedures and equations for obtaining nth-order Runge-K utta methods
can be found in the works by Carnahan et al. (1969), Luther (1966), and Yu et
al. (1992). A fourth-order Runge-K utta method, attributed to K utta, is given by
Llt
u(1) =u" +-Rn
2
Step 1:
Step 2:
Step 3:
Llt
U(2) =u" +-R(l)
2
U(3) =u" +Llt R(2)
Llt
u
n
+
1
=u" +-(Rn +2R(1) +2R(2) +R(3
6
Step 4:
where R() =- cu~) for the linear wave equation. If second-order accurate
spatial differences are inserted into this algorithm, the resulting scheme will
have aT.E. of O[(Llt)4, (Llx)2]. In order to obtain higher-order spatial accuracy,
it is convenient to employ compact difference schemes (Yu et aI., 1992) with the
Runge-K utta time stepping.
4.1.14 Additional Comments
The improved accuracy of higher-order methods is at the expense of added
computer time and additional complexity. These factors must be considered
carefully when choosing a scheme to solve a PDE. In general, second-order
accurate methods provide enough accuracy for most practical problems.
For the 1-Dlinear waveequation, the second-order accurate explicit schemes
such as the Lax-Wendroff scheme give excellent results with a minimum of
computational effort. An implicit scheme may not be the optimum choice in this
126 FUNDAMENTALS
case because the solution is unsteady and intermediate results are typically
desired at relatively small time intervals.
4.2 HEAT EQUATI ON
The 1-D heat equation (diffusion equation),
a u a
2
u
a t =a a x2
is a parabolic PDE. In its present form, it is the governing equation for heat
conduction or diffusion in a 1-Disotropic medium. It can be used to "model" in
arudimentary fashion the parabolic boundary-layer equations. The exact solution
of the heat equation for the initial condition
(4.71)
u(x,O) =f(x)
and boundary conditions
u(O, t) =u(1, t) =
is
00
u(x, t) = LAne-a k2t sin(kx)
n=l
(4.72)
where
An =2f f (x) sin(kx) dx
o
and k =n7T. Let us now examine some of the more important finite-difference
algorithms that can be used to solve the heat equation.
4.2.1 Simple Explicit Method
The following explicit one-step method,
Un+
1
- u" u" 2u
n
+u"
] ] j + 1 - j j-l
= a
~t (~ xt
is first-order accurate with T.E. of O[~t, (~t )2 J. At steady-state the accuracy is
O[ (~X)2J . As we have shown earlier, this scheme is stable whenever
(4.73)
(4.74)
where
a ~t
r=--
(~X)2
(4.75)
APPLICATION OF NUMERICAL METHODS TO SELECTED MODEL EQUATIONS 127
The modified equation is given by
[
1 2 a C: 1x)2]
u, - a uxx = -"2
a
I1t + 12 u
xxxx
[
1 2 1 2 1 4]
+ -a
3
(l 1t ) - -a
2
I1t (l 1x) +-a (l 1x) U +...
3 12 360 xxxxxx
(4.76)
We note that if r =-t, the T.E. becomes of O[ (l 1t )2,(l 1x)4] . It is also interesting
to note that no odd derivative terms appear in the T.E. As a consequence, this
scheme, as well as almost all other schemes for the heat equation, has no
dispersive error. This fact can also be ascertained by examining the amplification
factor for this scheme:
G =1 +2r(cos f3 - 1) (4.77)
which has no imaginary part and hence no phase shift. The amplification factor
is plotted in Fig. 4.12 for two values of r and is compared with the exact
amplification factor of the solution. The exact amplification (decay) factor is
obtained by substituting the elemental solution
G
~ S I MP l E E X P L I C I T
- E X AC T r=}
-1.00L- ~ ~ ~ ~------~--~~~~----~
0.00 0.50 1.50
B
3.00
Figure 4.12 Amplification factor for simple explicit method.
128 FUNDAMENTALS
into
G =
e
u(t +a t)
u(t)
which gives
(4.78)
or
(4.79)
where f3 =k
m
a x. Hence the amplitude of the exact solution decreases by the
factor e-
r
{ 32 during one time step, assuming no boundary condition influence.
In Fig. 4.12, we observe that the simple explicit method is highly dissipative
for large values of f3 when r =4-. As expected, the amplification factor is in
closer agreement with the exact decay factor when r =i .
The present simple explicit scheme marches the solution outward from the
initial data line in much the same manner as the explicit schemes of the
previous section. This is illustrated in Fig. 4.13. In this figure we see that the
unknown u can be calculated at point P without any knowledge of the boundary
conditions along AB and CD. We know, however, that point P should depend on
the boundary conditions along AB and CD, since the parabolic heat equation
has the characteristic t =const. From this we conclude that the present explicit
scheme (with a finite a t) does not properly model the physical behavior of the
parabolic PDE. It would appear that an implicit method would be the more
appropriate method for solving a parabolic PDE, since an implicit method
normally assimilates information from all grid points located on or below the
characteristic t =const. On the other hand, explicit schemes seem to provide a
t
:::
~
~
~
~
;.
~
:::
t.t
~
~
/~
~
;. t:.x
~D
~
~~~ L
~
~
~
"
::;.
::::
T T
~
~
I'"
~I
J
'\~
" '\
~
;..
~
U,
I
..II
~
~
1 J
~L
::;.
;;.
I'
t" \
E : c
-,.. -...
-~
x
B
A
I NI TI AL DATA
LI NE
Figure 4.13 Zone of influence of simple explicit scheme.
APPLICATION OF NUMERICAL MElHODS TO SELECfED MODEL EQUATIONS 129
more natural finite-difference approximation for hyperbolic PDEs that possess
limited zones of influence.
Exa mp l e 4.2 Suppose the simple explicit method is used to solve the heat
equation (a =0.05) with the initial condition
u(x,O) =sin(27 TX) O~x~l
and periodic boundary conditions. Determine the amplitude error after 10steps
if a t =0.1 and a x =0.1.
Solution A unique value of f3 can be determined in this problem for the same
reason that was given in Example 4.1. Thus the value of f3 becomes
f3 =k
m
a x =(27 T )(0.1) =0.27T
After computing r ,
a a t (0.05)(0.1)
r =-- =----::-- =0.5
(a x)2 (0.1)2
the amplification factor for the simple explicit method is given by
G =1+2r(cos f3 - 1) =0.809017
while the exact amplification factor becomes
G
e
=e-
r
{ 32 =0.820869
As a result, the amplitude error is
AoIG;o - GlO l =(1)(0.1389 - 0.1201) =0.0188
Using Eq. (4.72), the exact solution after 10steps (t =1.0) is given by
u(x,1) =e-
a 41T2
sin(27 TX) =0.1389sin(27 TX)
which can be compared to the numerical solution:
u(x,1) =0.1201sin(27 TX)
4.2.2Richardson's Method
Richardson (1910) proposed the following explicit one-step three-time-level
scheme for solving the heat equation:
u
n
+
1
_ ur : '
] ]
2 a t
Uj+l - 2uj +Uj-l
= a
(a x)2
(4.80)
This scheme is second-order accurate with T.E. of O[(a t)2, (a x)2]. Unfor-
tunately, this method proves to be unconditionally unstable and cannot be used
to solve the heat equation. It is presented here for historic purposes only.
130 FUNDAMENTALS
4.2.3 Simple Implicit (Laasonen) Method
A simple implicit scheme for the heat equation was proposed by Laasonen
(1949). The algorithm for this scheme is
U
n
+
1
_ u" u
n
+
1
_ 2u
n
+
1
+u
n
+
1
J J =a J + 1 J J -1
a t (a x)2
If we make use of the central-difference operator
s;,2 n n 2 n + n
u
x
u
j
= u
j
+
1
- u
j
u
j
_
1
(4.81)
we can rewrite Eq. (4.81) in the simpler form:
(4.82)
This scheme has first-order accuracy with a T.E. of O[ a t ,(a x)2] and is
unconditionally stable. Upon examining Eq. (4.82), it is apparent that a
tridiagonal system of linear algebraic equations must be solved at each time
level n +1.
The modified equation for this scheme is
[
1 a (a x)2]
u, - a uxx = " 2a
2
a t + 12 uxxxx
[
1 3()2 1 2 ( 2 1 4]
+ "3
a
a t + 12a a t a x) + 360a (a x) uxxxxxx +...
(4.83)
It is interesting to observe that in this modified equation, the terms in the
coefficient of u
xxxx
are of the same sign, whereas they are of opposite sign in
the modified equation for the simple explicit scheme, Eq. (4.76). This observation
can explain why the simple explicit scheme is generally more accurate than the
simple implicit scheme when used within the appropriate stability limits. The
amplification factor for the simple implicit scheme,
G =[1 +2r(1 - cos {3)r
1
(4.84)
is plotted in Fig. 4.14for r =t and is compared with the exact decay factor.
4.2.4 Crank-Nicolson Method
Crank and Nicolson (1947) used the following implicit algorithm to solve the
heat equation:
(4.85)
APPLICATION OF NUMERICAL MElHODS TO SELECTED MODEL EQUATIONS 131
r = i
{
- o- SI MPLE I MPLI CI T
- - CRANK- NI COLSON
~ DUFORT- FRANKEL
- EX ACT
G
0. 50
Figur e 4.14 Amplification factors for several methods.
This unconditionally stable algorithm has become very well known and is
referred to as the Crank-Nicolson scheme. This scheme makes use of trapezoidal
differencing to achieve second-order accuracy with a T.E. of O[ (a t )2,(a x)2] .
Once again, atridiagonal systemof linear algebraic equations must be solved at
each time level n +1. The modified equation for the Crank-Nicolson method is
a (a x)2 [1 3 2 1 4]
u, - a uxx = 12 uxxxx + 12a (a t ) +360a (a x) uxxxxxx +...
(4.86)
The amplification factor
1 - r(1 - cos f3)
G=--------, -
1+r(1 - cos f3)
(4.87)
is plotted in Fig. 4.14for r =t.