Consider the one-dimensional scalar hyperbolic conservation law
Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: THEend8_
MSDM5004 Final Project
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) = 12u
2, Eq. (1) belongs to one of the nonlinear scalar
conservation laws and can be written as
∂u
∂t
+
∂( 12u
2)
∂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.1
π 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+ 12 − fi− 12 ). (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+ 12 = f
+
i+ 12
+ f−
i+ 12
. (8)
As for the positive flux f+,
f+
i+ 12
=
2∑
k=0
wkf
+
k,i+ 12
, (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∑2
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+ 12
=
2∑
k=0
wkf
−
k,i+ 12
, (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∑2
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.
Algorithm 1 The numerical algorithms to compute the exact solution of Eq. (1)
when f(u) = 12u
2
Input: The spatial and time coordinate (x, t); the total number of spatial grid
points n.
Output: The exact solution at (x, t).
1: xt0 = x− 0.5t;
2: if xt0 > 1.0 then
3: xt0 = xt0 − 2floor(0.5(xt0 + 1));
4: end if
5: if xt0 < −1.0 then
6: xt0 = xt0 + 2floor(0.5(−xt0 + 1));
7: end if
8: ay = abs(xt0);
9: for j = 0; j < n; + + j do
10: xj = 2 · j/n;
11: qm = 0.5 + sin(π · xj);
12: xt = xj + qm · t;
13: if ay < xt then
14: us = 0.5 + sin(π · xj);
15: un = us;
16: for k = 0; k! = 10000;+ + k do
17: us = un;
18: x0 = ay − us · t;
19: un = us − (us − sin(π · x0))/(1.0 + π · t · cos(π · x0));
20: if abs(un − us) < 10−15 then
21: u = 0.5 + un · sign(xt0);
22: break;
23: end if
24: end for
25: end if
26: end for
27: return u;
Question 2 Computerized Tomography (CT)
Fourier analysis has widespread applicability in many areas. In this project we shall study
an example --- one that has revolutionized modern medicine.
When Wilhelm Roentgen discovered x-rays in 1895, its importance to the practice of
medicine was quickly recognized. Within months, the first medical x-ray pictures were
taken. For his efforts, Roentgen was awarded the first Nobel Prize in Physics in 1901. In
the 1960s and 1970s, another revolution occurred as x-rays were employed to obtain
detailed internal images of the body, using computerized tomography (CT). Today, CT
scans are an important and commonly used tool in the practice of medicine. In 1979,
Cormack and Hounsfield were awarded the Nobel Prize in Medicine for their efforts. Our
interest in computerized tomography is simple --- CT is fundamentally based upon
Fourier transforms.
As a beam of x-rays passes through a uniform slab of material, the intensity of the beam
decreases in a regular pattern. If we measure that intensity, we find that the intensity to be
= 0
− , (1)
where 0 is the initial intensity, d is the thickness of the slab, and is an appropriate
coefficient. is often called an absorption coefficient, but there are actually many
processes occurring that diminish the intensity of the beam.
In particular, the detector is arranged to measure the intensity of the beam that passes
straight through the object, so that any x-rays that are scattered will not be detected.
(Furthermore, the detector is constructed so that x-rays entering at an angle, and hence
scattered from some other region of the object, are rejected.) The reduction of the signal
in CT is primarily due to scattering out of the forward direction. We'll refer to the loss of
intensity, to whatever cause, as attenuation.
In CT, a two-dimensional slice through an object is imaged. In its simplest configuration,
x-rays are directed in a series of thin pencil-like parallel beams through the image plane,
and the attenuation of each beam measured. Since the object is non-uniform, the
attenuation of the beam varies from point to point, and the total attenuation is given as an
integral along the path of the beam. That is, if (, ) is the two-dimensional attenuation
function, the intensity of each of the beams will be of the form
= 0
−∫(,) , (2)
where is an element along the path. Each parallel beam, offset a distance from the
origin, traverses a different portion of the object and is thus attenuated to a different
extent. From the collection of beams a profile of the object can be built. More to the
point, the projection at can be obtained by evaluating the logarithm of the measured
0⁄ ratio,