hw8s - CME302 Homework 8 Computer Exercise A.M.B. November...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

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 ||rn||A-1 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 A-1 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 A-1 , but I wanted to. In practice only rn 2 can be computed well, though there are approximations for rn A-1 .) 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.

Ask a homework question - tutors are online