Function [x , normrn] = gmres(A,b,maxiter)



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);

[pic]

[pic]

% 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);

[pic] [pic]

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download