Mini-project – Numerical Linear Algebra – 3NMNLA(4)
Instructions: Please upload your handwritten solutions as a single PDF document by 22 March 5pm. Please upload to Canvas by the same deadline the file eulermini.m as described in Q3.
Plagiarism check: Your submission should be your own work and should not be identical or substantially similar to other submissions. A check for plagiarism will be performed on all submissions.
Assessment: This sheet is assessed, with a maximum contribution to your final mark of 25%.
Background
I. Let Ω ⊂ R
2 be an open, bounded domain with boundary Γ and consider the following reaction-diffusion problem:
where ε > 0. The finite element method (FEM) applied to (BVP) yields a linear system of equations of the form.
where A, M ∈ RN×N are symmetric and positive-definite matrices. In particular,
• A corresponds to a FEM discretisation of the negative Laplacian operator (−∆).
• M represents a FEM discretisation of the identity operator (it is also the matrix arising in FEM discreti-sations of generalised eigenvalue problems).
• The right hand side f ∈ RN is related to the data f(x) and g(x).
Examples of matrices A, M, but also of right hand side vectors f can be generated using poisson2dgen.m for the case Ω = (0, 1)2
, f(x) = 1, g(x) = 0.
A related problem that will be referred to is the generalised eigenvalue problem
II. We associate with (BVP) the following time-dependent (or initial value) problem:
Note that if v(x, t) = v∞(x), then v∞(x) is called the steady-state solution of (IBVP).
A certain FEM discretisation of (IBVP) yields the following system of initial value problems (IVP)
where, unlike the vector u in (1), the vector v has entries which depend on time:
III. Problem (3) is solved numerically by using finite difference approximations for the time derivative of v. In the case of one-step methods, this results in recurrences of the form.
where Vk ≈ v(tk) and where {tk = kh, k = 1, . . . , Nt} is a discrete set of Nt points in time (with h > 0 the time-step). Different choices of Φk correspond to different (one-step) methods of numerical integration.
Outline The aim of this mini-project is to work with the systems of equations (1) and (3) and use NLA to derive a set of theoretical and practical results. In the following, Q1 contains a set of results necessary in Q2, while Q3 suggests a certain practical approach and investigation via Matlab.
1. Let G ∈ R
n×n
satisfy ρ(G) < 1. Define
(a) Show that I − G is invertible.
(b) Prove that
(c) Let v0, g ∈ R
n be given and consider the recurrence
Define wk = g for all k ∈ N and Using (4) and the definition of wk, show that
where H ∈ R2n×2n is a 2-by-2 block matrix that you should identify.
(d) Using the above recurrence for zk, show by induction that for all k ∈ N
Deduce that
2. Let M, A ∈ Rn×n be large, sparse, symmetric and positive-definite matrices. Define K = εA + M, where ε > 0 is a parameter to be specified. Let f ∈ R
n
. The solution v(t) to the system of initial value problems
is approximated at time tk := kh,(k ∈ N) by a sequence Vk generated by approximating the time derivative using a finite difference method:
The parameter h > 0 is called the time-step and the above recurrence is known as the explicit Euler’s method.
(a) Write recurrence (5) as a two-term recurrence in Vk in the form. (4), for a matrix G and vector g which you should specify in terms of M, K,f and h.
(b) Consider the generalised eigenvalue problem (2). By Q18, Examples sheet 4, the eigenpairs of (2) are real.
i. Show that µi are positive for all i.
ii. Show further that ρ(G) < 1 if and only if h < hmax for a value hmax which you should specify in terms of ε and (one of) the eigenvalues µi in (2).
(c) Let h < hmax. Using the results from Q1, find the limit Vk. [10]
3. The implementation of recurrence (5) requires
• the estimation of a suitable time-step h;
• the efficient solution of a certain linear system at every time-step.
We describe below how both tasks should be achieved.
(a) As found in Q2(b), the estimation of a suitable time-step h is related to the eigenvalues of (2). The result of Q18, Examples sheet 4 points to an efficient implementation which avoids the use of complex quantities. Let us write (2) in the symmetric form.
Since the eigenpairs are real, we can write down a suitable single-vector iteration method for the above standard eigenvalue problem (in terms of hatted quantities: etc.), with all the variables arising being real. Re-write this method in terms of the original quantities: A, M, xi
. Implement the resulting algorithm in a function file eigsymgen. Use this function as a local function that provides a suitable estimate for hmax (which was described in Q2(b)).
(b) Let us re-write (5) in the form.
where is a function which you should specify in the comments for your code described below. The above linear system needs to be solved for every k; we do this using the preconditioned Steepest Descent method (PSD) with Vk−1 as initial guess and a tridiagonal preconditioner P comprising the main, sub-and super-diagonals of M. Given that P is large, sparse and s.p.d., identify an efficient solution of linear systems of the form. Pz = r to implement as part of the PSD algorithm. Implement the solution of this system in a local function precond.
(c) Write a Matlab function eulermini to implement explicit Euler’s method in the form. (6), using a fixed time-step h = 0.8hmax, where hmax is approximated by using the suitable eigenvalue solver of single iteration type derived in part (a) above. System (6) should be solved at every time step k as described in part (b). The input to eulermini should be A, M, ε, v0,f and the number of time-steps Nt (in this exact order), while the output should be
• a column vector v containing the (final) vector VNt;
• the time-step h employed;
• a column vector k containing the number of PSD iterations required at every time step in order to bring the norm of the current residual divided by the norm of the initial residual below a tolerance τ = 10−6
.
Your file should contain all necessary files included as local functions (e.g., eigsymgen, psd, precond).
(d) Generate test problems using the file poisson2gen, which takes as input the number n of nodes in both x and y directions and generates the matrices A, M and the right-hand side vector f. This file is posted on the assignment page on Canvas. You should run your file eulermini and record in the table below the total number of PSD iterations (i.e., the sum of the entries in the output vector k) for the range of ε and of input parameters n indicated. You should take v0 = 1 and Nt = 100.
Table 1: SD/PSD performance comparison: total iteration counts.
(e) Comment on the results you obtained in part (d). [25]
Notes on the marking procedure
Questions 1,2 You should ensure that you state clearly any results you use; you should also reference these results, by indicating where they can be found (e.g., typed lecture notes, example sheets, Matlab tasks etc).
Question 3 The following issues with the file eulermini.m may attract penalties:
• file does not use the suggested input and output as described in the question;
• file does not run due to incorrect syntax;
• file runs but output is incorrect;
• file does not contain the necessary local files;
• implementation does not follow a least complexity approach/principle (e.g., generating a dense matrix inverse in order to solve a linear system etc.); performing unnecessary operations etc;
• file does not suppress unnecessary output;
• file does not contain comments/descriptions;
• file is partially or entirely identical to other submitted files.