This preview shows page 1. Sign up to view the full content.
Unformatted text preview: CME302 Homework 8 Computer Exercise A.M.B. November 16, 2009
10
0 (T&B 38.6) The code for this problem follows:
function [x rc ra] = cg(A,b,nits) n = length(b); x = zeros(n,1); r = b; p = r; rc = zeros(1,nits+1); ra = rc; rc(1) = norm(b); ra = rc(1); for(i = 1:nits) Ap = A*p; r1 = dot(r,r); alpha = r1/dot(p,Ap); x = x + alpha*p; r = r  alpha*Ap; r2 = dot(r,r); rc(i+1) = sqrt(r2); ra(i+1) = norm(b  A*x); beta = r2/r1; p = r + beta*p; end function [x r] = steep(A,b,nits) n = length(b); x = zeros(n,1); p = b; r = zeros(1,nits+1); r(1) = norm(b); for(i = 1:nits) alpha = dot(p,p)/dot(p,A*p); x = x + alpha*p; p = b  A*x; r(i+1) = norm(p); end 10 5 10 10 5 10 15 20 10 15 Computed CG rn Actual CG r 
n Actual SD rn Actual CG rnA1 Theorem 38.5 bound 10 20 30 40 50 60 Iteration n 70 80 90 100 Figure 1: Problem 38.6. When r is small, so is p . But x is not. At some point in the final iterations, when x is updated, numerically it remains the same and so the actual residual is unchanged; in contrast, the computed r numerically changes. In steepest descent, as in CG, rn 2 does not decrease monotonically; rather 1 than quantity 2 xT Ax  bT x does. Theorem 38.5 bounds en A , not rn 2 . Nonetheless, if the bound is multipled by b , it gives an upper bound for rn 2 in this problem. Steepest descent converges much more slowly than CG more slowly than even the upper bound given by Thm. 38.5. Figure 1 shows the convergence results for the problem Ax = b. The subplot zooms in on the first several iterations. Comments: In CG, en A = rn A1 decreases monotonically, but rn 2 does not. This is evident in our problem. The green and blue curves increase during approximately iterations 12 to 18, while the cyan curve is monotonically decreasing. (You were not asked to plot rn A1 , but I wanted to. In practice only rn 2 can be computed well, though there are approximations for rn A1 .) The actual residual norm levels off, while the computed one decreases. Consider the code
x = x + alpha*p; r = r  alpha*Ap; 1 ...
View
Full
Document
This document was uploaded on 06/17/2010.
 Fall '09

Click to edit the document details