Exercise sheet 4 – NLAA– Numerical Linear Algebra with Applications
Instructions: Please submit your scripts by 1pm Monday 29 April using the white box labeled D. Loghin on the first floor of the Watson Building. Please upload to Canvas by the same deadline the file eulertest.m corresponding to question 2(f). 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.
Let Ω ⊂ R2 be an open, bounded domain with boundary Γ and consider the following reaction-diffusion problem −ε∆u(x) + c(x)u(x) = f(x) x ∈ Ω;
(BVP): u(x) =g(x) x∈Γ; where ε > 0 and 0 < c1 ≤ c(x) ≤ c2 for all x ∈ Ω. The variational formulation reads
Findu∈V suchthatforallv∈V, (VF) : k(u, v) := εa(u, v) + mc(u, v) = l(v) ,
where V is a suitable function space and
a(u,v)= ∇u(x)·∇v(x)dx, mc(u,v)= c(x)u(x)v(x)dx, l(v)= f(x)v(x)dx. ΩΩΩ
Let
be a finite dimensional space of functions defined on a subdivision of Ω into triangles of maximum size h. The
Vh = span{φ1(x),φ2(x),··· ,φN(x)} ⊂ V
finite element method applied to (VF) with V replaced by Vh yields a linear system of the form
u(x, 0) = u0(x),
x ∈ Ω.
Ku = (εA + Mc)u = f,
A := [a(φj , φi)], Mc := [mc(φj , φi)], fi = l(φi).
where
The time-dependent problem corresponding to (BVP) is included below
ut(x,t)−ε∆u(x,t)+c(x)u(x,t) =f(x) x∈Ω,t>0; u(x,t) =g(x) x∈Γ,t>0;
(IBVP):
A finite element discretisation of (IBVP) yields the system of ODE
Mdu+Ku=f, u(0)=u0, dt
(1)
where
M := [m(φj , φi)], m(u, v) := u(x)v(x)dx.
Ω
Notation: in the following questions, the notation KN , MN , AN , MNc is sometimes used to indicate that the matrices K=:KN,M=:MN,Mc=:McN,A=:AN arematricesofsizeN×N.
1. Let PN,QN ∈ RN×N be symmetric and positive-definite matrices. We say PN is spectrally equivalent to QN uniformly with respect to N (or uniformly spectrally equivalent) if there exist positive constants γ,Γ independent of N such that for all u ∈ RN \ {0}
γ ≤ uT PN u ≤ Γ. uTQNu
We write PN ≈ QN . The constants γ, Γ are called constants of equivalence.
(a)
(b) (c) (d)
(e) (f )
Consider the set of linear systems
KNuN =fN,
where KN ∈ RN×N with N ∈ N, N ≤ Nmax. Let PN denote a set of preconditioners for KN such that PN ≈ KN . What are the practical implications of uniform spectral equivalence of PN to KN , assuming an iterative solver such as the preconditioned Steepest Descent method is used?
ShowthatifPN ≈QN,thenQN ≈PN. ShowthatifwefurtherhaveQN ≈RN thenPN ≈RN. Show that M ≈ Mc (i.e., MN ≈ McN ), indicating the constants of equivalence.
It is given that for any function uh ∈ Vh there holds
2 η2 2
η1 uh(x)dx ≤ ∇uh(x) · ∇uh(x)dx ≤ h2 uh(x)dx.
ΩΩΩ
Use these inequalities to show that if h > √ε, then MN ≈ εAN .
Use parts (c–d) to show that KN ≈ MN provided h > √ε, indicating the constants of equivalence. Problem (BVP) is solved on the L-shaped domain Ω = (0, 2)2 \ [1, 2)2 using the following data:
f(x)=(1+x+y)(x+y), g(x)=x+y, ε=10−5, c(x)=1+x+y.
Use femsol with quasi-uniform subdivisions (set fem.Hmax=[]) and values of fem.level ranging between 1 and 5 to investigate the performance of the preconditioned Steepest Descent method. As preconditioner you should use the (sparse) diagonal matrix D containing the main diagonal of M. For each linear system write down the value of h2 and the least number of iterations k for which ∥f − Kuk∥2 ≤ 10−6∥f∥2. Comment on your results in view of your statement from part (a). Using parts (b) and (e), discuss what your experiments imply about the spectral equivalence of D and M.
2. (a)
Let v0, g ∈ Rn, G ∈ Rn×n be given and assume ρ(G) < 1. Consider the recurrence vk+1 = Gvk + g.
Show by induction that for k > 0
and therefore that Hence, show that
(b)
(c)
(d) (e) (f)
k→∞
Write down explicit Euler’s method for the system (1) in the form
vk=Gkv0+ Gj g, j=0
vk =Gkv0+(I−G)−1(I−Gk)g. lim vk = (I − G)−1g.
Uk+1 = Uk + htF(Uk), (3) where ht is some constant time-step and F is a function which you should specify in terms of K,M and f.
Write recurrence (3) in the form (2), 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 the eigenvalues of G are real and positive.
Show that if 0 < ht < 2/λ1(M−1K) then ρ(G) < 1.
Write a script file eulertest.m to verify the result in part (c). In particular, your file should
• generate the matrices K, M and f using femsol given the data
Ω=(0,1)2, f(x)=1, g(x)=0, ε=10−2, c(x)=1+x+y;
• use u0 = 1;
• estimate λ1(M−1K) using the power method and find a suitable value of ht;
• use Euler’s method (3) to compute an approximation U ̃ to u(x, T ), where T = 5 on three subdivisions
of Ω corresponding to fem.level equal to 4, 5 and 6;
• evaluate and display the norm ∥U ̃ − U∞∥2 for each of the three subdivisions above.
k−1
(2)