MSDM5004 Final Project
Due time: 31th May 2022, 9:00 AM
May 15, 2022
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,
(; ) = −ln

0
= ∫(, ).
(3)
The projection can be thought of as a function of obtained at the particular angle .
Our goal is to reconstruct (, ) from the various projections (; ) obtained at
different angles . Consider the coordinate systems in Figure 1. Associated with the
object being scanned are the "fixed" xy-coordinates. We have a second coordinate system
aligned with the x-ray beams, with coordinate along the beam path and perpendicular
to them. The angle is simply the angle between + and + (and between + and +).
These coordinates are related by the familiar relations
= cos + sin (4)
= − sin + cos. (5)




FIGURE 1 CT scan geometry. The path of a typical beam is shown by the dotted
line; its attenuation is represented by a single point in the projection (; ).

Note that Jacobian matrix
=
(, )
(, )
=
[











]




= [
cos sin
− sin cos
]
(6)
is just the rotation matrix and the Jacobian determinant
det() = |
(, )
(, )
| = |
cos sin
−sin cos
| = 1
(7)
Let's look at the Fourier transform of (, ), which is simply
ℱ[ (, )] = (, ) =
1
2
∫ ∫ (, )(+)
∞
−∞

∞
−∞
.
(8)
Define (, ) = ((, ), (, )) and
= cos + sin (9)
= − sin + cos , (10)
which represents a rotation by the same angle in the transformed space, as shown in
Figure 2.










FIGURE 2 The relevant coordinate systems in transform space.
It can be shown that
+ = + (11)
Then
( cos − sin , sin + cos )
=
1
2
∫ ∫ (, )(+)det(−1)
∞
−∞

∞
−∞

=
1
2
∫ ∫ (, )(+)
∞
−∞

∞
−∞

(12)
For fixed , we do not care how the attenuation varies with because the measured
intensity is due to the total attenuation along the path, which is only determined by the
zero “frequency” component. Setting = 0. we have
( cos , sin) =
1
2
∫ ∫ (, )
∞
−∞

∞
−∞

=
1
2
∫ (∫ (, )
∞
−∞
)
∞
−∞
=
1
2
∫ (; )
∞
−∞
.
(13)
This is known as the projection theorem, and simply states that the Fourier transform of
the projection at is the Fourier transform of the attenuation function along in
transform space.
Measuring the projections (; ) at several different angles will yield the transform of
the attenuation function throughout the two-dimensional plane, one radial line at a time.
In principle, the attenuation function could be obtained by evaluating the inverse Fourier
transform. There is, however, a substantial obstacle to this straightforward approach ---
the points at which the transform is known are not the same as those needed for the
inverse transform.
In building the Fourier transform from the projections, knowledge of the transform was
acquired along the radial coordinate at different angles, as illustrated in Figure 3. But the
grid points of a polar coordinate system do not coincide with the grid points of a
Cartesian system, although these are precisely where they're needed if the Fourier
transform is to be evaluated by the FFT. Increasing the number of grid points in a
projection, or the number of projections, does not help solve the problem since the points
still will not be where they are needed. Ultimately, interpolation must be used. A high
order interpolation will require extensive calculations, in addition to the two-dimensional
FFT needed to obtain the attenuation function in coordinate space. Any interpolation
introduces additional error, and interpolation in transform space is particularly tricky ---
the error at one grid point, when transformed back to coordinate space, is spread over the
entire plane. The image thus obtained is the sum of the true image plus error
contributions from every interpolated point in transform space. As a result of these
difficulties, direct inversion is rarely used. Rather, so-called back projection methods are
employed. In effect, these methods perform the (necessary!) interpolation in coordinate
space rather than in transform space.











To see how this works, let's begin by writing the attenuation function as
(, ) =
1
2
∫ ∫ (, )
−(+)
∞
−∞

∞
−∞
.
(14)
Writing x and y in terms of the polar coordinates r and (see Figure 4), and similarly,
and in terms of polar coordinates and , we find
+ = cos cos + sin sin = cos( − ), (15)
so that (, ) becomes
( cos , sin ) =
1
2
∫ ∫ ( cos , sin)− cos(−)
2
0

∞
0

FIGURE 3 The Fourier transform is known on the polar grid, indicated by the black
dots, while the straightforward FFT requires the data to be known on a Cartesian grid.
=
1
2
∫ ∫ ( cos , sin)− cos(−)

0

∞
0
+
1
2
∫ ∫ ( cos , sin )− cos(−)
2


∞
−∞

=
1
2
∫ ∫ ( cos , sin)− cos(−)

0

∞
0
+
1
2
∫ ∫ (− cos ,− sin)+ cos(−)

0

∞
0

=
1
2
∫ (∫ ( cos , sin)− cos(−)
∞
0

0
+ ∫ (− cos ,− sin)+ cos(−)
∞
0
).
(16)











The above integral means we are integrating over infinite straight lines with angle
between the lines and the axis varying from 0 to . From Figure 2 and Figure 4, we
can see that if we let = , then the line is just the axis.
In the last integral, we might be tempted to make the substitution → −. But is the
radial polar coordinate and is necessarily positive --- to avoid confusion, we should leave
it positive. However, if the integral is done last, the angle is fixed while performing the
integral. And at fixed , lies along . Thus, in the next-to-last integral we make the
substitution → , while in the last we use → − and reverse the limits of
FIGURE 4 Geometry relating (, ) and (, ) to .
integration . By our definition of these coordinate systems, we also have that
cos( − ) = . Thus we find
( cos , sin )
=
1
2
∫ (∫ ( cos , sin )
−
∞
0

0
+ ∫ ( cos , sin)
−(−)
0
−∞
)
=
1
2
∫ ∫ ( cos , sin)
−||
∞
−∞

0

= ∫ (; )

0
,
(17)
where
(; ) =
1
2
∫ ( cos , sin)
−||
∞
−∞

(18)
is the modified projection at .
Let's review what we've done. Taking the Fourier transform of the projections, we know
the Fourier transform of the attenuation function on a polar grid in transform space. To
recover the function in coordinate space, we need to perform the inverse transform, but
this is easily done only in Cartesian coordinates. To avoid interpolation in transform
space, the inverse transform is written in terms of polar coordinates. Then, instead of
integrating from 0 to 2, we integrate from 0 to and extend the "radial integration" to
negative values. In this way, we're still integrating over the entire two-dimensional plane,
but the "radial" integral now has the limits of a Cartesian coordinate. The price we pay is
the presence of the factor || in the integrand --- but the integral has the form of an
inverse Fourier transform and can be evaluated by the FFT!
What, then, are the actual steps to be taken to reconstruct the image? One possible
scheme is the following:
0. Initialize everything. Zero the entries of the image array.
1. Compute/import the projection (; ). The projection will be known at , .
2. Use FFT to obtain ( cos , sin) by Equation (13).
3. Evaluate (; ) by Equation (18) with inverse FFT. This modified
projection will be known at the same , as the original projection.
4. Evaluate the contribution to the image from this one projection.
4a. For entry (, ) of the image array, determine by = cos( − ).
4b. Approximate the modified projection at by interpolation of the (; ).
(Linear interpolation is usually sufficient.)
4b. Add this contribution to the integral to the image array, Equation (17).
4c. Repeat steps 4a, 4b, and 4c for all (, ) in the image array.
5. Display the image obtained thus far. (Optional.)
6. Repeat steps 1 – 5 for all angles in the integration.
7. Display the image obtained and compare it with the actual test object.
With this outline as a guide, construction of a computer program is made easier. But first,
we need to obtain the projections (Step 1). In practice, of course, this is the data that
would be obtained from the CT scanning device. In this project, instead, we will compute
the projections ourselves by using a square test object with side length of 1, as shown in
Figure 5.
















FIGURE 5 A test image to explore computerized tomography.
y



x
When the center of the object is at the origin of the xy plane, its attenuation function is
given by
(, ) = {
−64(−0.25)2(−0.25)2 If − 0.5 ≤ , ≤ 0.5
0 Otherwise

Note that the object is inside a unit circle centered at the origin.
Compute the projections (; ) on a discrete set of (; ) by any numerical integration
method you have learned in this course.
It's not necessary to display the image after every additional projection (Step 5).
However, it's interesting to see the image build in this way.
If the quality of the image you obtained is not good, try to vary the following three
parameters:
 The number of points in each projection
 The range of in each projection
 The number of projections

Note that the first two parameters determine Δ, which is inversely related to the range of
. The third parameter is somewhat independent --- you can always add more
projections until the image quality stabilizes.