Exercise sheet 4 – NLAA– Numerical Linear Algebra with Applications
Instructions: Please submit your scripts by 12noon Friday 27 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 fsolve.m, partd.m, psd.m, eulertest.m.
Let Ω = (0, 1)2 and consider the classical formulation of the following reaction-diffusion −ε∆u(x)+cu(x) =f(x) x∈Ω;
(BVP): u(x) =0 x∈∂Ω; (ε,c∈R+). The variational formulation reads
with A a discrete Laplacian matrix and M the mass matrix. The time-dependent version of (BVP) is included below
(IBVP):
A finite element discretisation yields the system of ODE
u(x, 0) = u0(x),
(VF) :
Findu∈V suchthatforallv∈V,
k(u, v) := εa(u, v) + cm(u, v) = l(v) ,
where V is a suitable function space and
a(u,v)= ∇u·∇v dx, m(u,v)= uv dx, l(v)= fv dx. ΩΩΩ
The Ritz-Galerkin finite element method applied to (VF) then yields a linear system of the form Ku = (εA + cM)u = f,
ut(x,t)−ε∆u(x,t)+cu(x,t) =f(x) u(x,t) =0
x∈Ω,t>0; x∈∂Ω,t>0; x ∈ Ω.
Mdu+Ku=f, u(0)=u0. dt
In the following we will assume that Th denotes a uniform subdivision of Ω into (n − 1) × (n − 1) squares of side h; we also let N = (n − 2)2.
1. (a)
A certain choice of finite element space yields a system matrix K in block-tridiagonal form: K = tridiag[E,T,E],
where E, T ∈ R(n−2)×(n−2) are the symmetric Toeplitz tridiagonal matrices:
ch2 ε 4ch2 ε ch2 ε 4ch2 ε 16ch2 8ε 4ch2 ε
E = tridiag 36 − 3, 36 − 3, 36 − 3 , T = tridiag 36 − 3, 36 + 3 , 36 − 3 .
(b)
Write down A,M in block-tridiagonal form, where the blocks are tridiagonal matrices which you should specify. Explain why E,T are simultaneously diagonalizable with eigenvector matrix V which you should specify. Hence, write down a factorisation of K involving V .
Use the decomposition from part (a) to derive a fast solver for the linear system Ku = f. Write a function file fsolve.m with input ε,c and f ∈ RN and output the solution u of the linear system Ku = f.
(d)
2. (a)
bounds for the preconditioned matrix K−1K. Comment on the optimality of K as preconditioner for K.
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 partd.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 fsolve.m derived in part (b). In each case write down the number of iterations. Comment on your results in relation to part (c).
(c)
The constant c in problem (BVP) 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) := and the resulting linear system will employ the coefficient matrix
(b)
(c)
(d) (e)
k→∞
Write down explicit Euler’s method for the system of ODE in the form
K = εA + M.
α1k(u, u) ≤ k(u, u) ≤ α2k(u, u),
Show that
for some positive constants α1,α2 which you should specify in terms of c,γ,Γ. Hence, derive eigenvalue
Let v0, g ∈ Rn, G ∈ Rn×n be given and assume ρ(G) < 1. Consider the recurrence vk+1 = Gvk + g.
Show by induction that
and therefore that Hence, show that
(1)
vk=Gkv0+ Gj g, j=0
vk =Gkv0+(I−G)−1(I−Gk)g. lim vk = (I − G)−1g.
uk+1 =uk +htΦ(tk,uk), (2) where ht is some constant time-step and where Φ is the increment function which you should specify in
terms of K,M and f.
Write recurrence (2) in the form (1), for a matrix G and vector g which you should specify. Hence find
the limit
u∞ := lim uk. k→∞
Comment on this result, in view of the problems (IBVP), (BVP) described above.
Show that ρ(G) < 1 provided 0 < ht < 2/λ1(M−1K).
Write a script file eulertest.m to verify the result in part (c). In particular, your file should
- definef=M1anduseε=c=1;
- estimate λ1(M−1K) using the power method and find a suitable value of ht;
- use Euler’s method (2) to compute an approximation u ̃ to u(x, T ), where T = 10 on subdivisions Th
with parameters n = 8, 16, 32;
- evaluate the norm ∥u ̃ − u∞∥2 for each value of n.
k−1
Ω
c(x)u(x)v(x)dx,