Exercise sheet 4 – NLAA– Numerical Linear Algebra with Applications
Instructions: Please submit your scripts by 2pm Friday 28 April using the white box labeled D. Loghin on the first
floor of the Watson Building. Please upload to Canvas by the same deadline a zip archive containing the files fpsq.m,
part-d.m, psd.m, part-e.m, shiftinvert.m.
1. Let Ω = (0, 1)2 and consider the classical formulation of the following reaction-diffusion
(CF) :
{
−ε∆u(x) + cu(x) = f(x) x ∈ Ω;
u(x) = 0 x ∈ ∂Ω; (ε, c ∈ R+).
The weak formulation reads
(WF) :
{
Find u ∈ V such that for all v ∈ V,
k(u, v) := εa(u, v) + cm(u, v) = `(v)
,
where V is a suitable function space and
a(u, v) =
∫
Ω
∇u · ∇v dx, m(u, v) =
∫
Ω
uv dx, `(v) =
∫
Ω
fv dx.
The Ritz-Galerkin finite element method applied to (WF) then yields a linear system of the form
Ku = (εA+ cM)u = f ,
with A a discrete Laplacian matrix and M the mass matrix.
(a) Let Th denote a uniform subdivision of Ω into (n− 1)× (n− 1) squares of side h and let N = (n− 2)2. A certain
choice of finite element space yields the following elemental matrices associated with any square in Th and which
are used in the assembly process of A,M ∈ RN×N , respectively:
Ae =
1
6
4 −1 −2 −1
−1 4 −1 −2
−2 −1 4 −1
−1 −2 −1 4
, Me = h236
4 2 1 2
2 4 2 1
1 2 4 2
2 1 2 4
.
Find the stencil for K := εA+ cM and hence show that the structure of K is block tridiagonal
K = tridiag [E, T,E] ,
where E, T ∈ R(n−2)×(n−2) are the symmetric Toeplitz tridiagonal matrices:
E = tridiag
[
ch2
36
−
ε
3
,
4ch2
36
−
ε
3
,
ch2
36
−
ε
3
]
, T = tridiag
[
4ch2
36
−
ε
3
,
16ch2
36
+
8ε
3
,
4ch2
36
−
ε
3
]
.
(b) Since E, T are simultaneously diagonalizable, a fast solver can be derived for the linear system Ku = f . Write
a function file fpsq.m with input ε, c and f ∈ RN and output the solution u of the linear system Ku = f .
(c) The constant c in problem (CF) is replaced with a variable coefficient c(x) satisfying 0 < γ ≤ c(x) ≤ Γ for all x ∈ Ω. Correspondingly, the weak formulation will employ a modified bilinear form: k̂(u, v) := εa(u, v) + m̂(u, v), m̂(u, v) := ∫ Ω c(x)u(x)v(x)dx, and the resulting linear system becomes K̂u = εA+ M̂. Show that α1k(u, u) ≤ k̂(u, u) ≤ α2k(u, u), for some positive constants α1, α2 depending on c, γ,Γ. Explain the implication of this result with regard to preconditioning the system K̂u = f with K. (d) Let ε = 10−2, c = 1. Download from Canvas the file nrdiff.mat containing three FEM discretizations K̂u = f of the problem in part (c). Write a matlab script part-d.m to solve these three linear systems by using the Preconditioned Steepest Descent method with preconditioner K = εA+ cM . You should modify psd.m so that the preconditioning step is implemented using the fast solver fpsq.m derived in part (b). In each case write down the number of iterations. Comment on your results in relation to part (c). (e) Let f(x) = λu(x), so that (CF) becomes an operator eigenvalue problem. Explain why one the interior eigen- values is λ = 1 + π2/10. The above FEM discretization of (CF) is used to find numerical approximations λ̃ to λ, using the Shift-and-Invert method with a tolerance of 10−8 and shift λ. Modify the script shiftinvert.m to include the fast solver implemented in fpsq.m for the solution of the linear system arising in the Shift-and- Invert method. Write a script part-e.m to compute approximations λ̃ for n = 64, 128, 256 and to display the absolute errors ∣∣∣λ− λ̃∣∣∣ in each case, together with an estimate of the order of convergence for the sequence of approximations. [50]