MA30051/MA50178 Assessed Coursework 2020/21
Set: 4p.m. Thursday 19th November 2020.
Due: 12 Noon Thursday 3rd December 2020. Submit a single zip file of all of your work
to the Coursework section of the moodle site.
Expected completion time: About 12 hours in total, including writing–up.
This assignment is worth 25% of the mark for MA30051 and 20% for MA50178. Figures
in square brackets, e.g. [2], represent the marks available for each question, out of a total
of 25.
Remarks
(i) Any form of cheating is a serious breach of University regulations and will result in
disciplinary proceedings. You should not discuss the details of your work with anyone else.
The work which you hand in must be your own. You should be prepared to explain to
an examiner anything which you write, if asked to do so. In particular, if it is discovered
that all or part of your code or your written work has been copied, both parties involved
risk a penalty to be determined by the Dean.
(ii) Results from MA30051, as well as Year 1 and Year 2 lecture courses, may be quoted
without proof. Results from other sources must be stated and proved in full.
(iii) Partial answers to questions will receive partial credit, particularly if they are well
justified.
(iv) Codes relating to Questions 2, 3 and 4 are available in the Coursework section of the
course moodle site.
(v) If exceptional circumstances, e.g. disability or illness, prevent you from completing
this assignment by the deadline, contact the lecturer or the D.o.S. team at an early
stage. E-mail masath or [email protected].
1
Accelerated convergence for the Power method
The diffusion Dirichlet problem: Let d : [0, 1] −→ (0, ∞) be a function describing
the diffusion coefficient, as it varies over the domain [0, 1]. The corresponding diffusion
problem is to find smooth u : [0, 1] −→ R such that
−d(y)u′′(y) = f(y), u(0) = u(1) = 0, (1)
for a given function f : [0, 1] −→ R.
The discretisation matrix A: As in Example 1.2.3, problem (1) is discretised on a grid
yj = jh, 0 ≤ j ≤ N + 1, h = 1/(N + 1), to yield the matrix–vector equation
Ax = b, bj = h
2f(yj), 1 ≤ j ≤ N. (2)
Due to the coefficient d(y), the discretisation matrix A ∈ RN×N differs from (1.14):
A =

2d1 −d1 0 · · · · · · · · · 0
−d2 2d2 −d2 0 · · · · · · 0
0 −d3 2d3 −d3 0 . . . 0
...
. . . . . . . . . . . . . . .
...
0 · · · 0 −dN−2 2dN−2 −dN−2 0
0 · · · · · · 0 −dN−1 2dN−1 −dN−1
0 · · · · · · · · · 0 −dN 2dN

, (3)
where dj = d(yj) > 0, 1 ≤ j ≤ N .
This coursework is concerned with finding λ1, the largest eigenvalue of A.
I Approximate eigenvalue location: Initially, it is not clear that σ(A) ⊂ R, as A
is not symmetric if the dj are not equal. The first question uses theory to locate the
spectrum of A approximately.
Q1(a) Use Gerschgorin’s theorem to prove that if λ ∈ σ(A), then∣∣λ− 2‖d‖∞∣∣ ≤ 2‖d‖∞, d := [d1, . . . , dN ]T ∈ RN . [2 marks]
(b) Prove that 0 6∈σ(A). [1 mark]
(c) Prove that A is similar to a symmetric tridiagonal matrix. [3 marks]
(d) Deduce that A has a full set of eigenpairs {λj, vj}Nj=1 ∈ R × RN \ {0}, such that
σ(A) ⊂ (0, 4‖d‖∞]. [1 mark]
Related fact: Analysis of the characteristic equations of the principal submatrices of A
shows that the eigenvalues of A are distinct, and therefore strictly ordered. This fact is
assumed in the sequel.
2
II The Power method
Consider the matrix (3) for the case N = 100 and d = [1, 1, . . . , 1]T ∈ RN .
Q2(a) Run the code Power1.m to implement the Power method for this matrix.
Determine the approximately straight-line region of the graph of log10 ‖rk‖2 versus
k, and use suitably-spaced values in that region to estimate α ∈ (0, 1) such that
‖rk+1‖2/‖rk‖2 ≈ α. Write down α to 11 s.f. (significant figures). [1 mark]
(b) State a theoretical formula for α, valid for large k. (Theoretical justification of this
answer in not required.) Evaluate this formula with the aid of the MATLAB instruction
eig(A). Briefly comment on the difference between your answer and the computational
α found in part (a). [1 mark]
(c) Find the number of iterations required to reduce ‖r0‖2 by a factor of 108, under the
idealised assumption that ‖rk‖2 = ‖r0‖2αk, k ∈ N, for the theoretical value of α found
in part (b). Justify your working mathematically. [1 mark]
III The Two-step Power method: Here, we continue to consider the matrix A ∈ RN×N
given by (3) with d = [1, 1, . . . , 1]T ∈ RN for N = 100.
A linear mapping of A: Label the distinct eigenvalues of SPD A ∈ RN×N as
λ1 > b := λ2 > · · · > λN =: a > 0.
Thus, σ(A) \ {λ1} ⊂ [a, b]. For shift s = (a + b)/2 and half-diameter df = (b − a)/2,
consider the map
λ 7→ λ̂, λ̂ := (λ− s)/df,
which maps [a, b] to [−1, 1]. Let  := (A− sI)/df . Then, σ(Â) = {λ̂j}Nj=1, such that
λ̂1 > 1 = λ̂2 > · · · > λ̂N = −1.
Thus, σ(Â) \ {λ̂1} ⊂ [−1, 1] and λ̂1 > 1.
Pseudocode for the two-step power method:
Initialisation: Choose random z0, x1 ∈ RN \ {0} where ‖x1‖2 = 1 and tol > 0.
Set y1 = Âx1, µ̂1 = x
T
1 y1, r1 = y1 − µ̂1x1 and evaluate ‖r1‖2.
Iteration: For k = 2, 3, . . .
wk = 2yk−1 − zk−2,
xk = wk/‖wk‖2,
yk = Âxk,
µ̂k = x
T
k yk,
rk = yk − µ̂kxk; Stop if ‖rk‖2 ≤ tol,
zk−1 = xk−1/‖wk‖2,
end of for loop.
3
Outputs: µ = dfµ̂k0 + s ≈ λ1, where µ̂k0 is the final iterate, and a graph of log10 ‖rk‖2
versus k.
Q3(a) Show a relationship between the two-step power method algorithm described above
and the standard power method for the matrix
B =
[
0 I
−I 2Â
]
∈ R2N×2N . [2 marks]
(b) Show that ζ ∈ σ(B) if and only if ζ ∈ C is the solution of a quadratic
ζ2 − 2λ̂ζ + 1 = 0,
for some λ̂ ∈ σ(Â). [1 mark]
(c) Show that there exist ζ1, ζ2 ∈ σ(B) such that
ζ1 > 1 > ζ2 > 0.
Furthermore, show that |ζ| = 1 for ζ ∈ σ(B) \ {ζ1, ζ2}. [2 marks]
(d) Write down, without justification, a theoretical value for lim
k→∞
‖rk‖2
‖rk−1‖2 . (For this
part of the question it is assumed that the two-step power method is implemented
without stopping criterion to allow k −→∞.) [1 mark]
(e) Run the code twostep.m, which implements the two-step power method. Observe
that the computed estimate of ‖rk‖2/‖rk−1‖2 is larger than the theoretical asymptotic
value found in part (d); i.e. the convergence is slower in practice. Give an extended
mathematical analysis of the behaviour of the two-step power method sufficient to
explain why this is the case. [4 marks]
IV A relaxation of the Two-step Power method: Here, we continue to consider the
matrix A ∈ RN×N given by (3) with d = [1, 1, . . . , 1]T ∈ RN for N = 100.
Let n ∈ N, n > pi/ cos−1(λ̂−11 ), be fixed in advance of implementing the algorithm. In a
modified version of the two-step power method, the pseudocode line wk = 2yk−1− zk−2 is
changed to
wk = 2 cos(pi/n)yk−1 − zk−2.
The modified pseudocode is implemented in the MATLAB code twostepr.m.
Q4(a) Find, with brief justification, a theoretical value for lim
k→∞
‖rk‖2
‖rk−1‖2 . (As in 3(d),
assume that the method is implemented without stopping criterion.) How does your
answer compare to the one in Q3 (d)? [1 mark]
(b) Run the code twostepr.m for n = 389 as well as the code twostep.m. Give an
extended mathematical analysis to explain the two convergence graphs of log10 ‖rk‖2
for k ∈ [150, 389]. [4 marks]
4

欢迎咨询51作业君