MA3227 Numerical Analysis 2 AY 2009/2010 Sem 1 NATIONAL UNIVERSITY OF SINGAPORE MATHEMATICS SOCIETY PAST YEAR PAPER SOLUTIONS with credits to Wang Xingyin MA3227 Numerical Analysis 2 AY 2009/2010 Sem 1 Question 1 (a) Let A = L + D + U D: a diagonal matrix with dagonal entries equal to those of A’s L: a lower triangular matrix with 0 as diagonal entries and the same entries below its diagonal as those of A’s U : an upper triangular matrix with 0 as diagonal entries and the same entries above its diagonal as those of A’s For Jacobi iteration: P = D, and for Gauss-Seidel iteration: P = L + D (b) P x(k) = (P − A)x(k−1) + b P x = (P − A)x + b Taking the difference, we have P (x(k) − x) = (P − A)(x(k−1) − x) x(k) − x = (I − P −1 A)(x(k−1) − x) Inductively, x(k) − x = (I − P −1 A)k (x(0) − x) (k) x − x = (I − P −1 A)k (x(0) − x) ≤ (I − P −1 A)k (x(0) − x) Given ρ(I − P −1 A) < 1, limk→∞ (I − P −1 A)k = 0, and (0) x − x is a constant independent of k so limk→∞ x(k) − x = 0. Question 2 (a) Consider backward Euler method to solve y 0 (t) = λy(t) We have, y (n+1) − y (n) = λy (n+1) ∆t y (n) y (n+1) = 1 − λ∆t NUS Math LaTeXify Proj Team Page: 1 of 5 NUS Mathematics Society MA3227 Numerical Analysis 2 AY 2009/2010 Sem 1 Inductively, 1 y (n) = ( )n y (0) 1 − λ∆t n (0) 1 (n) y y = 1 − λ∆t 1 n Given λ < 0, ∆t > 0 ∀n ∈ N , 1−λ∆t < 1 so y (n) ≤ y (n) . (b) Consider forward Euler method to solve y 0 (t) = λy(t) We have, y (n+1) − y (n) = λy (n) ∆t y (n+1) = (1 + λ∆t)y (n) Inductively, y (n) = (1 + λ∆t)n y (0) (n) n y = |(1 + λ∆t)| y (0) Given 0 ≤ ∆t ≤ 2 −λ and λ < 0, so |1 + λ∆t| ≤ 1 and y (n) ≤ y (n) . Question 3 (a) Assuming y(t) ∈ C 2 [0, T ], τn = y(tn+1 ) − y(tn ) − ∆tf (tn , y(tn )) = [y(tn ) + y 0 (tn )∆t + O(∆t2 )] − y(tn ) − ∆ty 0 (tn ) = O(∆t2 ) So |τn | ≤ C∆t2 for some constant C. (b) |en+1 | = |y(tn+1 ) − yn+1 | = |y(tn ) + ∆tf (tn , y(tn )) + τn − (yn + ∆tf (tn , yn )| ≤ |y(tn − yn )| + ∆t |f (tn , y(tn ) − f (tn , yn )| + |τn | ≤ |en | + ∆tL |y(tn ) − yn | + C∆t2 = (1 + ∆tL) |en | + C∆t2 Inductively, for 0 ≤ n ≤ T ∆t , |en | ≤ (1 + ∆tL)n |e0 | + C∆t2 n−1 X (1 + ∆tL)i i=0 (1 + ∆tL)n − 1 = 0 + C∆t2 (1 + ∆tL) − 1 n∆tL e −1 ≤ C∆t2 ∆tL eLT − 1 =C ∆t L NUS Math LaTeXify Proj Team Page: 2 of 5 NUS Mathematics Society MA3227 Numerical Analysis 2 AY 2009/2010 Sem 1 Question 4 (a) kH wk ~ 22 = (H w) ~ T (H w) ~ =w ~ T HT Hw ~ =w ~T w ~ = kwk ~ 22 So kH wk ~ 2 = kwk ~ 2 (b) Let ~v = sign(w1 ) kwk ~ 2 ~e1 + w ~ ~v~v ∗ Then H = I − 2 ∗ . ~v ~v Question 5 (a) Consider y 0 = iλy, y(0) = y0 , then f (tn , yn ) = iλyn . Applying (A), we have, ∆t (iλyn + iλ[yn + ∆t(iλyn )]) 2 ∆t [iλyn + iλyn − ∆tλ2 yn ] = yn + 2 1 = yn + iλ∆tyn − λ2 ∆t2 yn 2 1 = (1 + iλ∆t − λ2 ∆t2 )yn 2 yn+1 = yn + Applying (B), we have, ∆t iλyn ] 2 1 = yn + iλ∆tyn − λ2 ∆t2 yn 2 1 2 2 = (1 + iλ∆t − λ ∆t )yn 2 yn+1 = yn + ∆tiλ[yn + We obtain the same relation between yn and yn+1 when (A) and (B) are applied. (b) 1 yn+1 = (1 + iλ∆t − λ2 ∆t2 )yn 2 Inductively, n 1 2 2 |yn | = 1 + iλ∆t − λ ∆t |y0 | 2 2 1 + iλ∆t − 1 λ2 ∆t2 = (1 − 1 λ2 ∆t2 )2 + (λ∆t)2 2 2 1 = 1 − λ2 ∆t2 + λ4 ∆t4 + λ2 ∆t2 4 1 4 4 = 1 + λ ∆t 4 > 1 for any fixed ∆t Hence |yn | → ∞ as n → ∞. NUS Math LaTeXify Proj Team Page: 3 of 5 NUS Mathematics Society MA3227 Numerical Analysis 2 AY 2009/2010 Sem 1 (c) Applying (C), we have, ∆t [iλyn + iλyn+1 ] 2 iλ∆t = (1 + )yn 2 yn+1 = yn + (1 − iλ∆t )yn+1 2 = 1 + iλ∆t 6= 0, |yn+1 | = |yn |. Since 1 − iλ∆t 2 2 Inductively, |yn | = |y0 |. Question 6 (a) Suppose λ is an arbitrary eigenvalue of A(k−1) , then λI − A(k−1) = 0 (k) (k) (k) λI − A = λI − R Q = λQ(k),T Q(k) − Q(k),T Q(k) R(k) Q(k) = Q(k),T λI − A(k−1) Q(k) =0 So λ is an eigenvalue of A(k) Similarly, we can prove any eigenvalue of A(k) is also eigenvalue of A(k−1) , so A(k) and A(k−1) have the same eigenvalues. (b) . . . • A = ~u~v T = [v1 ~u .. v2 ~u .. · · · .. vn ~u] Suppose q~1 = u ˆ 6=~0 and span{q~1 , q~2 , . . . q~n } = <n T q~i q~j = 0 if i 6= j and q~i T q~i = 1 Choose r~jk = vk k~uk if j = 1 and rjk = 0 if j 6= 1 Then A = QR. • QR = ~u~v T Comparing the first column r11 q~1 = v1 ~u So ~u = k q~1 for some k ∈ < and q~i T ~u = 0 for i = 2, 3, . . . , n RQ = QT QRQ = QT ~u~v T Q = (QT ~u)(QT ~v )T , where . . . QT ~v = [q~1 T ~v .. q~2 T ~v .. · · · .. q~n T ~v ]T . . . QT ~u = [q~1 T ~u .. ~0 .. · · · .. ~0]T So RQ is an upper triangular matrix with 0s below its first row, and the (1, 1) entry is (q~1 T ~u)(q~1 T ~v ) = k(q~1 T q~1 )(q~1 T ~v ) = ~uT ~v Eigenvalues of RQ are ~uT ~v and 0 Eigenvalues of A = QR are ~uT ~v and 0 Question 7 NUS Math LaTeXify Proj Team Page: 4 of 5 NUS Mathematics Society MA3227 Numerical Analysis 2 AY 2009/2010 Sem 1 (a) Z tn+1 f (t, y(t))dt y(tn+1 ) = y(tn ) + tn Using first order polynomial approximation, Z tn+1 t − tn−1 t − tn [ fn + fn−1 ]dt tn − tn−1 tn−1 − tn tn Z tn+1 fn fn−1 = [ (t − tn−1 ) − (t − tn )]dt ∆t ∆t tn ∆t = (3fn − fn−1 ) 2 A second order AB2 scheme is yn+1 = yn + ∆t 2 (3fn − fn−1 ) or yn+1 −yn ∆t = 23 fn − 12 fn−1 . (b) Consider applying AB − 2 to y 0 (t) = λy(t) yn+1 − yn = λ( 32 yn − 21 yn−1 ) ∆t Replacing yn by z n zn+1 − zn = λ( 32 zn − 21 zn−1 ) ∆t λ∆t = z−1 z n+1 − z n 1 n−1 = 1 1 3 n 2z − 2z 2 (3 − z ) At the boundary of the stability region, |z| = 1 or z = eiθ , θ ∈ <. In the programme, r = z − 1 and s = 12 (3 − z1 ). The programme varies θ and plots r/s, which is the boundary of the stability region. NUS Math LaTeXify Proj Team Page: 5 of 5 NUS Mathematics Society