School of Mathematics
MTH4321, MTH5321, MTH5530
2021 [1]
Assignment 1
Due date: Thursday 1 April, 2020, 6pm (submit the electronic copy of your assignment and the
code via Moodle). Late submissions will be subject to late penalties (see the unit guide for full
details).
This assignment contains 7 questions for a total of 100 marks (equal to 10% of the final mark for this unit).
The relevant output of the programming exercises are to be submitted as part of your electronic submission.
You can either include the outputs as a part of your pdf file (if you can edit the pdf document), or include
them in a separate word document (you need to label the figures and outputs consistently with the question
number). In addition, collect all your Matlab m-files in one zip file and submit this file through Moodle.
(Note: use of Matlab is suggested, but you are free to use another programming language.)
Part I: Theoretical Work
1. Eigenvectors and eigenvalues for finite difference operators (10 marks).
Consider the finite difference operator for the 1D elliptic problem
Ah =
1
h2

−2 1
1 −2 1
1 −2 1
. . .
. . .
. . .
1 −2 1
1 −2

∈ Rm×m,
where
h =
1
m+ 1
.
This matrix can be considered a discrete version of the continuous operator d2/dx2 that acts upon a
function u(x).
(a) Show that the m eigenvectors of Ah are given by the vectors ~x(p) (p = 1, . . . ,m) with components
x
(p)
j = sin(ppijh),
and with eigenvalues
λp =
2
h2
(cos(ppih)− 1).
(Hint: verify for a general component (Ah ~x(p))j , using trigonometric identities.)
(b) Verify that the functions u(p)(x) = sin(ppix) (p = 1, 2, . . .) are eigenfunctions of the continuous
differential operator d2/dx2 on domain [0, 1] with boundary conditions u(0) = 0 and u(1) = 0. (Here,
u(x) is called an eigenfunction of d2/dx2 with eigenvalue λ if d2u(x)/dx2 = λu(x), and u(x) satisfies
the boundary conditions.) What are the eigenvalues? Compare the eigenvectors and eigenvalues for
the discrete and the continuous operators and comment (we would expect the discrete and continuous
eigenvalues and eigenvectors to be similar, if the discrete operator is a useful approximation of the
continuous operator . . . ). Are the discrete and continuous eigenvalues similar for small values of h p?
2. Truncation error of Leapfrog scheme for linear advection equation (10 marks).
Show that the truncation error tnj for the Leapfrog scheme for the linear advection equation
vn+1j − vn−1j
2 ∆t
+ a
vnj+1 − vnj−1
2 ∆x
= 0
is O(∆t2) +O(∆x2). (Recall that the truncation error is obtained by plugging an exact solution into the
finite difference equation.)
3. Numerical stability of Crank-Nicolson for linear advection equation (10 marks).
Determine under which conditions the Crank-Nicolson scheme for the linear advection equation
vn+1j − vnj
∆t
+
a
2
[
vnj+1 − vnj−1
2∆x
+
vn+1j+1 − vn+1j−1
2∆x
]
= 0
is numerically stable.
4. General solution of 2D linear advection equation (10 marks).
We have shown in class that u(x, t) = f(x− a t) is a solution of the linear advection equation,
ut + a ux = 0.
Find an equivalent expression for a general solution u(x, y, t) = g(?, ?) of the 2D linear advection equation
ut + a ux + b uy = 0 (show that the expression you propose is always a solution of the PDE).
5. Numerical stability of FD method for heat equation (10 marks).
Consider finite difference method
vn+1j − vnj
∆t
=
vnj−1 − vn+1j − vnj + vnj+1
∆x2
for heat equation
ut = uxx.
For which values of ∆t is this FD method stable?
2
Part II: Computational Work
6. Finite Difference Method for 2D Elliptic PDEs (20 marks).
(a) Implement a Matlab function with header
function A=build laplace 2D(m)
that constructs a sparse matrix A containing the 2D Laplacian matrix as defined in class and in the
lecture notes. Here, m is the number of interior points per direction, i.e., the number of rows and
columns of A is n = m2.
Test your program using A=build laplace 2D(5) and spy(A), and submit the spy plot as part of
your hardcopy submission.
Programming notes:
– Use a sparse matrix to store A. For example, some of the Matlab commands sparse, speye,
spdiags could be useful. (Use the help command to learn about them.)
– Don’t use the Matlab Kronecker product command kron here, but use more elementary Matlab
commands to put the m×m blocks into A.
(b) Write a matlab script laplace solve.m to solve the PDE boundary value problem
BVP

∂2u
∂x2
+
∂2u
∂y2
= (2− 12x2)(y2 − y4) + (2− 12y2)(x2 − x4)
(x, y) ∈ Ω = [0, 1]× [0, 1]
u(x, y) = 2 on ∂Ω.
The exact solution of this problem can be found in closed form and is given by
u(x, y) = 2 + (x2 − x4)(y2 − y4).
You can use your Matlab function build laplace 2D.m from part (a), or you can download and
use the simpler function build laplace 2D kron.m from Moodle (which generates the 2D Laplacian
matrix using a Kronecker product approach), if you did not get part (a) to work.
Test your program using m = 32 interior points in each direction, and submit mesh plots of the
following grid functions: the numerical approximation V h, the exact solution Uh, and the error
Eh = Uh − V h.
Programming notes:
– Since the boundary conditions are nonzero, the boundary values need to be added to the RHS
of the linear system.
– Make sure to consistently use row lexicographic ordering, e.g., when putting the RHS of the
PDE into the RHS vector of the linear system (don’t forget the h2 factor), or when taking the
solution vector of the linear system and turning it into a matrix corresponding to the grid points
for plotting with the mesh command. (You can also use the reshape command to reshape a
lexicographically ordered vector into a matrix.)
– Once you have set up the matrix and the RHS of the linear system, you can use Matlab’s built-in
backslash operator to solve the linear system: solve the system A~v = ~f using v=A\f.
– Unknowns for the boundary points, with value equal to 2, are not included in your linear systems,
but it is best to add these values to the arrays that you use to plot the results, such that the
boundary values appear in the result plots.
3
– To generate the grids for plotting, for example, the exact solution, you can use this sequence of
commands: x=(0:h:1); y=(0:h:1); [X Y]=meshgrid(x,y);
u=2+(X.^2-X.^4).*(Y.^2-Y.^4); mesh(X,Y,u) (after defining h = 1/(m+ 1)). (Here X and Y
are matrices of size (m + 2) × (m + 2) that contain the x and y coordinates of the grid points,
respectively.)
(c) Perform a grid convergence study, with m = 2, 4, 8, 16, 32, 64, 128. Give a table with ‖Eh‖2,h (the
grid norm) for the different values of m, and ratios between ‖Eh‖2,h for successive values of m. Which
ratio do you expect asymptotically? Do you observe the theoretically predicted O(h2) convergence?
Submit your Matlab script laplace solve.m (and possibly other files with additional functions) to
the Moodle drop box. The code submitted should produce the grid convergence study described
above. Also, include a copy of the code in your assignment package to be submitted in class.
7. Finite Difference method for the 1D advection-diffusion equation (30 marks)
Consider the IVP
∂u
∂t
+ a
∂u
∂x
= η
∂2u
∂x2
on Ω,
Ω : x ∈ R, t ∈ [0.01, 1],
u(x, 0.01) =
1√
0.01
exp
( −x2
4 η 0.01
)
,
with η = 2 and a = 1. Consider the following two numerical methods for this IVP:
(1) The explicit method:
vn+1i − vni
∆t
+ a
vni − vni−1
∆x
= η
vni+1 − 2vni + vni−1
∆x2
.
(2) The semi-implicit method:
vn+1i − vni
∆t
+ a
vni − vni−1
∆x
= η
vn+1i+1 − 2vn+1i + vn+1i−1
∆x2
.
Define two relevant time step sizes as
∆tadv =
∆x
a
and ∆tdiff =
∆x2
2 η
.
It can be shown that Method 1 is stable under the restrictive condition
∆t1 ≤ 1
2
min(∆tadv,∆tdiff ).
Similarly, it can be shown that Method 2 is stable when
∆t2 ≤ ∆tadv.
Therefore, you can use
∆t1 = 0.9
1
2
min(∆tadv,∆tdiff ),
and
∆t2 = 0.9 ∆tadv,
for simulations.
4
(a) Write a matlab m-file Parabolic.m that implements methods (1) and (2) for the given IVP. Truncate
the spatial domain to [−10, 10] (no boundary conditions are needed for this choice of computational
domain since the solution remains 0 at the boundaries for all times of interest.) Make plots of the
solutions with the two methods for n = 81 grid points at t = 1. You can choose your grid points
such that the middle grid point lies at x = 0 and the first and last grid points lie at x = −10 and
x = 10. To compute the solution at t = 1 using an integer number of steps, you can choose
4tnew = 1
ceil
(4t−1) ,
which is slightly below the time step 4t recommended in the question. Submit your matlab file
Parabolic.m to the Moodle drop box on the course webpage. Also, include a copy of the code and
a printout of the plots in your assignment package to be submitted in class.
(b) Give a table of execution times and required number of time steps for number of grid points n =
11, 21, 41, 81, 161, 321 and 641. The table should also list the time step sizes taken by the two methods
for the various values of n, and the advection and dissipation time step sizes defined above. As the
number of grid points is doubled, what happens with the execution time for the two methods, for
large n? Explain in terms of the time step sizes. Which of the two methods is most efficient for large
n, and why?
You are allowed to discuss theoretical and computational assignment problems and solution
strategies with your classmates at a general level (not step-by-step). You are not allowed to
show others your solutions or computer code, or to copy solutions or computer code from others.
You are not allowed to copy solutions or computer code from online resources. All assignment
materials you submit (including written documents, program code and graphical output) must
be strictly your own work. Compliance will be actively monitored.
5

欢迎咨询51作业君