Mini-project – NLAA– Numerical Linear Algebra with Applications
Instructions: Please hand in solutions to Q1, and Q2 by Friday 5pm week 11 (27 March). Make sure to include your name and student ID on the front page. The ESO will set up in week 11 a box labeled 3NLA/3NLA4 outside their office (Poynting S11). You should also upload to Canvas by the same deadline a single file: ipowergenfem.m. Please name your file as suggested in the question.
Plagiarism check: Your handwritten and electronic submissions should be your own work and should not be identical or similar to other submissions. A check for plagiarism will be performed on all submissions.
1. Let A ∈ Rn×n denote a symmetric matrix with distinct eigenvalues. Consider the generalised eigenvalue problem
GEVP: Av = λMv, where M ∈ Rn×n is symmetric and positive definite.
(a) By considering a suitable symmetric factorisation of M, show that the eigenvalues of GEVP are real.
(b) Show that {vi}1≤i≤n can be chosen such that viTMvj =δij,
viTAvj =λiδij,
where (λi,vi) denotes the ith eigensolution (i.e., eigenpair) of GEVP.
(c) Let i denoted a fixed integer with 1 ≤ i ≤ n and define c := αMvi, where α ∈ R\{0}. Consider the modified generalised eigenvalue problem
A c M 0 GEVP-mod: cT 0 u=τ 0T 1 u.
Using part (a), or otherwise, show that GEVP-mod has eigenvectors of the form v
u=a,
where v ∈ {v1, v2, . . . , vn} and a is a real constant which you should specify for each eigenvector
u. Deduce that the eigenvalues τk corresponding to GEVP-mod are non-zero; in particular, show
that two of them satisfy the equation
α2
τ =τ−λi.
Hence, write down all the solutions {(τk,uk)}1≤k≤n+1 of GEVP-mod.
2. Let Ω ⊂ R2 denote a bounded domain with boundary Γ. Consider the following Laplacian eigenvalue
problem:
−∆v(x,y) = λv(x,y) (x,y) ∈ Ω, EVP: n·∇v(x,y) =0 (x,y)∈Γ,
where n denotes the unit outward normal to the boundary of Ω. The boundary condition is known as a Neumann boundary condition and the negative Laplacian operator equipped with this boundary condition is known as a Neumann-Laplacian. It is known that the eigenvalues of EVP are real and non-negative. One can also check that v(x,y) = 1 is an eigenfunction of EVP with corresponding eigenvalue λ = 0. These properties are preserved by the following finite element discretisation of EVP :
FEM-EVP: Av = λMv.
The matrix M is symmetric and positive definite, while A is symmetric and positive-semidefinite. Moreover, there holds A1 = 0, so that (0, 1) is an eigenpair of A. Therefore, an attempt to target the small magnitude eigenvalues using the inverse iteration method is bound to fail since A is singular. To circumvent this issue, we can use an approach based on the result of Q1.
(a) We would like to set up a generalised eigenvalue problem of the form given in Q1(b) such that all its eigenvalues are non-zero, while its smallest eigenvalue is equal to the second smallest eigenvalue of FEM-EVP. Write down a suitable vector c, specifying the choice of vi and α. In particular, you should comment on your choice of α and how you could estimate this parameter in practice in order to guarantee that the inverse iteration method targets the desired eigenpair.
(b) Assuming the setup from part (a), the inverse iteration method requires at every step the solution of linear systems of the form
A c
cT 0z=w.
We would like to solve this problem using a preconditioned iterative solver. Since the above coeffi- cient matrix is indefinite, one can use a dedicated iterative method, such as the Minimum Residual Method (MINRES) or the Symmetric LQ method (SYMMLQ). These methods are implemented in Matlab – type help minres, help symmlq, for details on usage. Both have the requirement that the preconditioner is a positive-definite matrix. Suggest a suitable preconditioner for the above linear system. Prove that your choice of preconditioner is indeed positive-definite.
(c) Modify the file sipowergen.m to provide an implementation of the procedure corresponding to parts (a) and (b). You should name your file ipowergenfem.m and observe the following guidelines:
• your input should be the matrices A,M, together with the tolerance and the maximum steps allowed for the inverse iteration method;
• your output should be the second smallest eigenvalue of FEM-EVP;
• you should include a local file that performs the estimation of the parameter α;
• you should include a local file that generates your choice of preconditioner;
• you should use a preconditioned method (either MINRES or SYMMLQ) at every step of the
inverse iteration method where a linear system needs to be solved.
[A zip file was provided on the assignment page for you to test your code. The file contains several
pairs of matrices A, M arising from finite element discretisations.]
Generic marking scheme. Full marks can be obtained provided you observe the following guidelines:
• you either provide proofs or quote suitable results from lectures, homework sheets or example
sheets;
• when you provide proofs, you justify your steps clearly and rigourously;
• you provide comments, when required (e.g., Q2 (a), (b));
• your work is your own and does not resemble closely other submitted work.
Matlab generic marking scheme. Full marks can be obtained for your Matlab submission if the recommendations in Q2(c) are followed. Penalties will be applied in the following cases:
• file does not run due to incorrect syntax;
• file runs, but output is incorrect due to incorrect code;
• file runs correctly, but does not supress output, i.e., it displays the values of intermediate variables and/or output;
• file is partially or entirely identical to other files.