1. Let f (x) = log(x + 2) - x/2.

(a)  Show that f has a unique root α in the interval (0, o).

(b) Prove that [2, 4] is a valid starting interval for the bisection method, and determine

the number of iterations required to find α with an absolute error less than or equal to 10-6 .

(c) Write down the Newton iteration explicitly for this problem.  Does the iteration converge linearly or quadratically to α? Justify your answer.

(d)  Show that φ1 (x) = 2 log(x + 2) and φ2 (x) = eα/2 - 2 both provide consistent fixed point formulations for f (x) = 0.

Prove that one of these formulations converges to α for any initial guess in [2, 4], while the other does not converge to α, regardless of how good the initial guess is, unless one of the iterates equals α itself.

2.   (a)  Consider the iterative solution of the linear system Ax = b, where A = ╱ 3(ξ)   4(1) 、,

ξ  0 is a real parameter, and x, b e R2 .

Without computing eigenvalues, give a sucient condition on ξ for the Jacobi and Gauss-Seidel methods to be convergent.

(b) Write down the iteration matrices for the Jacobi and the Gauss-Seidel methods

for the matrix in part (a), and determine a necessary and sucient condition on ξ for the Jacobi and Gauss-Seidel methods to be convergent.

(c) In the gradient method for the solution of a general linear system Ax = b, the relative step length αk  in the increment xk+1  - xk  = αk rk   (with rk  = b - Axk denoting the residual) is chosen so as to minimise the magnitude of the error vector ek+1  = x - xk+1  with respect to a certain A-weighted norm. Explain why it is not in general possible in a practical implementation to choose αk  so as to minimise the magnitude of ek+1  with respect to the Euclidean norm.

3.  (a)  Show that if X, Y e Rnxn  are symmetric then ρ(XY) s ρ(X)ρ(Y).

(b) Let A e Rnxn  be an invertible matrix, and let b e Rn . Suppose that A = A1 + A2 for some symmetric positive definite matrices A1 , A2  e Rnxn .  Given a1 , a2  > 0,

consider the following iteration, where I e Rnxn  denotes the identity matrix:

 }       k = 0, 1, 2, . . . .

Show that the iteration can be written in the form xk+1   = Bxk  + c for some B e Rnxn  and c e Rn  that you should specify.

(c)  Show that if

R := λ(1σ(a)1 )  │ λ(2σ(a)2 )  │ < 1

then the iteration converges to the exact solution of the linear system Ax = b.   You may assume without proof that the iteration is consistent, and that if X, Y e Rnxn  then ρ(XY) = ρ(YX).

(d) Now suppose that σ(A1 ), σ(A2 ) c [λ- , λ+] for some 0 < λ-  s λ+  < o.

Show that if a1  = a2  = a > 0 the quantity R defined in (c) satisfies

R s ╱max {  │ ,  │}、 2 .

By minimising the right-hand side of this inequality as a function of a, show that a can be chosen so as to give

|xk+1 - x|2  s \2 |xk  - x|2 ,        k = 0, 1, 2, . . . ,

assuming the matrix B in (b) is symmetric.

4.  Consider the following IVP involving a system of ODEs for two unknowns S and I:

,

ì ì ì ì ì

(

ì ì ì ì ì

(

dS

(t) = -βI(t)S(t),

dt

dI

(t) = βI(t)S(t) - γI(t), dt

S(0) = S0 ,        I(0) = I0 ,

0 < t < T,

0 < t < T,

where T, β , γ , S0  and I0  are all positive real numbers.

(a) Write the system in the form

, y\ (t) = f(t, y(t)),    0 < t < T,

defining y, y0  and f, and show that if u, v e [0, 1] × [0, 1] then

|f(t, u) - f(t, v)|o  s L|u - v|o

for some constant L > 0 that you should specify.

(b) Define the truncation error Tn  for the backward Euler scheme applied to this IVP

system with a uniform time step, and show that, under appropriate smoothness assumptions,

|Tn |o  s Ch,

where h is the time step length and C > 0 is a constant you should specify. State the smoothness assumptions on S(t) and I(t) under which your analysis is valid.

(c)  Show that the backward Euler solution un  considered in part (b) satisfies |y(tn ) - un |o  s C* h,        0 < h < h* ,

for some C*  > 0 and h*  > 0 you should specify.

(d) Write down the system (S) of nonlinear equations that must be solved at each time step in the backward Euler method for this IVP system.

Show that if the system (S) has a solution in [0, 1] × [0, 1] and if γ > β then the Newton method for the solution of (S) converges quadratically to that solution, assuming a sufficiently good initial guess (that you need not specify).  If you use any results from lectures, quote them carefully.

1.  (a) The function f is continuous on [0, ∞) and satisfies f(0) = log 2 > 0 and f(x) →

f\ (x) = 1/(x + 2) − 1/2, which is negative on (0, ∞), so f is strictly decreasing.

a point between them where f\  vanishes, which is impossible.

(b)  [2, 4] is a valid starting interval for the bisection method because f(2) = log 4− 1 >

0 and f(4) = log 6 − 2 < 0 (since e < 4 and e2  > 6), so f(2)f(4) < 0.

|xk  − α| ≤ (b − a)/2k+1  = 2/2k+1  = 2 −k .

Hence to ensure the error is less than or equal to 10 −6  it suffices to ensure that 2 −k  ≤ 10 −6, i.e.

log 10

(c) The Newton iteration is

xk+1  = xk  −  = xk  − l1(o)/(g)x(x)k(k) 2(2) 1(x)/(k)2(/)2 ,        k = 0, 1, . . . .

which is non-zero since α > 0.

(d) Both formulations are consistent with f(x) = 0 because for x > 0

2log(x + 2) = x ⇔ log(x + 2) = x/2 ⇔ x + 2 = ex/2  ⇔ x = ex/2 − 2.

The fixed point iteration xk+1   = ϕ 1 (xk ) converges to α for any initial guess in [2, 4] by the contraction mapping theorem.  To justify this we note that ϕ(x) = 2/(x + 2) ≤ 2/(2 + 2) = 1/2 < 1 for all x ∈ [2, 4], so ϕ 1  is a strict contraction on [2, 4].  Also, ϕ 1 (x) ∈ [2, 4] for x ∈ [2, 4] since ϕ(x) > 0 for x ∈ [2, 4], so ϕ 1  is strictly increasing on [2, 4], and ϕ 1 (2) = 2log 4 > 2, ϕ 1 (4) = 2log 6 < 4.                The fixed point iteration xk+1   =  ϕ2 (xk ) does not converge to α,  even for an arbitrarily good initial guess, unless one of the iterates equals α itself.  To prove this we note that ϕ2   is differentiable with ϕ(x)  = ex/2/2.   From part  (b), or from the analysis of the ϕ 1  iteration, we know that α ≥ 2, which implies that ϕ(α) = eα/2/2 ≥ e/2 > 1.  Moreover, by the continuity of ϕ there exists δ > 0

such that ϕ(x)  >  1 for all x  ∈ I6   :=  (α − δ,α + δ).   Suppose that xk   → α .

2.   (a) The Jacobi and Gauss-Seidel methods are convergent provided A is strictly di- agonally dominant by rows. This occurs when |ξ| > 1.

(b) Let

D = ╱ 0(ξ)   4(0) 、,    E = ╱ 3   0(0) 、,    F = ╱ 0(0)   1 、.

Then the iteration matrices for the Jacobi and Gauss-Seidel methods are

BJ  = D−1(E + F) = ╱ 10(/ξ)   1/(0)4 、╱ 3   1 、= ╱ −3(0)/4   − 0(1)/ξ 、

and

BGS  = (D − E)−1F = ╱ 3(ξ)   4(0) 、−1  ╱ 0(0)   1 、    =  ╱ 3   ξ(0) 、╱ 0(0)   1 、 =  ╱ 0(0)   4 、

= ╱ 0(0)

 、.

To obtain necessary and sufficient conditions for convergence we determine when the spectral radii of the iteration matrices are smaller than one.

For BJ  the eigenvalues satisfy

λ2 −  = 0,

so that ρ(BJ ) = !  .

It follows that the Jacobi method converges if and only if |ξ| > 3/4.

For Gauss-Seidel the eigenvalues satisfy

λ ╱λ − 4ξ(3)、 = 0,

so that ρ(BGS ) = 4| .

Once again the condition for convergence is |ξ| > 3/4.

(c) We are considering non-stationary iterations of the form

xk+1 − rk  = αk rk ,        k = 0, 1, 2, . . . .

where rk  = b − Axk .

We recall that the residual and the error ek  = x − xk  are related by rk  = Aek ,        k = 0, 1, 2, . . . ,

and that the error updates by

ek+1 − ek  = −αk rk ,        k = 0, 1, 2, . . . .

In the gradient method, given xk  we choose αk  so as to minimise ∥ek+1∥A , where ∥v∥A(2)   =  (Av, v) for v  ∈ Rn .   In the current question we are considering the possibility of instead choosing αk  so as to minimise ∥ek+1∥2 . Note that

∥ek+1 ∥2(2)  = (ek+1 , ek+1)

= (ek  − αk rk , ek  − αk rk )

= ∥ek ∥ − 2αk (rk , ek ) + α ∥rk ∥ .

This is a quadratic in αk , and by differentiating with respect to αk  and setting the resulting expression equal to zero, we find that the unique minimum is at

(rk , ek )       ∥ek ∥A(2)

αk  =   ∥rk ∥2(2)    =  ∥rk ∥2(2) .

However, the problem is that, because ek   = x − xk  is unknown (because x is

3.  (a) Using the fact that ρ(A) ≤ ∥A∥M   and  ∥AB∥M   ≤ ∥A∥M ∥B∥M  for any induced matrix norm ∥ · ∥M , and using the equality between the spectral radius and the 2-norm for symmetric matrices,

ρ(XY) ≤ ∥XY ∥2  ≤ ∥X∥2 ∥Y ∥2  = ρ(X)ρ(Y).

(b) By eliminating xk+1/2  we can write the iteration in the desired form with

B = (I + a2 A2 )−1(I − a2 A1 )(I + a1 A1 )−1(I − a1 A2 ), c = (I + a2 A2 )−1  ╱a1 (I − a2 A1 )(I + a1 A1 )−1 + a2 I、b,

where we used the fact that, for both i = 1, 2, I+ai Ai is SPD, and hence invertible, because Ai  is SPD, ai  > 0 and σ(I + ai Ai ) = 1 + ai σ(Ai ).

(c) A sufficient condition for convergence of the iteration is that ρ(B) < 1. The stated result then follows from the fact that

ρ(B) ≤ λ(1σ(a)1 ) │ 1(1)  │ λ(2σ(a)2 ) │ 1(1)  │ .

To prove (*), we first use part (a) and the hint to see that

ρ(B) = ρ ╱(I − a2 A1 )(I + a1 A1 )−1(I − a1 A2 )(I + a2 A2 )−1、     ≤ ρ ╱(I − a2 A1 )(I + a1 A1 )−1、ρ ╱(I − a1 A2 )(I + a2 A2 )−1、.

(*)

(**)

Hence F(aλ) ≤ G(a) := max{F(aλ− ),F(aλ+ )}, which proves the first inequality in (d).

We now want to minimise G(a) as a function of a. If λ− = λ+  then the minimiser is clearly at a = 1/λ− = 1/λ+ .  If λ−  λ+ , then on (0, 1/λ+ ) both F(aλ+ ) and F(aλ− ) are decreasing functions of a with F(aλ+ ) ≤ F(aλ− ), and on (1/λ− , ∞) both F(aλ+ ) and F(aλ− ) are increasing functions of a with F(aλ− ) ≤ F(aλ+ ). Hence the value of a minimising G(a) lies in [1/λ+ , 1/λ−], where

G(a) = max {1(1)  − 1(1)  } .

The minimium occurs where

1 − aλ−         1 − aλ+

which after manipulation gives a2 λ+ λ− = 1, i.e.

a = a∗ := 1/^λ+ λ− ,

with

1 − !λ−/λ+

1 + !λ−/λ+ .

Combining this result with the first inequality in (d), and recalling that ρ(B) ≤ R, and assuming that B is symmetric, we find that

∥B∥2  = ρ(B) ≤ │  1(1)  !(!)\2 ,

from which the claimed linear convergence estimate follows by a result from lec- tures.

4.  (a) The system can be written in the form

y(0) = y0 ,

with y(t) = (y1 (t),y2 (t)) = (S(t),I(t)), y0  = (S0 ,I0 ) and

f(t,y) = (−βy2y1 ,βy2y1 − γy2 ) = ( −βIS,βIS − γI).

For the claimed Lipschitz continuity, students should be given credit for any valid analysis. For instance, noting that

|u2u1 −v2v1 | = |u2u1 −u1v2 +u1v2 −v2v1 | ≤ |u1 ||u2 −v2 |+|v2 ||u1 −v1 | ≤ 2∥u−v∥∞ ,

where we used the fact that |u1 |, |v2 | ≤ 1, we find that

∥f(t,u) − f(t,v)∥∞  = max{|β(u2u1 − v2v1 )|, |β(u2u1 − v2v1 ) − γ(u2 − v2 )|} ≤ max{2β,2β + γ}∥u − v∥∞

= (2β + γ)∥u − v∥∞ ,

which proves the required estimate with L = 2β + γ .

Alternatively, students might apply the lemma from lectures (based on the MVT) that gives

∥f(t,u) − f(t,v)∥∞  ≤ max ∥Dy f(t,(1 − ξ)u + ξv)∥∞ ∥u − v∥∞

where

Dy f(t,w) = ╱ 2     ββ11γ 、

denotes the partial Jacobian of f with respect to y. Assuming that w ∈ [0, 1] × [0, 1]

we have −γ ≤ βw1  − γ ≤ β, and then, using a formula for the matrix ∞-norm

∥Dy f(t,w)∥∞  = max{|βw2 | + |βw1 |, |βw2 | + |βw1 − γ|}

≤ max{2β,β + max{β,γ}},

= max{2β,β + γ},

which proves the required estimate with L = max{2β,β + γ}.  (This is slightly sharper than the estimate obtained above.)

(b) We define the (vector-valued) truncation error as

Tn  =  (yn+1 − (yn + hf(tn+1 , yn+1)) ,

and note that, by the ODE system,

Tn  =  ╱yn+1 −╱yn + hy+1、、,

where yn  = y(tn ) and y = y\ (tn ) etc., with y denoting the exact solution. Taylor-expanding each component separately gives

yn+1  = yn + hy +  ╱ y(y) 、

and

y+1  = y + h ╱ y(y)ξ(ξ) 、,

for some ξ1 ,ξ2 ,ξ ,ξ between tn  and tn+1 .

Hence

Tn  = h ╱ y(y) ξ(ξ)2(1)  y(y) ξ(ξ) 、,

so that

∥Tn ∥∞ ≤  max{ξ [0(u)] |y(ξ)|, ξ [0(u)] |y(ξ)|} =  ξ [0(u)] ∥y\\ (ξ)∥∞ ,

proving the desired estimate with C =  supξ∈[0,T] ∥y\\ (ξ)∥∞ .

This analysis is valid provided that the exact solution y = (S,I) is twice differen- tiable, with y\\ = (S\\ ,I\\) bounded on [0,T].

(c) Let wn  = yn − un , where un  denotes the numerical solution at time step tn  = nh.

wn+1  = wn + h(f(tn+1 , yn+1) − f(tn+1 , un+1)) + hTn+1 ,

so that by the Lipschitz property in (a) and the triangle inequality ∥wn+1∥∞ ≤ ∥wn ∥∞ + hL∥wn+1∥∞ + h∥Tn+1∥∞ ,

which, after rearranging and applying the truncation error estimate in (b), gives

∥wn+1∥∞ ≤ a∥wn ∥∞ + b,        0 < h < 1/L,

where a = 1/(1 − hL) = 1 + hL/(1 − hL) and b = Ch2 /(1 − hL).

By induction on n it follows that

∥wn ∥∞ ≤ an ∥w0 ∥∞ + 11 b = (an − 1)  ,

0 < h < 1/L,

where we used the fact that u0  = y0 , so w0  = 0.

Finally, noting that 1 + x ≤ ex  for x > 0, we obtain that

∥wn ∥∞ ≤ ╱enhL/(1−hL) − 1、 ,

0 < h < 1/L,

which proves the desired estimate with e.g.

h∗ = 1/(2L)     and     C∗ = ╱e2TL − 1、 .