Systems of Linear Algebraic Equations - Numerical Methods in Engineering with MATLAB - page 89

the function Av(v) which returns the product Av . Thisfunctionmust besuppliedby
the user (see Example 2.18).Wemust also supply the starting vector x and the constant
(right-hand-side) vector b .

function [x,numIter] = conjGrad(func,x,b,epsilon)
%SolvesAx=bbyconjugategradientmethod.
% USAGE: [x,numIter] = conjGrad(func,x,b,epsilon)
% INPUT:
% func = handle of function that returns the vector A*v
% x = starting solution vector
% b = constant vector in A*x = b
% epsilon = error tolerance (default = 1.0e-9)
% OUTPUT:
% x = solution vector
% numIter = number of iterations carried out

if nargin == 3; epsilon = 1.0e-9; end
n = length(b);
r=b-feval(func,x);s=r;
for numIter = 1:n
u = feval(func,s);
alpha = dot(s,r)/dot(s,u);
x=x+alpha*s;
r=b-feval(func,x);
if sqrt(dot(r,r)) < epsilon
return

else

beta = -dot(r,u)/dot(s,u);
s=r+beta*s;

end

end
error('Too many iterations')

EXAMPLE 2.15
Solve the equations

⎡
⎣

⎤
⎦

⎡
⎣

⎤
⎦ =

⎡
⎣

⎤
⎦

4 − 11 x 1

x 2
x 3

12
−

− 1 4 − 2 1

5 1 − 2 4

by the Gauss-Seidel methodwithout relaxation.