Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: daixieit

MATH0033

Answer all questions .

1.        (a) Consider the function f : [0, ∞) → R defined by f(x) = sinx − x cos x.

We are interested in finding positive roots of f . Consider the following two fixed point iterations

xk+1  = ϕ 1 (xk ),                        ϕ1 (x) = tanx

k+1  = ϕ2 (k ),                        ϕ2 (x) = x  +

where x0 and 0 are initial guesses. For each of the two fixed point iterations above, determine

(i) if the method is consistent for any positive root α > 0 of f;           (ii) if the method is locally convergent near a positive root α > 0 of f .

If either of these methods is locally convergent, determine also the corre- sponding order of convergence of that method.

(b) Provide an original example of a real-valued function f of a single real vari- able and a corresponding real root α of f such that Newton’s method has local convergence of order equal to three, and no higher than three. Justify your answer by proving that your example fulfills the required condition.  Your  answer  should  consist  of an  original  example  of your  own  creation that is not found already in the course materials or other references . (25 marks)

2.  Let A Rn ×n  be the symmetric tridiagonal matrix given by                                                           

A =           . . .     . . .     . . .                                        

(               )

where a and b are both positive real numbers, and where all other entries of A not shown above are zero.

(a) Show that the eigenvalues λk  of A are given by

λk  = a + 2bcos ( ) ,

for k = 1, . . . ,n. What is the maximum and minimum eigenvalue of A? Hint:  check that the vectors vk  given by

vj(k)  = sin ( )    j = 1, . . . ,n

are  eigenvectors for k = 1, . . . ,n .

(b) Considering still the matrix A given above, show that the Jacobi method for the system Ax = b is convergent if and only if A is positive definite.

(c) For a given initial guess x0 , let xk  , k ≥ 1, denote the k-th iterate of Jacobi’s method for the linear system Ax = b. Let ϵ > 0 be a given tolerance. Find

a sufficient condition on k in terms of n, a, b and ϵ to ensure that ∥xk  − x∥2  ≤ ϵ .

Use this to describe the behaviour of the convergence of Jacobi’s method as n becomes large.

Hint:  your answer should distinguish  the  cases  where a > 2b, a = 2b and a < 2b . (25 marks)

3.  Let A  ∈ Rn ×n   be a symmetric positive definite matrix, let b  Rn .   For a given 北0  ∈ Rn , consider the Conjugate Gradient method given by the sequence of

iterations

k+1  = 北k + αkpk ,       αk  =

pk+1  = rk+1 βkpk ,   βk  = (Apk , rk+1)

(rk , rj ) = 0   Aj = 0, . . . ,k − 1.

Recall any results from the lectures that you use as part of your answer. (b) Show that the coefficients αk  can be equivalently written as

rk 2(2)

αk  =

What can be said about the error  − 北k  if αk  = 0 for some value of k?

(c) Construct an original example of a matrix A and vector b to show why the Conjugate Gradient method cannot generally be applied in the case where A is not symmetric positive definite.   Your answer should show clearly the reason why the Conjugate Gradient method is not applicable to your example.

Your answer should consist an  original example  of your own  creation that is not found already in the course materials or other references . (25 marks)

y\ (t) = y2 (t)   t > 0,

y(0) = 1

and consider it’s approximation by the Backward Euler method with step length h given by

(1)                                    un+1  = un + hun(2)+1 ,    u0  = 1.

(a) Show that the equation (1) admits two real solutions for un+1  in terms of un , provided that un  satisfies a condition that you should specify as part of your answer.  Explain why one of these solutions is not appropriate for approximating the solution y(t) and can thus be discarded.

(b) Considering only the appropriate solution found in part (a), find a recur- rence relation for the quantity zn  = 1 − 4hun  of the form zn+1  = ϕ(zn ) for some function ϕ that you should determine.

(c) Using your answers in part (a) and (b), prove that for any h > 0, there exists an n ∈ N such that no real solution un+1  exists. What feature of the solution y(t) of the original problem does this numerical scheme reproduce? Hint:   you  may  use  without proof the fact  that  the  exact  solution  of the initial value problem is y(t) = (1 − t)−1 (25 marks)

MATH0033 Numerical Methods

Solutions Exam 2021-2022

1.   (a)   (i) Suppose α > 0 is such that sin α α cos α = 0.  Then it follows from implies that sin α    0.   Hence we can divide by cos α  and find that α = tan α = ϕ 1 (α). So the first fixed point iteration is consistent. Since α > 0 is nonzero we also clearly get  =  and thus α = ϕ2 (α), so the

second method is also consistent.

(ii) We compute ϕ(α) to find

ϕ(α) = 1 + tan2 α

From above we have sin α  0 and cos α  0 for a root, and thus ϕ(α) > 1 for any positive root α . Therefore the first method is not locally conver-

gent (in any neighbourhood) of any positive root α .

We next compute ϕ(α)

 

ϕ(α) = 1 +      −                  =      −            = 0.

Since |ϕ(α)| < 1 the second method is locally convergent, and moreover since ϕ is smooth in a neighbourhood of any positive root α, the conver- gence is at least quadratic. We check that it is not higher than quadratic since

kl   =  =  (  + ) =   0

owing to sin α  0 and cos α  0 from above.

(b)  This question requires an original example for an unseen situation, for which

there is no unique  correct solution.  Marking will take into  account the  origi- nality and interest of the example, and also the justifications provided with the example .  In particular, the answer must pay particular attention to justifying why the  order of convergence is  exactly  equal to  but no higher than three .  A model solution is given below.

Consider f(x) = x x3 /3, which has a root α = 0. Newton’s method for this

x x3 /3         2x3      

the lectures that Newton’s method has local convergence of order at least 2.

However, in fact the order of convergence is exactly three since, supposing that x0   0 is sufficiently close to zero (e.g |x0 | < 1/^3 but x0   0, we get xk   0 for all k by induction, and

xk+1 α        xk+1          2     1           2

(xk  α)3           xk(3)           3 1 − xk(2)         3

since xk  → 0 as k → ∞ by local convergence. Therefore the order of conver- gence is at least three.  To see if it is at most three, note that any order of convergence q > 3 would require necessarily that

kl   = kl  (xk α)q3 = 0.

However in this example, limk→∞     0, so the order of convergence

cannot be higher than three.

2.   (a) We can extend the formula of vj(k)  for j = 0 and j = n + 1 to obtain v0(k)  = vn(k)+1 = 0 for all k = 1, . . . ,n. If we set ωk  = , we can write

(Avk )j  = avj(k) + b(vj(k)1 + vj(k)+1)

= asin(ωkj) + b(sin(ωk(j − 1)) + sin(ωk(j + 1)) = asin(ωkj) + 2bcos(ωk)sin(ωkj)

= (a + 2bcos(ωk))vj(k)  = λkvj(k)

for all j = 1, . . . ,n, where we note that the special cases j = 1 and j = n use the fact v0(k)  = vn(k)+1  = 0 as shown above. Therefore we find that A has n distinct eigenvalues {λk}k(n)=1 .

If a and b are both positive, then it is easy to see that the maximum eigenvalue

of A is

λmax = a + 2bcos ( )

whereas the minimum eigenvalue is

λmin = a + 2bcos ( ) = a 2bcos ( ) .

(b) The iteration matrix of the Jacobi method is B  = I A, and thus the

 

µk  = 1  =  cos ( )

The Jacobi method is convergent if and only if ρ(B) < 1, i.e. |µk| < 1 for all k, or equivalently

1 <  cos ( ) < 1  ⇐⇒ 0 < a ± 2bcos ( ) .

Since

λn+1k  = a 2bcos ( ) ,    λk  = a + 2bcos ( )

we see that |µk| < 1 for all k iff all eigenvalues of A are positive, i.e. that A is s.p.d. Therefore the Jacobi method is convergent with ρ(B) < 1 if and only if A is s.p.d.

(c) Since B is symmetric, the iteration xk  satisfy the bound xk − x∥2  ≤ ρ(B)k ∥x0 − x∥2 ,

according to Corollary 4.6.3 in the lectures.  Provided that ρ(B) ∈ (0, 1), a

sufficient condition to obtain xx0 2  ϵ is then

ρ(B)k  ≤  ⇐⇒ k ≥ log ϵ o(l)g(o)B(x) x0 2

We compute

ρ(B) = mk(a)x|µk| =  k a,x.,n 'cos ( ) ' =  cos ( )

Therefore, the sufficient condition above is equivalent to

log ϵ  logx x0 2

|log (  cos( n   1 )) |

Note that

2b

n→∞                   a .

If a > 2b then the spectral radius remains uniformly bounded away from 1 for all n and log(ρ(B)) → log( ) which remains bounded away from 0 for all n.  Hence the sufficient condition to ensure that the error reaches a desired tolerance can be taken to be independent of n.  However if a = 2b then the spectral radius ρ(B) < 1 for all n, but tends to 1 as n becomes large, and the right-hand side in the condition above tends to + ∞ as n → ∞.  This shows that we can expect the convergence to become very slow for large n. If a < 2b then there would then exist n sufficiently large such that A is not s.p.d. — in this case the condition above is no longer applicable as the Jacobi method is no longer convergent.

3.   (a) From the lectures, it is shown that rk  is orthogonal in the Euclidean inner-

product to the k-th Krylov subspace Vk . It is also shown that rk  ∈ Vk+1  and

that the spaces are nested i.e. V1  V2  ⊂ · · · ⊂ Vk .  Therefore, if j k 1,

(b) From the lectures, we have rk   ⊥ Vk  in the Euclidean inner-product, and pk  ∈ Vk+1  for all k ≥ 0.  We have also seen in the lectures that since A is SPD, that ∥pk ∥A(2) = (Apk , pk ). Next, writing

pk  = rβk1p1

we use the linearity of the inner-product

(rk , p ) = (rkk , r  )k βk1(rk , p1) = (rk , r  )k

where the term (rk , pk1) = 0 since rk  ⊥ p1  ∈ Vk .  Then, we use the fact that (rk , rk ) = ∥rk ∥2(2)  to find that

rk 2(2)

αk  =

We have αk  = 0 for some k if and only if rk  = 0 and thus the residual is exactly zero, which is equivalent to xk  being equal to the exact solution x.

(c)   The answer should include an original example not seen elsewhere .  The grade should reflect the originality and interest of the example, i. e . an overly trivial example  (e .g.   A  being  the  zero  matrix)  does  not get full marks .   The  grade should  also  reflect  the  quality  of the justification  as  to  why  the  Conjugate Gradient  method  is  not  applicable:  for full  marks  the  answer  should  detail where  and how the  algorithm fails .  A reasonable  example for instance would be as follows:  Consider

A = (0(1)   1(2)) , b = ( 1) , x0 = 0

The matrix A is nonsingular, but it is not symmetric positive definite.  In this case, we get p0  = r0  = b.  But then α0  is not well-defined since in the denominator we have

(Ap0 , p0 ) = b = ( 1   1) (0(1)   1(2))( 1) = ( 1   1) ( 1(1)) = 0.

This means that the first iterate x1 of the Conjugate Gradient method is not well-defined.

Other suitable examples could include symmetric but indefinite matrices, etc .

4.   (a) The BE method gives

un+1 = un + hun(2)+1   ⇐⇒ (un+1 )2 = 1  

which has real solutions

un+1 =  (1 ± ^1 4hun)

if and only if 4hun  ≤ 1 or equivalently un  ≤  .

To find which of the roots of the quadratic above is the physically relevant one,  notice that if we choose the positive root at the first timestep,  i.e.

u1   =   ( 1 + 1 4h), then we would have u1   → ∞ as h  0, whereas

is not appropriate for approximating y(t) and can be discarded (the same argument applies to any later timestep as well).  On the other hand, for the second root we get by Taylor expansion

un+1 =  (1 ^1 4hun) (2hun+ 2h2un(2)) = un + hun(2)

as h → 0, so this approaches the Forward Euler solution. We conclude that the physically relevant solution at each timestep is then

un+1 =  (1 −^1 4hun)

(b) Let zn = 1 4hun, therefore

zn+1 = 1 − 4hun+1 = 1 − 2(1 − zn) = 2zn − 1,

so we have the recurrence relation zn+1 = ϕ(zn) with the function zn  defined

by ϕ(z) = 21.

(c)  Recall from (a) that a real solution un+1  exists if and only if zn  ≥ 0. We now

show by contradiction that there is an n ∈ N such that zn  < 0.  Suppose to

the contrary that zn  0 for all n N. Next, note that z0 = 1 4h < 1. Note

implies that the sequence zn+1 = ϕ(zn) ≤ zn  for all n ∈ N, so the sequence is

decreasing.

Then the sequence zn  is decreasing and bounded from below, so it must have

a limit 0. Moreover since z0  = 1 4h < 1, we must have 1 4h < 1.

point ℓ < 1.  Therefore, we conclude that there must exist an n ∈ N such that zn   <  0 which is equivalent to un   >    thus the numerical solutions blow up.  This is similar to the behaviour of the continuous solution since

y(t) = (1 t)1  blows up as t 1.