1 Question 1
Consider the one-dimensional scalar hyperbolic conservation law
ut + f (u)x = 0, (1)
(I) when set f (u) = 2u, Eq. (1) can be written as
∂u
∂t + 2 ∂u
∂x = 0. (2)
Solve it on the domain X × T = [0, 2] × (0, +∞) with the periodic boundary
condition. And the initial condition is given as
u0(x) = 0.5 + sin(πx). (3)
The space and time domain can be discretized into uniform intervals of size
∆x = 0.01 with xi = i × ∆x, and ∆t with tj = j × ∆t, j = 0, 1, 2, · · · .
(a) Find the exact solution.
(b) Use the forward Euler method and the fourth-order central difference
scheme for the time integration and the spatial discretization, respectively. Ex-
amine the stability property of this numerical method. Write a program to
compute the numerical solution at t = 1 with a suitable timestep and compare
with the exact solution, if the numerical method can keep stable; otherwise,
clarify why the numerical method cannot keep stable.
(c) Use the Crank–Nicolson method and the second-order central difference
scheme for the time integration and the spatial discretization, respectively. Ex-
amine the stability property of this numerical method. Write a program to
compute the numerical solution at t = 1 with a suitable timestep and compare
with the exact solution, if the numerical method can keep stable; otherwise,
clarify why the numerical method cannot keep stable.
(II) When we set f (u) = 1
2 u2, Eq. (1) belongs to one of the nonlinear scalar
conservation laws and can be written as
∂u
∂t + ∂( 1
2 u2)
∂x = 0. (4)
Solve this equation with the same boundary and initial conditions as problem
(I). Note that this equation is known as the inviscid nonlinear Burgers equation.
Different from Eq. (2), discontinuities will appear in the solution of Burger equa-
tion with the time evolution, even if the initial condition is sufficiently smooth.
With the initial condition of this problem, the discontinuity will appear since
t = 1
π . The numerical algorithms to compute the analytical reference solution
can be found in Algorithm 1.
(a) Use the fourth-order Runge-Kutta method (RK4) for the time integration,
and the second-order central difference scheme for the spatial discretization, re-
spectively. Write a program to compute the solutions at t = 0.5
π and 1.1
π with
the timestep ∆t = 0.00001 and resolution N = 200, and compare the numerical
solutions with the analytical references.
(b) Use the fourth-order Runge-Kutta method (RK4) for the time integration,
and the fifth-order WENO5-JS scheme for the spatial discretization, respec-
tively. Write a program to compute the solutions at t = 0.5
π and 1.1https://weibo.com/u/7916053997
π with the
timestep ∆t = 0.00001 and resolution N = 200, and compare the numerical
solutions with the analytical references.
Remarks:
1. In the problem (I.c), if the method can keep stable, use the SOR method to
solve the implicit equation, and matrix calculation is not permitted in your
code. In addition, the parameter ω of the SOR method is recommended
to be set as 1.25.
2. The fifth-order oscillation-free WENO5-JS scheme is widely used in solving
the nonlinear hyperbolic conservation laws. The main steps of WENO5-
JS are as follow:
(a) Eq. (1) can be approximated by a conservative finite-difference scheme
as
dui
dt ≈ − 1
∆x (fi+ 1
2 − fi− 1
2 ). (5)
For a general flux function f (u), one can split it into two parts (i.e., the
positive and negative flux) to ensure the computational stability,
f (u) = f +(u) + f −(u), (6)
where f ±(u) can be calculated by the following equation
f ±(u) = 1
2 (f (u) ± αu), (7)
and the parameter α = max ∂f
∂u .
(b) In the specific discretization, the cell interface numerical flux can be
split as
fi+ 1
2 = f +
i+ 1
2
+ f −
i+ 1
2
. (8)
As for the positive flux f +,
f +
i+ 1
2
=
2X
k=0
wkf +
k,i+ 1
2
, (9)
where,
f +
0,i+1/2 = 1
6
2f +
i−2 − 7f +
i−1 + 11f +
i

,
f +
1,i+1/2 = 1
6
−f +
i−1 + 5f +
i + 2f +
i+1

,
f +
2,i+1/2 = 1
6
2f +
i + 5f +
i+1 − f +
i+2

,
(10)
wk = αk
P2
k=0 αk
, αk = dk
(βk + ε)q , k = 0, 1, 2, (11)
β0 = 1
4
f +
i−2 − 4f +
i−1 + 3f +
i

2 + 13
12
f +
i−2 − 2f +
i−1 + f +
i

2 ,
β1 = 1
4
f +
i−1 − f +
i+1

2 + 13
12
f +
i−1 − 2f +
i + f +
i+1

2 ,
β2 = 1
4
3f +
i − 4f +
i+1 + f +
i+2

2 + 13
12
f +
i − 2f +
i+1 + f +
i+2

2 .
(12)
As for the negative flux f −,
f −
i+ 1
2
=
2X
k=0
wkf −
k,i+ 1
2
, (13)
where,
f −
0,i+1/2 = 1
6
2f −
i+3 − 7f −
i+2 + 11f −
i+1

,
f −
1,i+1/2 = 1
6
−f −
i+2 + 5f −
i+1 + 2f −
i

,
f −
2,i+1/2 = 1
6
2f −
i+1 + 5f −
i − f −
i−1

,
(14)
wk = αk
P2
k=0 αk
, αk = dk
(βk + ε)q , k = 0, 1, 2, (15)
β0 = 1
4
f −
i+3 − 4f −
i+2 + 3f −
i+1

2 + 13
12
f −
i+3 − 2f −
i+2 + f −
i+1

2 ,
β1 = 1
4
f −
i+2 − f −
i

2 + 13
12
f −
i+2 − 2f −
i+1 + f −
i

2 ,
β2 = 1
4
3f −
i+1 − 4f −
i + f −
i−1

2 + 13
12
f −
i+1 − 2f −
i + f −
i−1

2 .
(16)
Here, the parameter q = 2, ε = 10−6, d0 = 0.1, d1 = 0.6, d2 = 0.3. For
more details about the WENO5-JS scheme, you can refer to the reference
paper “Efficient implementation of weighted ENO schemes”, Journal of
computational physics, Volume 126, Issue 1, June 1996, Pages 202-228.