H

Hyperbolic
b li EEquation
i Advection Equation
Wave equation (Typical hyperbolic equation) ∂φ ∂φ Analytic Solution ：
+u =0
∂ 2φ 2 ∂ φ
2 ∂t ∂x φ = F ( x − ut )
−u =0
∂t 2 ∂x 2 F(x) is an arbitrary function of x .
C fi i f (t,x
Confirming t,x)) = F(t-Δt, x-uΔt), we can understand
d t d th the
By factorizing as ⎛⎜ ∂ − u ∂ ⎞⎟⎛⎜ ∂ + u ∂ ⎞⎟φ = 0 Profile F(x) moves with the speed u.
⎝ ∂t ∂x ⎠⎝ ∂t ∂x ⎠

∂φ ∂φ ● Let’s solve the advection equation

+u =0 as an initial boundary problem.

∂t ∂x Initial
Condition φ =1

The advection equation is the simplest hyperbolic φ=0

φ(t , x) = F (t − Δt , x − uΔt )
equation. x=0 x =1
1 2

Simple
p
Finite Difference Method CFL Number
At the time t = t n and the position x = xi n +1 φin+1 − φin−1
φ = φ − uΔt
n

2Δx
i i
∂φ φn+1 − φin
( )
■Time derivative → i （Forward Difference) 1
∂t Δt
= φin − C φin+1 − φin−1
2
∂φ φn − φn
■Spatial derivative → i +1 i −1 (Center Difference)
∂x 2Δx C = uΔt/Δx is called CFL (Courant Friedrich--Levy) number.
(Courant--Friedrich

Discretized φin +1 − φin φin+1 − φin−1 Physical traveling distance uΔt

form： = −u C=
Δt 2Δx Numerical traveling distance Δx
3 4
 Time and Space
p Cone for
Information travel S
Sample
l PProgram 3
#include "program1.h"
#define N 101
int main() {
t n +1
double f[N], fn[N], x[N], cfl = 0.5;
tn i t j,
int j icnt
i t = 0;
0
i−1 i i+1
n −1
while(icnt < 100) {
t f (j 1 j < N-1;
for(j=1; N 1 j++) {
fn[j] = f[j] − cfl*0.5*(f[j+1] − f[j-1]);
}
If C > 1, the scheme does not have enough information to obtain the for(j=0; j < N; j++) f[j] = fn[j]; /* updating */
value of the next time }
Requirment:
}
Source Code
The condition C ( CFL number ) ≤ 1 must be statisfied. 5 6

For the advection equation For the advection equation

(1/2）
Stability Analysis (1/2） (2/2）
Stability Analysis (2/2）
Neumann’s Method
von Neumann’
φn +1 δφ
δφ
1
(
φn = 1 − C eikΔx − e −ikΔx
2
)
ik ⋅ jΔx
Assuming the perturbation φ = δφ e
n n
j
= 1 − C i sin kΔx
Substituting into φ j = φ j − C (φ j +1 − φ j −1 )
n +1 n 1 n n

2 δφ n +1 δφn = 1 + C 2 sin 2 kΔx

C : CFL number
n +1
Analytically δφ δφ = 1
n

( )
1 δφn +1 δφn ≥ 1
n +1 ik ⋅ jΔx n ik ⋅ jΔx
δφ e = δφ e − C δφn e ik ⋅( j +1) Δx − δφ n e ik ⋅( j −1) Δx The above formula shows
2 It is unstable to apply the forward difference in time and
the center difference in space to the advection equation.
7 8
 Lax Scheme （1/2
1/2）） Sample Program 4
Sample
#include "program1.h"
Advection Equation φ n +1
j
n
j
1
2
(
= φ − C φ nj +1 − φ nj −1 ) #define N 101
int main() {
double f[N], fn[N], x[N], cfl = 0.5;
φ nj →
2
(
1 n
φ j +1 + φ nj −1 ) i t j,
int j icnt
i t = 0;
0
while(icnt < 100) {
f (j 1 j < N-1;
for(j=1; N 1 j++) {
φ nj +1 =
1 n
2
( 1
) (
φ j +1 + φ nj −1 − C φ nj +1 − φ nj −1
2
) fn[j] = 0.5*(f[j+1] – f[j-1])
- cfl*0.5*(f[j+1]
cfl*0 5*(f[j+1] - f[j-1]);
f[j 1]);
φ nj is replaced by the average of φ nj+1 and φnj −1 }
for(j=0;
(j ; j < N;; jj++)) f[j] [j]; /* updating
[j] = fn[j]; p g */
The result is expected to be diffusive. } Source Code
9
} 10

Lax Scheme （2/2

2/2）） 1st order upwind scheme
von Neumann’s stability analysis: Since the solution of the advection equation is the
δφ n +1
δφ = cos kΔx − iC sin kΔx
n profile moving
p g with the velocityy u from the upwind
p to
the down direction, it is natural that we apply the
backward finite difference to the advection term.
δφn +1 δφn = cos 2 kΔx + C 2 sin 2 kΔx

C ≤1
It is understood that the scheme is stable for
n +1
φ nj − φ nj −1
φ = φ −u
n
Δt
Rewriting to 1
(
φ nj +1 = φ nj − C φ nj +1 − φ nj −1 ) j j
Δx
( )
2
= φ nj − C φ nj − φ nj −1
1
(
+ φ nj +1 − 2φ nj + φ nj −1
2
)
11 12
 Sample Program 5
Sample 1st order upwind scheme
#include "program1.h"
#define N 101 (
φ nj +1 = φ nj − C φ nj − φ nj −1 )
− (φ ) ( )
int main() { C C n
= φ nj n
j +1 − φ nj −1 + φ j +1 − 2φ nj −1 + φ nj −1
double f[N], fn[N], x[N], cfl = 0.5; 2 2
i t j,
int j icnt
i t = 0;
0
while(icnt < 100) {
f (j 1 j < N-1;
for(j=1; N 1 j++) { von Neumann’s
N ’ stability
bili analysis:
l i
fn[j] = f[j] – cfl*(f[j] − f[j-1]);
}
for(j=0; j < N; j++) f[j] = fn[j]; /* updating */ δφ nj +1 δφ nj = (1 − C + C cos kΔx )2 + C 2 sin 2 kΔx
} ≤1
}
Source Code
13 14

Cubic Semi
Semi--Lagrangian Method
Other Interpretation
Higher-order Interpolation Fc ( x) = a( x − x j ) 3 + b( x − x j ) 2 + c( x − x j ) + f jn
Linear interpolation for the upwind region x j −1 ≤ x ≤ x j For upwind region,
0≤u f j +1
φ nj − φ nj −1 f j −2
F ( x) = φ +
n n
(x − x j ) Fc ( x j −1 ) = f jn−1
x j − x j −1
j fj
Fc ( x j − 2 ) = f n
j −2
f j −1

Fc ( x j +1 ) = f jn+1
uΔ t Δx
Δx Δx
Δx Δx
Δx
x j −2 x j −1 xj x j +1

f jn+1 − 3 f jn + 3 f jn−1 − f jn− 2 f jn+1 − 2 f jn + f jn−1 2 f jn+1 + 3 f jn − 6 f jn−1 + f jn− 2

a= b= c=
j−1 j 6Δx 3 2Δx 2 6Δx

f jn+1 = Fcn ( x j − uΔt ) = aξ3 + bξ2 + cξ + f jn

At the time t n+Δt , fj n+1 = F n (xj−uΔt) x = x j − uΔt (ξ = −uΔt ) Source Code
15 16
 H
Hermite
it Interpolation
I t l ti CIP Method
t = t n +1
FCIP ( x) = a ( x − xi ) + b( x − xi ) + f x ,i ( x − xi ) + f i
3 2

fi
f jn+1 = FCIP
n
( x j − uΔt ) = aξ3 + bξ 2 + f xn, j ξ + f jn
ui > 0 f x ,i

M t hi Conditions
Matching C diti :
n +1 ∂FCIP
n
F ( xi ) = f i f = ( x j − uΔt ) = 3aξ 2 + 2bξ + f xn, j
∂x
f i −1 x, j
Fx ( xi ) = f x ,i
F ( xi −1 ) = f i −1
f x ,i −1 where ξ = −uΔt
Fx ( xi −1 ) = f x ,i −1
xi −1 Δx xi

a=
1
( f x,i + f x,i −1 ) − 2 ( fi − fi−1 )　, 　b = 1 (2 f　
x ,i + f x ,i −1 ) −
3
( f i − f i −1 )
Source Code
Δx 2 Δx Δx Δx 2 18