CA AA214 Ch6

AA214: NUMERICAL METHODS FOR COMPRESSIBLE FLOWS 1 / 75
AA214: NUMERICAL METHODS FOR

COMPRESSIBLE FLOWS
The Finite Difference Method
These slides are based on the recommended textbook: Culbert B. Laney. “Computational
Gas Dynamics,” CAMBRIDGE UNIVERSITY PRESS, ISBN 0-521-62558-0
1 / 75
Outline
1 Conservative Finite Difference Methods in One Dimension
2 Forward, Backward, and Central Time Methods
3 Domain of Dependence and CFL Condition
4 Linear Stability Analysis
5 Formal, Global, and Local Order of Accuracy
6 Upwind Schemes in One Dimension
7 Nonlinear Stability Analysis
8 Multidimensional Extensions
2 / 75
Note: The material covered in this chapter equally applies to scalar

conservation laws and to the Euler equations, in one and multiple
dimensions. In order to keep things as simple as possible, it is presented
in most cases for scalar conservation laws: first in one dimension, then in
multiple dimensions. It is also presented in two dimensions for the Euler
equations.
3 / 75
Conservative Finite Difference Methods in One Dimension
Recall that scalar conservation laws are simple scalar models of the Euler equations that can
be written in strong conservation form as
∂u ∂f (u)
+ =0 (1)
∂t ∂x
Suppose that a 1D space is divided into grid points xi and “cells” [xi−1/2 , xi+1/2 ], where
xi+1/2 is called a “cell edge”
Also suppose that time is divided into time-intervals [t n , t n+1 ]

The conservation form of a finite difference method applied to the numerical solution of
equation (1) is defined as follows
\
n
∂u n n
∆t = −λ(fî+1/2 − fî−1/2 ) (2)
∂t i
where the subscript i designates the point xi , the superscript n designates the time t n , a
“hat” designates a time-approximation, and
∆t n+1 n
λ= , ∆t = t −t , ∆x = xi+1/2 − xi−1/2
∆x
4 / 75
One interpretation of the finite difference approach (2) and the

conservation form label is the approximation of the following integral
form of equation (1)
Z x !
1 i+1/2
n+1 n
[u(x, t ) dx − u(x, t )] dx
∆x xi−1/2
Z t n+1 !
∆t 1
= − [f u(xi+1/2 , t) − f u(xi−1/2 , t) ] dt
∆x ∆t tn
which clearly describes a conservation law (hint: recall the mean

value theorem)
5 / 75
Not every finite difference method can be written in conservation

form: Those which can are called conservative and their associated
quantities fˆn
i+1/2are called conservative numerical fluxes
finite difference methods derived from the conservation form of the
Euler equations or scalar conservation laws tend to be conservative
finite difference methods derived from other differential forms (for
example, primitive or characteristic forms) of the aforementioned
equations tend not to be conservative
conservative finite differencing implies correct shock and contact
locations
6 / 75
Like many approximation methods, conservative finite difference

methods can be divided into implicit and explicit methods
in a typical implicit method
\n \
∂u ∂u n n n+1 n+1
= (ui−K 1
, ..., ui+K 2
; ui−L1
, ..., ui+L 2
)
∂t i ∂t
fî+1/2
n
= fˆ(ui−K
n n n+1
+1 , ..., ui+K ; ui−L +1 , ..., ui+L )
1 2 1
n+1
2
(3)
so that from (2) one has
uin+1 = u(ui−K
n
1
n
, ..., ui+K2
n+1
; ui−L 1
, ..., uin+1 , ..., ui+L
n+1
2
) (4)
=⇒ the solution of a system of equations is required at each

time-step
7 / 75

in a typical implicit method
\n \
∂u ∂u n n n+1 n+1
= (ui−K 1
, ..., ui+K 2
; ui−L1
, ..., ui+L 2
)
∂t i ∂t
fî+1/2
n
= fˆ(ui−K
n n n+1
+1 , ..., ui+K ; ui−L +1 , ..., ui+L )
1 2 1
n+1
2
(3)
uin+1 = u(ui−K
n
1
n
, ..., ui+K2
n+1
; ui−L 1
, ..., uin+1 , ..., ui+L
n+1
2
) (4)
=⇒ the solution of a system of equations is required at each

time-step
n n n n
Note: if ui−K 1 +1
(ui−L 1 +1
) in (3) were written as ui−K 1
(ui−L 1
), one
would get the less convenient notation
uin+1 = u(ui−K
n
1 −1
n
, ..., ui+K 2
n+1
; ui−L 1 −1
n+1
, ..., ui+L 2
) instead of (4)
7 / 75

(continue)
in a typical explicit method
\n \
∂u ∂u n n
= (ui−K 1
, ..., ui+K2
; uin+1 )
∂t i ∂t
fî+1/2
n
= fˆ(ui−K
n n
+1 , ..., ui+K )
1 2
uin+1 = u(ui−K
n
1
n
, ..., ui+K 2
)
=⇒ only function evaluations are incurred at each time-step

n n n+1 n+1
(ui−K1
, ..., ui+K 2
)
and (ui−L 1
, ..., ui+L 2
) are called the stencil or direct
numerical domain of dependence of uin+1
K1 + K2 + 1 and L1 + L2 + 1 are called the stencil widths
8 / 75
Summary: typical stencil diagram
9 / 75
Like any proper numerical approximation, proper finite difference

approximation becomes perfect in the limit ∆x → 0 and ∆t → 0
an approximate equation is consistent if it equals the true equation in the limit
∆x → 0 and ∆t → 0
a solution to an approximate equation is convergent if it equals the true solution of
the true equation in the limit ∆x → 0 and ∆t → 0
Hence, a conservative (finite difference) approximation is consistent when
fˆ(u, ..., u) = f (u)
=⇒ in this case, the conservative numerical flux fˆ is said to be consistent

with the physical flux
A conservative numerical method — and therefore a conservative finite
difference method — automatically locates shocks correctly (however, it
does not necessarily reproduce their shapes correctly)
Methods that explicitly enforce the Rankine-Hugoniot relations are called
shock-tracking or shock-fitting methods
Methods that do not explicitly enforce the Ranking-Hugoniot relations are
called shock-capturing methods: They must be conservative and are the
subject of this course
10 / 75
Forward, Backward, and Central Time Methods
Forward Time Methods
Forward Time (FT) conservative finite difference methods

correspond to the choices
\n
∂u
∆t = uin+1 − uin and fî+1/2
n
= fˆ(ui−K
n
1 +1
n
, ..., ui+K )
∂t i 2
∂u
with Forward Space (FS) approximation of the term (xi , t n ), this
∂x
leads to the FTFS scheme
uin+1 = uin − λ(fî+1/2
n
− fî−1/2
n
), with fî+1/2
n n
= f (ui+1 )
∂u
with Backward Space (BS) approximation of the term (xi , t n ),
∂x
this leads to the FTBS scheme
n
− fî−1/2
n
), with fî+1/2
n
= f (uin )
∂u
with Central Space (CS) approximation of the term (xi , t n ), this
∂x
leads to the FTCS scheme
1
n
− fî−1/2
n
), with fî+1/2
n n
= (f (ui+1 ) + f (uin ))
2
11 / 75
Backward Time Methods
Backward Time (BT) conservative finite difference methods

correspond to the choices
\n
∂u
∆t = uin+1 − uin and fî+1/2
n
= fˆ(ui−K
n+1
1 +1
n+1
, ..., ui+K )
∂t i 2
∂u
∂x
leads to the BTFS scheme
n
− fî−1/2
n
), with fî+1/2
n n+1
= f (ui+1 )
∂u
∂x
this leads to the BTBS scheme
n
− fî−1/2
n
), with fî+1/2
n
= f (uin+1 )
∂u
∂x
leads to the BTCS scheme
1
n
− fî−1/2
n
), with fî+1/2
n
= n+1
f (ui+1 ) + f (uin+1 )
2
12 / 75
Central Time Methods
Central Time (CT) conservative finite difference methods correspond

to the choices
\n
∂u 1
∆t = (uin+1 − uin−1 ) and fî+1/2
n
= fˆ(ui−K
n
1 +1
n
, ..., ui+K )
∂t i 2 2
∂u
∂x
leads to the CTFS scheme
uin+1 = uin−1 − 2λ(fî+1/2
n
− fî−1/2
n
), with fî+1/2
n n
= f (ui+1 )
∂u
∂x
this leads to the CTBS scheme
uin+1 = uin−1 − 2λ(fî+1/2
n
− fî−1/2
n
), with fî+1/2
n
= f (uin )
∂u
∂x
leads to the CTCS scheme
1
uin+1 = uin−1 − 2λ(fî+1/2
n
− fî−1/2
n
), with fî+1/2
n n
= (f (ui+1 ) + f (uin ))
2
13 / 75
Domain of Dependence and CFL Condition
Numerical and Physical Domains of Dependence
Recall the theory of characteristics: A point in the x − t plane is

influenced only by points in a finite domain of dependence and
influences only points in a finite range of influence
Hence, the physical domain of dependence and physical range of

influence are bounded on the right and left by the waves with the
highest and lowest speeds
In a well-posed problem, the range of influence of the initial and
boundary conditions should exactly encompass the entire flow in the
x − t plane
14 / 75
The direct numerical domain of dependence of a finite difference

method is its stencil: For example, if the solution approximated by
an implicit finite difference method can be written as
uin+1 = u(ui−K
n
1
n
, ..., ui+K 2
n+1
; ui−L 1
n+1
, ..., ui+L 2
)
its direct numerical domain of dependence is the region of the x − t
n n n+1 n+1
plane covered by the points (ui−K 1
, ..., ui+K 2
; ui−L 1
, ..., ui+L2
)
Similarly, if the solution approximated by an explicit finite difference
method can be written as
uin+1 = u(ui−K
n
1
n
, ..., ui+K 2
)
its direct numerical domain of dependence is the region of the x − t
n n
plane covered by the points (ui−K 1
, ..., ui+K 2
)
The full (or complete) numerical domain of dependence of a finite
difference method consists of the union of its direct numerical
domain of dependence and the domain covered by the points of the
x − t plane upon which the numerical values in the direct numerical
domain of dependence depend upon
15 / 75
The Courant-Friedrichs-Lewy or (CFL) condition

The full numerical domain of dependence must contain the physical
domain of dependence
Any numerical method that violates the CFL condition misses

information affecting the exact solution and may blow up to infinity:
For this reason, the CFL condition is necessary but not sufficient for
numerical stability
16 / 75
Scalar Conservation Laws
Consider first the linear advection problem

∂u ∂u
+a = 0
∂t ∂x
1 if x < 0
u(x, 0) =
0 if x ≥ 0
Assume that a > 0: The exact solution is

1 if x − at < 0
u(x, t) = u(x − at, 0) =
0 if x − at ≥ 0
The FTFS approximation with ∆x = cst is
uin+1 = (1 + λa)uin − λaui+1
n

1 if i ≤ −1
ui0 = u(i∆x, 0) =
0 if i ≥ 0
∆t
where as before, λ =
∆x
17 / 75
Then

 1 if i ≤ −2
ui1 = 1 + λa if i = −1
0 if i ≥ 0



 1 if i ≤ −3
(1 + λa)(1 − λa) if i = −2

ui2 =
 (1 + λa)(1 + λa) if i = −1

0 if i ≥ 0

and so forth
The first two time-steps reveal that FTFS moves the jump in the
wrong direction (left rather than right!) and produces spurious
oscillations and overshoots
Furthermore, the exact solution yields u(0, ∆t) = 1, but FTFS yields
u01 = 0!
18 / 75
This is because for a > 0, FTFS violates the CFL condition
uin+1 = (1 + λa)uin − λaui+1

n
19 / 75
1
FTCS satisfies the CFL condition if λ ≤ ⇒ −1 ≤ λa ≤ 1
|a|
20 / 75
1
|a|
However, it almost always blow up (as will be seen in a homework):

This illustrates the fact that the CFL condition is a necessary but
not sufficient condition for numerical stability
20 / 75
1
|a|
However, it almost always blow up (as will be seen in a homework):

This illustrates the fact that the CFL condition is a necessary but
not sufficient condition for numerical stability
You can also check that when applied to the solution of any scalar
conservation law, the BTCS method always satifies the CFL
condition
20 / 75
For scalar conservation laws, the CFL condition translates into a simple inequality restricting
the time-step size
linear advection equation and explicit forward-time method with
uin+1 = u(ui−K
n n
, ..., ui+K )
1 2
in the x − t plane, the physical domain of dependence is the line of slope 1/a
in the x − t plane, the full numerical domain of dependence of uin+1 is bounded
∆t λ
on the left by a line of slope = and on the right by a line of slope
K1 ∆x K1
∆t λ
− =−
K2 ∆x K2
21 / 75
the time-step size
uin+1 = u(ui−K
n n
, ..., ui+K )
1 2
in the x − t plane, the physical domain of dependence is the line of slope 1/a
in the x − t plane, the full numerical domain of dependence of uin+1 is bounded
∆t λ
on the left by a line of slope = and on the right by a line of slope
K1 ∆x K1
∆t λ
− =−
K2 ∆x K2
hence, the CFL condition is

K2 K1
− ≤a≤ ⇔ −K2 ≤ λa ≤ K1 ⇔ −K2 ∆x ≤ a∆t ≤ K1 ∆x
λ λ
which requires that waves travel no more than K1 points to the right or K2
points to the left during a single time-step
21 / 75
the time-step size (continue)
uin+1 = u(ui−K
n n
, ..., ui+K ) (continue)
1 2
if K1 = K2 = K , the previous CFL condition becomes
∆x
λ|a| ≤ K ⇔ ∆t ≤ K (5)
|a|
for this reason, λa is called the CFL number or the Courant number
keep in mind however that in general, a = a(u) and therefore the CFL number
depends in general on the solution’s range
22 / 75
For scalar conservation laws, the CFL condition translates into a

simple inequality restricting the time-step size (continue)
linear advection equation and implicit backward-time method with
uin+1 = u(ui−K
n
1
n
, ..., ui+K2
n+1
; ui−L 1
n+1
, ..., ui+L 2
)
if L1 > 0 and L2 = 0, the full numerical domain of dependence of
uin+1 includes everything to the left of x = xi and beneath t = t n+1
in the x − t plane
if L1 = 0 and L2 > 0, the full numerical domain of dependence of
uin+1 includes everything to the right of x = xi and beneath t = t n+1
if L1 > 0 and L2 > 0, the full numerical domain of dependence of
uin+1 includes everything in the entire x − t plane beneath t = t n+1
23 / 75
For scalar conservation laws, the CFL condition translates into a

simple inequality restricting the time-step size (continue)
linear advection equation and implicit backward-time method with
uin+1 = u(ui−K
n
1
n
, ..., ui+K2
n+1
; ui−L 1
n+1
, ..., ui+L 2
)
if L1 > 0 and L2 = 0, the full numerical domain of dependence of
uin+1 includes everything to the left of x = xi and beneath t = t n+1
if L1 = 0 and L2 > 0, the full numerical domain of dependence of
uin+1 includes everything to the right of x = xi and beneath t = t n+1
if L1 > 0 and L2 > 0, the full numerical domain of dependence of
uin+1 includes everything in the entire x − t plane beneath t = t n+1
conclusion: as long as their stencil includes one point to
the left and one to the right, implicit methods avoid CFL
restrictions by using the entire computational domain
(hence, this includes BTCS but not BTFS or BTBS)
23 / 75
The Euler Equations
In 1D, the Euler equations have three families of waves that define
the physical domain of dependence
For each family of waves, a CFL condition can be established for a
given numerical method as in the case of a scalar conservation law:
Then, the overall CFL condition is the most restrictive of all
established CFL conditions
For example, if K1 = K2 = K , A is the Jacobian matrix of the
conservative flux vector, and ρ(A) denotes
its spectral radius
ρ(A) = max (|vx − a|, |vx |, |vx + a| , the CFL condition
of an
explicit forward-time method becomes recall (5)
∆x
λρ(A) ≤ K ⇔ ∆t ≤ K
ρ(A)
λρ(A) is called the CFL number or the Courant number
24 / 75
The Euler Equations
λρ(A) ≤ K
For supersonic flows, all waves travel in the same direction, either
left or right ⇒ the minimum stencil allowed by the CFL condition
n
contains either Wi−1 and Win for right-running supersonic flow, or
n n
Wi and Wi+1 for left-running supersonic flow
For subsonic flows, waves travel in both directions, and the
n
minimum stencil should always contain Wi−1 , Win , and Wi+1
n
25 / 75
The Euler Equations
λρ(A) ≤ K
For supersonic flows, all waves travel in the same direction, either
left or right ⇒ the minimum stencil allowed by the CFL condition
n
contains either Wi−1 and Win for right-running supersonic flow, or
n n
Wi and Wi+1 for left-running supersonic flow
For subsonic flows, waves travel in both directions, and the
n
minimum stencil should always contain Wi−1 , Win , and Wi+1
n
Hence, a smart or adaptive stencil can be useful for the case

of the Euler equations!
25 / 75
Linear Stability Analysis
Unstable solutions exhibit significant spurious oscillations and/or

overshoots
Unstable solutions of linear problems exhibit unbounded spurious
oscillations: Their errors grow to infinity as t → ∞
Hence the concept of instability discussed here for the solution of
linear problems is that of ubounded growth
26 / 75
Unstable solutions exhibit significant spurious oscillations and/or

overshoots
Unstable solutions of linear problems exhibit unbounded spurious
oscillations: Their errors grow to infinity as t → ∞
Hence the concept of instability discussed here for the solution of
linear problems is that of ubounded growth
Since unstable solutions typically oscillate, it makes sense
to describe the solution of a linear problem such as a linear
advection problem as a Fourier series (sum of oscillatory
trigonometric functions)
26 / 75
von Neumann Analysis
The Fourier series for the continuous (in space) solution u(x, t n ) on
any spatial domain [a, b] is
∞ X ∞
X x −a x −a
u(x, t n ) = a0n + n
am cos 2πm + n
bm sin 2πm
m=1
b−a m=1
b−a
(6)
27 / 75
The Fourier series for the continuous (in space) solution u(x, t n ) on
any spatial domain [a, b] is
∞ X ∞
X x −a x −a
u(x, t n ) = a0n + n
am cos 2πm + n
bm sin 2πm
m=1
b−a m=1
b−a
(6)
For a discrete solution u(xi , t n ), the Fourier series is obtained by

sampling (6) as follows
∞ X ∞
n n
X
n xi − a n xi − a
u(xi , t ) = a0 + am cos 2πm + bm sin 2πm
m=1
b−a m=1
b−a
(7)
Assume xi+1 − xi = ∆x = cst, x0 = a, and xN = b ⇒ xi − a = i∆x
and b − a = N∆x: This transforms (7) into
∞
X i i
u(xi , t n ) = a0n + n
am cos 2πm + bmn
sin 2πm (8)
m=1
N N
27 / 75
Recall that samples can only support wavelengths of 2∆x or longer

– the Nyquist sampling theorem states that samples spaced apart by
∆x perfectly represent functions whose shortest wavelengths are
4∆x (N/m ≥ 2): Hence (8) is truncated as follows
N/2
n n
X
n 2πmi n 2πmi
u(xi , t ) ≈ a0 + am cos + bm sin (9)
m=1
N N
which is called a discrete Fourier series
28 / 75
Recall that samples can only support wavelengths of 2∆x or longer

– the Nyquist sampling theorem states that samples spaced apart by
∆x perfectly represent functions whose shortest wavelengths are
4∆x (N/m ≥ 2): Hence (8) is truncated as follows
N/2
n n
X
n 2πmi n 2πmi
u(xi , t ) ≈ a0 + am cos + bm sin (9)
m=1
N N
which is called a discrete Fourier series
An equivalent expression in the complex plane using I as the
notation for the pure imaginary number (I 2 = −1) is
N/2 N/2
2πmi
X X
u(xi , t n ) ≈ Cmn e I N = uinm (10)
m=−N/2 m=−N/2
From (9), (10), and Euler’s formula e I θ = cos θ + I sin θ it follows

that
an − Ibm
n
an + Ibm n
C0n = a0n , Cmn = m , C−m n
= m
2 2
28 / 75
Hence, each term of the Fourier series can be written as

uinm = Cmn e (I ) = C n e I φm i
2πmi
N
m
2πm
where φm = and m = −N/2, · · · , N/2
N
Because of linearity, the amplification factor
Cmn+1
Gm = = Gm (λ)
Cmn
does not depend on n: However, it depends on λ (since uin+1 and uin

are produced by the numerical scheme being analyzed), which itself
depends on ∆t
Hence, each term of the Fourier series can be expressed as
Cmn Cm2 Cm1 0 I φm i
uinm = n−1 · · · C 1 C 0 Cm e = Gm · · · Gm Cm0 e I φm i = Gmn Cm0 e I φm i
Cm m m
n
where Gmn = Gmn (λ) = (Gm (λ))
29 / 75
Finally, assume that Cm0 = 1 (for example): This leads to
uinm = Gmn (λ)e I φm i
Conclusions
the linear approximation is linearly stable if |Gm (λ)| < 1 for all m
it is neutrally linearly stable if |Gm (λ)| ≤ 1 for all m and |Gm (λ)| = 1
for some m
it is linearly unstable if |Gm (λ)| > 1 for some m
Each of the above conclusion can be re-written in terms of
λ = ∆t/∆x
Application (in class): apply the von Neumann analysis to determine
the stability of the FTFS scheme for the linear advection equation
30 / 75
Matrix Method
Shortcomings of the von Neumann stability analysis method

requires the solution to be periodic (u0n = uNn )
requires constant spacing ∆x
does not account for the boundary conditions
31 / 75
Matrix Method
Shortcomings of the von Neumann stability analysis method

requires the solution to be periodic (u0n = uNn )
requires constant spacing ∆x
does not account for the boundary conditions
Alternative method: so-called Matrix (eigenvalue analysis) Method
based on the fact that for a linear problem and a linear
approximation method, one has
u n+1 = M(λ)u n , where u n = (uon u1n · · · uNn )T
and M is an amplification matrix which depends on the stencil of

the approximation scheme and on λ
This implies
u n = M n (λ)u 0
31 / 75
Matrix Method
u n = M n (λ)u 0
Suppose that M is diagonalizable
M(λ) = Q −1 Λ(λ)Q, Λ = diag (λ1 (λ), · · · , λN (λ))
Let v = Qu
Then
M n (λ) = Q −1 Λn Q ⇒ u n = Q −1 Λn (λ)Qu 0 ⇒ Qu n = Λn (λ)Qu 0
⇒ v n = Λn (λ)v 0
Conclusions
the linear approximation is linearly stable if ρ (M(λ)) < 1
it is neutrally linearly stable if ρ (M(λ)) = 1
it is linearly unstable if ρ (M(λ)) > 1
32 / 75
Matrix Method
Advantages of the Matrix Method for (linear) stability analysis

does not require the solution to be periodic
does not require constant grid spacing
incorporates the effects of the boundary conditions
Shortcoming: in general, the analytical computation of ρ (M(λ)) is not trivial
33 / 75
Matrix Method
Advantages of the Matrix Method for (linear) stability analysis

does not require the solution to be periodic
does not require constant grid spacing
incorporates the effects of the boundary conditions
Shortcoming: in general, the analytical computation of ρ (M(λ)) is not trivial
However, the above shortcoming is a non issue when the objective is to prove the
unconditional stability of an (implicit) scheme
re-write the linear version of equation (2) before time-discretization in matrix form as
du
+ B(∆x)u = 0 (11)
dt
suppose that B is diagonalizable and transform equation (11) into the set of
independent scalar equations
dvm
+ µm (∆x)vm = 0, µm > 0, m = 1, · · · , N
dt
focus on one of the above equations and discretize it in time
apply the scalar form of the Matrix Method for stability analysis: if the conclusion
turns out to be independent of µm , the aforementioned shortcoming is a non issue
example (in class): apply the Matrix Method to determine the stability of a BT
scheme for the linear advection equation
Similarly, the above shortcoming is a non issue when only the maximum eigenvalue µmax m is
needed to conclude the stability analysis – example (in class): apply the Matrix Method to
determine the stability of an FT scheme for the same equation
33 / 75
Formal, Global, and Local Order of Accuracy
Formal order of accuracy measures the orders of accuracy of the

individual space and time approximations separately
Taylor series expansions
modified linear equations
However due to instability, formal order of accuracy may not be
indicative of the actual performance of a method: For example,
recall that a necessary stability condition is
λρ(A) ≤ K ⇒ ∆tρ(A) ≤ K ∆x
and observe that such a stability condition prevents, for example,

fixing ∆t and studying the order of accuracy of the individual space
approximation when ∆x → 0
34 / 75
Besides formal order of accuracy, one way to measure the order of

accuracy is to reduce ∆x and ∆t simultaneously while maintaining
∆t
λ= constant and fixing the initial and boundary conditions
∆x
In this case, a method is said to have global p-th order of accuracy
(in space and time) if
kek∞ ≤ Cx ∆x p = Ct ∆t p , ei = u(xi , t n ) − uin
Cx
for some constant Cx and the related constant Ct =
λp
Other error measures can be obtained by using the 1-norm, 2-norm,
or any vector norm, or if the error is measured pointwise
35 / 75
Determining analytically the global order of accuracy defined above

can be challenging: For this reason, it is usually predicted by
comparing two different numerical solutions obtained using the same
numerical method but two different values of ∆x
log(ke2 k∞ /ke1 k∞ ) log(ke2 k∞ /ke1 k∞ )

p= =
log(∆x2 /∆x1 ) log(∆t2 /∆t1 )
where eli = u(xi , t n ) − uin is the error for ∆xl and ∆tl = λ∆xl ,
l = 1, 2
36 / 75
Another way to measure the order of accuracy is to assume that the

solution is perfect at time t n — that is, uin = u(xi , t n ) ∀i, which is
usually true for n = 0 — and measure the local (in time) truncation
error induced by a single time-step
u(xi , t n+1 ) − uin+1

ēi =
∆t
∆t
Now, let ∆t → 0 and ∆x → 0 while maintaining λ = , and the
∆x
initial and boundary conditions fixed: Then, a method is said to
have local p-th order of accuracy (in space and time) if
kēk∞ ≤ Cx ∆x p = Ct ∆t p
Cx
for some constant Cx and the related constant Ct = p
λ
The local order of accuracy is easy to determine analytically
Example (in class): determine analytically the local order of
accuracy of FTFS for linear advection
37 / 75
Upwind Schemes in One Dimension
Consider a nonlinear scalar conservation law

In 1D, there are right-running waves and left-running waves: For
right-running waves, right is the downwind direction and left is the
upwind direction
Similarly for left-running waves, left is the downwind direction and
right is the upwind direction
Then, every numerical approximation to a scalar conservation law
can be described as
Centered: if its stencil contains equal numbers of points in both
directions
Upwind: if its stencil contains more points in the upwind direction
Downwind: if its stencil contains more points in the downwind
direction
Upwind and downwind stencils are adjustable or adaptive stencils:
Upwind and downwind methods test for wind direction and then,
based on the outcome of the tests, select either a right- or a
left-biased stencil
38 / 75
39 / 75
40 / 75
Upwinding ensures shock avoidance if the shock reverses the wind,

whereas central differencing does not
41 / 75
Upwinding does not ensure shock avoidance if the shock does not
reverse the wind
downwinding on the right above avoids the shock but violates the
CFL condition and thus would create larger errors than crossing the
shock would
42 / 75
General remarks
upwind methods are popular because of their excellent shock
capturing ability
among simple FT or BT methods, upwind methods outdo centered
methods: However, higher-order upwind methods often have no
special advantages over higher-order centered methods
Sample techniques for designing methods with upwind and adaptive
stencils
flux averaging methods
flux splitting methods?
wave speed splitting methods?
43 / 75
Introduction to Flux Splitting
f (u) = f + (u) + f − (u)

Flux splitting is defined as
df + df −
≥ 0, ≤0
du du
Hence, f + (u) is associated with a right-running wave and f − (u) is
associated with a left-running wave
Using flux splitting, the governing conservation law becomes
∂u ∂f + ∂f −
+ + =0
∂t ∂x ∂x
which is called the flux split form
∂f +
Then, can be discretized conservatively using at least one point
∂x
∂f −
to the left, and can be discretized conservatively using at least
∂x
one point to the right, thus obtaining conservation and satisfaction
of the CFL condition
44 / 75
Unfortunately, because in general

+ −
∂f + ∂f −

∂f ∂f
6= and 6=
∂u ∂u ∂u ∂u
flux splitting cannot describe the true connection between fluxes and
waves, unless all waves run in the same direction
if all waves are right-running, the unique physical flux splitting is
f + = f and f − = 0
if all waves are left-running, the unique physical flux splitting is
f − = f and f + = 0
All waves run in the same direction only for (nonlinear) scalar
conservation laws away from sonic points, and for the Euler
equations in the supersonic regime
45 / 75
∂f +
Assume that is discretized with a leftward bias so that the
∂x
approximation at x = xi is centered or biased towards x = xi−1/2
∆fî−1/2
+
+
∂f
≈ for some ∆fî−1/2
+
∂x i ∆x
∂f −
Assume that is discretized with a rightward bias so that the
∂x
approximation at x = xi is centered or biased towards x = xi+1/2
−
∆fî+1/2
−
∂f −
≈ for some ∆fî+1/2
∂x i ∆x
Using forward Euler to perform the time-discretization leads to
n n
u n+1 = uin − λ(∆fˆ+
i + ∆fˆ− )
i−1/2 i+1/2
which is called the flux split form of the numerical approximation

46 / 75
A method in flux split form is conservative if and only if

+n −n
∆fî+1/2 + ∆fî+1/2 n
= gi+1 − gin for some gin (12)
47 / 75

+n −n
∆fî+1/2 + ∆fî+1/2 n
Proof n n
−
uin+1 = uin − λ(∆fî−1/2
+
+ gin − gin + ∆fî+1/2 )
compare with the conservation form uin+1 = uin − λ(fî+1/2

n
− fî−1/2
n
)
n n
−
=⇒ fî+1/2
n
= ∆fî+1/2 + gin , fî−1/2
n
= −∆fî−1/2
+
+ gin
n n
−
=⇒ ∆fî+1/2 = fî+1/2
n
− gin , ∆fî−1/2
+
= −fî−1/2
n
+ gin
47 / 75

+n −n
∆fî+1/2 + ∆fî+1/2 n
Proof n n
−
uin+1 = uin − λ(∆fî−1/2
+
+ gin − gin + ∆fî+1/2 )
compare with the conservation form uin+1 = uin − λ(fî+1/2

n
− fî−1/2
n
)
n n
−
=⇒ fî+1/2
n
= ∆fî+1/2 + gin , fî−1/2
n
= −∆fî−1/2
+
+ gin
n n
−
=⇒ ∆fî+1/2 = fî+1/2
n
− gin , ∆fî−1/2
+
= −fî−1/2
n
+ gin
require now that

−n +n
fî+1/2
n
= fˆ(i+1)−1/2
n
⇒ ∆fî+1/2 + gin = −∆fî+1/2 n
+ gi+1
n n
−
=⇒ ∆fî+1/2
+
+ ∆fî+1/2 n
= gi+1 − gin
Since there are no restrictions on gin , every conservative method has
infinitely many flux split forms useful for nonlinear stability analysis
47 / 75
Example: Design a first-order upwind method for Burgers’ equation

using flux splitting then re-write it in conservation form
for Burgers’ equation, the unique physical flux splitting is
u2 u u
f (u) = = max(0, u) + min(0, u)
2 | {z 2
} | {z 2}
f + (u) f − (u)
48 / 75

using flux splitting then re-write it in conservation form (continue)
a flux split form of Burgers’ equation is
∂u 1 ∂ 1 ∂
+ (max(0, u)u) + (min(0, u)u) = 0
∂t 2 ∂x 2 ∂x
49 / 75

∂u 1 ∂ 1 ∂
+ (max(0, u)u) + (min(0, u)u) = 0
∂t 2 ∂x 2 ∂x
+
∂f
a backward-space approximation of gives
∂x
+n
∆fî−1/2
n
max(0, uin )uin − max(0, ui−1
n n

∂ )ui−1
(max(0, u)u) ≈ =
∂x i ∆x ∆x
49 / 75

∂u 1 ∂ 1 ∂
+ (max(0, u)u) + (min(0, u)u) = 0
∂t 2 ∂x 2 ∂x
+
∂f
∂x
+n
∆fî−1/2
n
n n

∂ )ui−1
(max(0, u)u) ≈ =
∂x i ∆x ∆x
∂f −
a forward-space approximation of gives
∂x
n ˆ −n
∆ f n n
− min(0, uin )uin

∂ i+1/2 min(0, ui+1 )ui+1
(min(0, u)u) ≈ =
∂x i ∆x ∆x
49 / 75

∂u 1 ∂ 1 ∂
+ (max(0, u)u) + (min(0, u)u) = 0
∂t 2 ∂x 2 ∂x
+
∂f
∂x
+n
∆fî−1/2
n
n n

∂ )ui−1
(max(0, u)u) ≈ =
∂x i ∆x ∆x
∂f −
a forward-space approximation of gives
∂x
n ˆ −n
∆ f n n
− min(0, uin )uin

∂ i+1/2 min(0, ui+1 )ui+1
(min(0, u)u) ≈ =
∂x i ∆x ∆x
combining these with an FT approximation yields
λ
uin+1 = uin − (max(0, uin )uin − max(0, ui−1
n n
)ui−1 )
2
λ n n
− (min(0, ui+1 )ui+1 − min(0, uin )uin )
2
49 / 75

the reader can check that the first-order upwind method described in
the previous page can be re-written in conservation form using
1 n 2
gin = (ui )
2
50 / 75
Introduction to Wave Speed Splitting
In contrast with flux splitting, wave speed splitting uses the

governing equations in non conservation form and tends to yield non
conservative approximations
51 / 75
In contrast with flux splitting, wave speed splitting uses the

governing equations in non conservation form and tends to yield non
conservative approximations
Hence in most cases, flux splitting is preferred over wave speed
splitting ... except when the flux function has the property
df
f (u) = u = a(u)u
du
which means that f (u) is a homogeneous function of degree 1
(recall Euler’s theorem which states that a differentiable function
f (u) is a homogeneous function of degree p if and only if
(df /du) u = pf (u)): This property makes flux splitting and wave
speed splitting closely related
51 / 75
For scalar conservation laws, wave speed splitting can be written as
a(u) = a+ (u) + a− (u)

a+ (u) ≥ 0, a− (u) ≤ 0
52 / 75
a(u) = a+ (u) + a− (u)

a+ (u) ≥ 0, a− (u) ≤ 0
Then, the scalar conservation law can be written as
∂u ∂u ∂u
+ a+ + a− =0
∂t ∂x ∂x
which is called the wave speed split form
52 / 75
a(u) = a+ (u) + a− (u)

a+ (u) ≥ 0, a− (u) ≤ 0
Then, the scalar conservation law can be written as
∂u ∂u ∂u
+ a+ + a− =0
∂t ∂x ∂x
which is called the wave speed split form

∂u
Then, a+ can be discretized conservatively using at least one
∂x
∂u
point to the left, and a− can be discretized conservatively using
∂x
at least one point to the right, thus obtaining satisfaction of the
CFL condition
52 / 75
For vector conservation laws such as the Euler equations, one can
split the Jacobian matrix as follows
A(W ) = A+ (W ) + A− (W )
where the eigenvalues of A+ are non negative and those of A− are
non positive
A+ ≥ 0, A− ≤ 0
(recall that A+ and A− are obtained by computing and splitting the
eigenvalues of A)
The wave speed split form of the Euler equations can then be
written as
∂W ∂W ∂W
+ A+ + A− =0
∂t ∂x ∂x
∂W
Again, A+ can then be discretized conservatively using at least
∂x
∂W
one point to the left, and A− using at least one point to the
∂x
right, thus obtaining satisfaction of the CFL condition
53 / 75
Back to scalar conservation laws

If f (u) is a homogeneous function of degree 1, then from Euler’s
theorem it follows that
f (u) = a(u)u ⇒ one may propose f ± (u) = a± (u)u
54 / 75

df +
However, the above splitting may or may not satisfy ≥ 0 and
du
df −
≤ 0, and therefore may or may not correspond to a flux
du
splitting
54 / 75

df +
However, the above splitting may or may not satisfy ≥ 0 and
du
df −
≤ 0, and therefore may or may not correspond to a flux
du
splitting
Why?
54 / 75
∂u
Assume that a+ is discretized with a leftward bias so that the
∂x
approximation at x = xi is centered or biased towards x = xi−1/2
n
uin − ui−1
n

+ ∂u +n
a ≈ ai−1/2
∂x i ∆x
∂u
Assume that a− is discretized with a rightward bias so that the
∂x
approximation at x = xi is centered or biased towards x = xi+1/2
n n
− uin

− ∂u −n ui+1
a ≈ ai+1/2
∂x i ∆x
Using forward Euler to perform the time-discretization leads to
n n
−
uin+1 = uin − λai+1/2 n
(ui+1 +
− uin ) − λai−1/2 (uin − ui−1
n
)
which is called the wave speed split form of the numerical

approximation 55 / 75
The flux split form and wave speed form are connected via
n n
± ±
∆fî+1/2 = ai+1/2 n
(ui+1 − uin )
From the above relation and equation (12), it follows that a method
in wave speed split form is conservative if and only if
n n
+ − n
(ai+1/2 + ai+1/2 )(ui+1 − uin ) = gi+1
n
− gin for some flux function gin
(13)
Hence, the transformation from conservation form to wave speed
form and vice versa is
−n +n
fî+1/2
n
= ai+1/2 n
(ui+1 − uin ) + gin , fî−1/2
n
= −ai−1/2 (uin − ui−1
n
) + gin
(14)
56 / 75
n n
−
(ui+1 +
− uin ) − λai−1/2 (uin − ui−1
n
)
The above notation for the wave speed split form is the standard
notation when wave speed splitting is used to derive new (forward
Euler) approximation methods
57 / 75
n n
−
(ui+1 +
− uin ) − λai−1/2 (uin − ui−1
n
)
Wave speed split form is also often used as a preliminary step in
nonlinear stability analysis, in which case the standard notation is
n n
−
uin+1 = uin + Ci+1/2
+ n
(ui+1 − uin ) − Ci−1/2 (uin − ui−1
n
)
57 / 75
n n
−
(ui+1 +
− uin ) − λai−1/2 (uin − ui−1
n
)
Wave speed split form is also often used as a preliminary step in
nonlinear stability analysis, in which case the standard notation is
n n
−
+ n
(ui+1 − uin ) − Ci−1/2 (uin − ui−1
n
)
Hence
n n n n n n
+ − − + − +
Ci+1/2 = −λai+1/2 , Ci−1/2 = λai−1/2 ⇔ Ci+1/2 = λai+1/2
(15)
57 / 75
From (15), it follows that if a method is derived using wave speed

splitting and not just written in wave speed split form, the splitting
underlying (13) can also be written as
λa(u) = C − (u) − C + (u), C + (u) ≥ 0, C − (u) ≥ 0
Then, the conservation condition (13) becomes

n n
− + n
(Ci+1/2 − Ci+1/2 )(ui+1 − uin ) = λ(gi+1
n
− gin )
And equations (14) become

n n
−
λfî+1/2
n +
= −Ci+1/2 n
(ui+1 −uin )+λgin , λfî−1/2
n
= −Ci−1/2 (uin −ui−1
n
)+λgin
58 / 75
Nonlinear Stability Analysis
Focus is set here on explicit FT difference approximations

Recall that unstable solutions exhibit significant spurious oscillations
and/or overshoots
Recall also that linear stability analysis focuses on these oscillations
and relies on the Fourier series representation of the numerical
solution: It requires only that this solution should not blow up, or
more specifically, that each component in its Fourier series
representation should not increase to infinity
because of linearity, this is equivalent to requiring that each
component in the Fourier series should shrink by the same amount or
stay constant at each time-step
59 / 75
Similarly, nonlinear stability analysis focuses on the spurious

oscillations of the numerical solution, but without representing it by
a Fourier series
1 The total variation of f defined on [a, b] ⊂ R is a measure of the one-dimensional arclength

of the curve with parametric equation x → f (x), for x ∈ [a, b])
60 / 75

a Fourier series
it can require that the overall amount of oscillation remains bounded,
which is known as the Total Variation 1 Bounded (TVB) condition

60 / 75

a Fourier series
it can also require that the overall amount of oscillation, as measured
by the total variation, either shrinks or remains constant at each
time-step this is known as the Total Variation Diminishing (TVD)
condition


60 / 75

a Fourier series
it can also require that the overall amount of oscillation, as measured
by the total variation, either shrinks or remains constant at each
time-step this is known as the Total Variation Diminishing (TVD)
condition

however, whereas not blowing up and shrinking are equivalent

notions for linear equations, these are different notions for nonlinear
equations: In particular, TVD implies TVB but TVB does not
necessarily imply TVD

60 / 75
Monotonicity Preservation
The exact solution of a scalar conservation law on an infinite spatial

domain is monotonicity preserving: If the initial condition is
monotone increasing (decreasing), the solution is monotone
increasing (decreasing) at all times
61 / 75

Suppose that a numerical approximation inherits this monotonicity
preservation property: Then, if the initial condition is monotone, the
numerical solution cannot exhibit a spurious oscillation
61 / 75

Suppose that a numerical approximation inherits this monotonicity
preservation property: Then, if the initial condition is monotone, the
numerical solution cannot exhibit a spurious oscillation
Monotonocity preservation was first suggested by the Russian
scientist Godunov in 1959: It is a nonlinear stability condition, but
not a great one for the following reasons:
it does not address the case of nonmonotone solutions
it is a too strong condition:
it does not allow even an insignificant spurious oscillation that does
not threaten numerical stability
attempting to purge all oscillatory errors, even the small ones, may
cause much larger nonoscillatory errors
61 / 75
Godunov’s theorem: For linear schemes, that is, schemes that can
be described by
XK2
ujn+1 = n
γm uj+m
m=−K1
(NOT to be confused with linear problems), monotonicity

preservation leads to first-order accuracy at best
62 / 75
Total Variation Diminishing
TVD was first proposed by the american applied mathematician

Amiram Harten in 1983 as a nonlinear stability condition
The total variation of the exact solution may be defined as follows
∞
X
TV (u(·, t)) = sup |u(xi+1 , t) − u(xi , t)|
all possible sets of samples xi i=−∞
63 / 75
TVD was first proposed by the american applied mathematician

Amiram Harten in 1983 as a nonlinear stability condition
The total variation of the exact solution may be defined as follows
∞
X
TV (u(·, t)) = sup |u(xi+1 , t) − u(xi , t)|
all possible sets of samples xi i=−∞
Laney and Caughey (1991):

The total variation of a function on an infinite domain is a sum of
extrema — maxima counted positively and minima counted negatively
— with the two infinite boundaries always treated as extrema and
counting each once, and every other extrema counting twice
63 / 75
64 / 75
Numerical effects that can cause the total variation to increase
65 / 75
The exact solution of a scalar conservation law is TVD

TV (u(·, t2 )) ≤ TV (u(·, t1 )) , ∀t2 ≥ t1
66 / 75

TV (u(·, t2 )) ≤ TV (u(·, t1 )) , ∀t2 ≥ t1
What about the numerical solution of a scalar conservation law?
66 / 75

TV (u(·, t2 )) ≤ TV (u(·, t1 )) , ∀t2 ≥ t1

The total variation of a numerical approximation at time t n may be
equally defined as
∞
X
TV (u n ) = n
|ui+1 − uin |
i=−∞
it is the sum of extrema — maxima counted positively and minima

counted negatively — with the two infinite boundaries always treated
as extrema and counting each once, and every other extrema
counting twice
66 / 75

TV (u(·, t2 )) ≤ TV (u(·, t1 )) , ∀t2 ≥ t1

The total variation of a numerical approximation at time t n may be
equally defined as
∞
X
TV (u n ) = n
|ui+1 − uin |
i=−∞
it is the sum of extrema — maxima counted positively and minima

counted negatively — with the two infinite boundaries always treated
as extrema and counting each once, and every other extrema
counting twice
Now, a numerical approximation inherits the TVD property if
∀n, TV u n+1 ≤ TV (u n )

66 / 75
Important result: TVD implies monotonicity preservation and

therefore implies nonlinear stability
Proof: Suppose that the initial condition is monotone
the TV of the initial condition is u∞ − u−∞ if it is monotone
increasing and u−∞ − u∞ if it is monotone decreasing
if the numerical solution remains monotone, TV = cst; otherwise, it
develops new maxima and minima causing the TV to increase
if the approximation method is TVD, this cannot happen and
therefore the numerical solution remains monotone
TVD can be a stronger nonlinear stability condition than the
monotonicity preserving condition
67 / 75
Drawback: Clipping phenomenon (illustrated with the linear

advection of a triangle-shaped initial condition)
68 / 75
Drawback: Clipping phenomenon (illustrated with the linear

advection of a triangle-shaped initial condition)
The TV should increase by ∆x between time-steps but a TVD

scheme will not allow this ⇒ clipping error here this error is O(∆x)
because it happens at a nonsmooth maximum, but for most smooth
extrema it is O(∆x 2 )
68 / 75
Summary of what should be known about TVD

in practice, most attempts at constructing a TVD scheme end up
enforcing stronger nonlinearity stability conditions such as the
positivity condition discussed next
TVD implies monotonicity preservation: This is desirable when
monotonicity preservation is too weak but less desirable when
monotonicity preservation is too strong given that TVD can be
stronger
TVD tends to cause clipping errors at extrema: In theory, clipping
does not need to occur at every extrema — since, for example, the
local maximum could increase provided that a local maximum
decreased or a local minimum increased or a local
maximum-minimum pair disappeared somewhere else — and may be
only moderate when it occurs: However, in practice, most TVD
schemes clip all extrema to between first- and second-order accuracy
in theory, TVD may allow large spurious oscillations but in practice it
rarely does — in any case, it does not allow the unbounded growth
type of instability
69 / 75
Positivity
Recall that the wave speed split form of a FT scheme is given by

n n
−
uin+1 = +
uin + Ci+1/2 n
(ui+1 − uin ) − Ci−1/2 (uin − ui−1
n
)
n n
+ −
Ci+1/2 ≥ 0 and Ci+1/2 ≥0
Suppose that a given FT numerical scheme can be written in wave

speed split form with
n n n n
+ − + −
Ci+1/2 ≥ 0, Ci+1/2 ≥0 and Ci+1/2 + Ci+1/2 ≤ 1 ∀i (16)
Condition (16) above is called the positivity condition (also proposed

first by Harten in 1983)
What is the connection between the positivity condition and the
nonlinear stability of a scheme?
70 / 75
Positivity
The answer is: The positivity condition implies TVD
71 / 75
Positivity
The answer is: The positivity condition implies TVD

Example: FTFS applied to the nonlinear advection equation
∂u ∂u n
+ a(u) = 0 is positive if −1 ≤ λai+1/2 ≤0
∂t ∂x
Proof:
FTFS can be written in wave speed split form for the purpose of
nonlinear stability analysis as follows
n n
−
+ n
(ui+1 − uin ) − Ci−1/2 (uin − ui−1
n
)
n n
+ n −
where Ci+1/2 = −λai+1/2 and Ci−1/2 =0
n
n + n
λai+1/2 ≤ 0 ⇒ Ci+1/2 = −λai+1/2 ≥0
n n
+ − n
Ci+1/2 + Ci+1/2 = −λai+1/2 and therefore the condition (16)
n
becomes in this case −1 ≤ λai+1/2 ≤0
also, note that the positivity condition is in this case equivalent to
the CFL condition
71 / 75
Multidimensional Extensions
The extension to multiple dimensions of the computational part of

the material covered in this chapter may be tedious in some cases
but is straightforward (except perhaps for the characteristic theory)
The expressions of the Euler equations in 2D and 3D can be
obtained from Chapter 2 (as particular cases of the expression of the
Navier-Stokes equations in 3D)
For simplicity, the focus is set here on the 2D Euler equations
∂W ∂Fx ∂Fy
+ (W ) + (W ) = 0
∂t ∂x ∂y
72 / 75
2D structured grid
73 / 75
∂W ∂Fx ∂Fy
+ (W ) + (W ) = 0
∂t ∂x ∂y
For the above 2D Euler equations, the equivalent of equation (2) on

a 2D structured grid is
\n
∂W
∆t = −λx (Fbxni+1/2,j − Fbxni−1/2,j ) − λy (Fbyni,j+1/2 − Fbyni,j−1/2 )
∂t i,j
where
∆t ∆t
λx = , λy =
∆xi ∆yj
∆xi = xi+1/2,j − xi−1/2,j ∀j, ∆yj = yi,j+1/2 − yi,j−1/2 ∀i
and Fbxi+1/2,j and Fbyi,j+1/2 are constructed exactly like fî+1/2 in 1D

74 / 75
For example, a 2D version of FTCS has the following conservative

numerical fluxes
1
Fbxni+1/2,j n n

= Fx (Wi+1,j ) + Fx (Wi,j )
2 
(ρvx )i+1,j + (ρvx )i,j
2
1 (ρv x )i+1,j + (ρvx2 )i,j + pi+1,j + pi,j 
=  
2  (ρvx vy )i+1,j + (ρvx vy )i,j 
(Evx )i+1,j + (Evx )i,j + (pvx )i+1,j + (pvx )i,j
1
Fbyni,j+1/2 n n

= Fy (Wi,j+1 ) + Fy (Wi,j )
2 
(ρvy )i,j+1 + (ρvy )i,j
1 (ρv x vy )i,j+1 + (ρvx vy )i,j

=  
2 (ρvy2 )i,j+1 + (ρvy2 )i,j + pi,j+1 + pi,j 
(Evy )i,j+1 + (Evy )i,j + (pvy )i,j+1 + (pvy )i,j
75 / 75

CA AA214 Ch6

Uploaded by

Copyright:

Available Formats

CA AA214 Ch6

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

CA AA214 Ch6

Uploaded by

Copyright:

Available Formats

AA214: NUMERICAL METHODS FOR COMPRESSIBLE FLOWS 1 / 75

AA214: NUMERICAL METHODS FOR

1 Conservative Finite Difference Methods in One Dimension

2 Forward, Backward, and Central Time Methods

3 Domain of Dependence and CFL Condition

4 Linear Stability Analysis

5 Formal, Global, and Local Order of Accuracy

6 Upwind Schemes in One Dimension

7 Nonlinear Stability Analysis

Note: The material covered in this chapter equally applies to scalar

Also suppose that time is divided into time-intervals [t n , t n+1 ]

One interpretation of the finite difference approach (2) and the

which clearly describes a conservation law (hint: recall the mean

Not every finite difference method can be written in conservation

Like many approximation methods, conservative finite difference

so that from (2) one has

=⇒ the solution of a system of equations is required at each

Like many approximation methods, conservative finite difference

so that from (2) one has

=⇒ the solution of a system of equations is required at each

Like many approximation methods, conservative finite difference

so that from (2) one has

=⇒ only function evaluations are incurred at each time-step

Summary: typical stencil diagram

Like any proper numerical approximation, proper finite difference

fˆ(u, ..., u) = f (u)

=⇒ in this case, the conservative numerical flux fˆ is said to be consistent

Forward Time (FT) conservative finite difference methods

Backward Time (BT) conservative finite difference methods

Central Time (CT) conservative finite difference methods correspond

Recall the theory of characteristics: A point in the x − t plane is

Hence, the physical domain of dependence and physical range of

The direct numerical domain of dependence of a finite difference

The Courant-Friedrichs-Lewy or (CFL) condition

Any numerical method that violates the CFL condition misses

Consider first the linear advection problem

This is because for a > 0, FTFS violates the CFL condition

uin+1 = (1 + λa)uin − λaui+1

However, it almost always blow up (as will be seen in a homework):

However, it almost always blow up (as will be seen in a homework):

hence, the CFL condition is

if K1 = K2 = K , the previous CFL condition becomes

For scalar conservation laws, the CFL condition translates into a

For scalar conservation laws, the CFL condition translates into a

λρ(A) is called the CFL number or the Courant number

Hence, a smart or adaptive stencil can be useful for the case

Unstable solutions exhibit significant spurious oscillations and/or

Unstable solutions exhibit significant spurious oscillations and/or

For a discrete solution u(xi , t n ), the Fourier series is obtained by

Recall that samples can only support wavelengths of 2∆x or longer

Recall that samples can only support wavelengths of 2∆x or longer

From (9), (10), and Euler’s formula e I θ = cos θ + I sin θ it follows

Hence, each term of the Fourier series can be written as

does not depend on n: However, it depends on λ (since uin+1 and uin

Finally, assume that Cm0 = 1 (for example): This leads to

uinm = Gmn (λ)e I φm i

Shortcomings of the von Neumann stability analysis method

Shortcomings of the von Neumann stability analysis method

u n+1 = M(λ)u n , where u n = (uon u1n · · · uNn )T