Docstoc

GMRES Algorithm

Document Sample
GMRES Algorithm Powered By Docstoc
					function [x , normrn] = gmres(A,b,maxiter) % solves A x = b using gmres % input: A - m by m matrix % b - m by 1 vector % output: x - approximate solution % normrn - norm(b-A*x) at each step of algorithm % Also % plot normrn/ norm(b) versus n (step number) % Example: % m = 200; % A = eigmat(2,1/2,m); % xtrue = randn(m,1); % b = A*xtrue; % [x , normrn] = gmres(A,b,20); % CAUTION: Matlab has a built-in function with the same % name (gmres). If you want to use the new version % of gmres make sure your active folder contains the % new version of gmres in file gmres.m. Q = []; H = 0; m = length(A); normb = norm(b); normrn=normb; Q(:,1) = b / normb; for n = 1:maxiter % Arnoldi step calculations (Algorithm 33.1 for Trefethen and Bau) v = A*Q(:,n); for j = 1:n H(j,n) = Q(:,j)'* v; v = v - H(j,n)*Q(:,j); end Hn = H(1:n,1:n); H(n+1,n) = norm(v); if H(n+1,n) == 0, break, end % breakdown so stop Q(:,n+1) = v / H(n+1,n); e1 = [1;zeros(n,1)]; y = H \ (normb*e1); % This can be done much more quickly % using Givens rotations. % For simplicity we just use Matlab's \ normrn = [normrn,norm(H*y-normb*e1)]; % remember residual norm end x = Q(:,1:n)*y; clf semilogy(0:n,normrn/normb,'--o'),shg xlabel('step number of gmres algorithm') ylabel(' norm( b - A * xn ) / norm(b)') title('convergence of residual in gmres') grid

% file example35p1.m % example 35.1 from Trefethen + Bau m=200; maxiter = 20; A=eigmat(2,.5,m); lam=eig(A); figure(1) clf plot(lam,'x'); axis square title('eigenvalues of A') ylabel('imag part') xlabel('real part') grid xtrue=randn(m,1); b = A*xtrue; figure(2) [x , normrn] = gmres(A,b,maxiter);
eigenvalues of A 0.5 0.4 0.3 0.2 0.1

imag part

0 -0.1 -0.2 -0.3 -0.4 -0.5 1.4

1.6

1.8

2 2.2 real part

2.4

2.6

2.8

10

0

convergence of residual in gmres

10

-2

norm( b - A * xn ) / norm(b)

10

-4

10

-6

10

-8

10

-10

10

-12

10

-14

0

2

4

6 8 10 12 14 step number of gmres algorithm

16

18

20

% file example35p2.m % example 35.2 from Trefethen + Bau m=200; maxiter = 20; th = (0:m-1)*pi / (m-1); d=(-2 +2* sin(th))+sqrt(-1)*cos(th); A=eigmat(2,.5,m)+diag(d); lam=eig(A); figure(1) clf plot(lam,'x'); axis square title('eigenvalues of A') ylabel('imag part') xlabel('real part') grid xtrue=randn(m,1); b = A*xtrue; figure(2) [x , normrn] = gmres(A,b,maxiter);

eigenvalues of A 1.5

1

0.5

imag part

0

-0.5

-1

-1.5 -0.5

0

0.5

1 real part

1.5

2

2.5

10

0

convergence of residual in gmres

norm( b - A * xn ) / norm(b)

10

-1

10

-2

10

-3

0

2

4

6 8 10 12 14 step number of gmres algorithm

16

18

20


				
DOCUMENT INFO
Shared By:
Stats:
views:3890
posted:8/26/2009
language:English
pages:3