Numerical methods for 1st order wave equation: ∂tf + a ∂xf = 0
A. Euler explicit (fn+1 - fn)/∆t = -a δxfn
1. Central (FTCS)
(fjn+1 - fjn) = -a ∆t/2∆x (fj+1n - fj-1n) = -C/2 (fj+1n - fj-1n)
fjn+1 = fjn -C/2 (fj+1n - fj-1n)
C = a ∆t/∆x is Courant-Fredrichs-Lewy (CFL) number.
VonNeuman stability: εjn+1 = εjn - C/2 (εj+1n - εj-1n). Error is of form εn = Aneiϴj :
An+1 = An[1 - C/2 (eiϴ - e-iϴ)] = An[1 - iC sinϴ]
Amplification factor:
G = ||An+1 /An || = ||1 - iC sinϴ || = √1+C2 sin2ϴ ≥ 1
Unconditionally unstable Max growth @ π/2
2. One sided difference for -a ∂xfn
Assume that C > 0
( (fjn+1 - fjn) = -C (fj+1n - fjn) or = -C (fjn - fj-1n)
( downwind upwind
( fjn+1 = fjn -C (fj+1n - fjn) ; fjn -C (fjn - fj-1n)
( An+1 = An[1 - C (eiϴ - 1) ] downwind
( An+1 = An[1 - C (1 - e-iϴ)] upwind
case 1 (downwind) G2 = (1+C-C cosϴ)2+ C2sin2ϴ
( ( ( = (1+ C)2-2(C+C2)cosϴ +C2
( ( ( = 1+2C+2C2-2(C+C2)cosϴ
( ( ( = 1+2(C2+C)(1-cosϴ) > 1
case 2 (upwind) G2 = (1 - C + C cosϴ)2+ C2sin2ϴ
( ( ( = (1 - C)2 + 2(C-C2)cosϴ+ C2
( ( ( = 1 - 2C + 2C2 - 2(C2 - C)cosϴ
( ( ( = 1+2(C2-C) (1-cosϴ)
( ( ( ( ( ( ≥0
G<1 if C2 - C< 0 → 0 < C < 1 ∴ Upwind
Upwind is stable if CFL<1; or ∆t < ∆x/a . Downwind is unstable! (C=1 has no error, but
this is not useful)
Animations: Convect_exact.gif, Convect_Central.gif, Convect_Up.gif,exact.giRK_central.gif
This study source was downloaded by 100000899610689 from CourseHero.com on 09-06-2025 01:40:10 GMT -05:00
1
https://www.coursehero.com/file/33324908/Lecture19pdf/
, AerE 546 ( Lecture 19
B. Physical interpretations:
a. Convection is from upwind, not down. If a<0
upwind is other direction (j+1).
b. In one time-step, particle canʼt move more than
the stencil (receiving grid point wonʼt have info):
numerical( ∆t a∆t < ∆x -> C = a∆t/∆x < 1.
physical c. Numerical domain of dependence must include
a∆t physical domain of dependence. Stable if
∆x physical ⊂ numerical. Unstable if numerical ⊂
physical.
∆x
|
a ∆t
1. Another view of upwinding: numerical viscosity: upwind = central + diffusion
δf/δt = -a/∆x(fj - fj-1) = -a/(2∆x)(fj+1 - fj-1) + a/(2∆x)(fj+1 + fj-1 - 2fj)
( ( ( = -aδxfc + ½a∆xδ2xf
½|a∆x| is the numerical diffusivity, say α.
comparison equation ∂tf+a∂xf = α ∂2xf
Dissipation of a sine wave: Let f = sin (k(x-at))A(t). Without diffusion, this is a
convected sine wave. Now:
∂t sin(k(x-at))A(t) +a∂x sin (k(x-at)) A(t) = α ∂2x sin (k(x-at)) A(t) =>
dtA = -αk2 A so sine wave damps with time f = sin(k(x-at))exp(-αk2t).
k=2π/λ so short waves damp rapidly: in one time step decay exponent is
( αk2∆t = 2π2 a∆t∆x/λ2 = 2π2 C (∆x/λ)2
so grid spacing must be small compared to wavelength for accuracy.
2. Automatic upwinding
In CFD a could be positive or negative.
1) Could use an IF statement:
IF(a > 0) δxf = fj - fj-1 / ∆x
ELSE ! δxf = fj+1 - fj / ∆x
2) Without IF statement: a∂xf = ½(a+|a|) ∂xf + ½(a-|a|) ∂xf called `splittingʼ
aδxf = ½(a+|a|) (fj - fj-1) /∆x + ½(a-|a|) ( fj+1 - fj )/∆x
combining terms shows that this is central + diffusion (1st order upwind)
= a/(2∆x)(fj+1 - fj-1) - ½|a∆x|(fj+1 + fj-1 - 2fj) / ∆x2
This study source was downloaded by 100000899610689 from CourseHero.com on 09-06-2025 01:40:10 GMT -05:00
2
https://www.coursehero.com/file/33324908/Lecture19pdf/