[SOLVED] Hive matlab assembly Exercise sheet 4 NLAA Numerical Linear Algebra with Applications

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 RNN , 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(n2)(n2) 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 allx . 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 becomesKu = A+ M.Show that1k(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 topreconditioning the system Ku = f with K.(d) Let = 102, c = 1. Download from Canvas the file nrdiff.mat containing three FEM discretizations Ku = fof the problem in part (c). Write a matlab script part-d.m to solve these three linear systems by using thePreconditioned Steepest Descent method with preconditioner K = A+ cM . You should modify psd.m so thatthe preconditioning step is implemented using the fast solver fpsq.m derived in part (b). In each case writedown 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 108 and shift . Modify the script shiftinvert.mto 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 theabsolute errors in each case, together with an estimate of the order of convergence for the sequence ofapproximations.[50]