Blog

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.