Page 244 - MATLAB an introduction with applications

Numerical Methods ——— 229

n = size(A,1);
for k = 1:n–2
u = A(k+1:n,k);
uMag = sqrt(dot(u,u));
if u(1)<0;uMag = –uMag; end
u(1) = u(1)+uMag;
A(k+1:n,k) = u;
H = dot(u,u)/2;
v = A(k+1:n,k+1:n)*u/H;
g = dot(u,v)/(2*H);
v = v–g*u;
A(k+1:n,k+1:n) = A(k+1:n,k+1:n)–v*u’–u*v’;
A(k,k+1) = –uMag;
d = diag(A);
c = diag(A,1);
end

Example E4.12: Solve the system of equations given below by Householder’s factorization method:
x
4  2 1 1   10

1

 210 2 1   20 


x
   =   
2
1  2 4 2   30 
x
3
 
x

1   2 4 8    40 
  
4
Solution:
>> A = [4 2 1 1;2 10 2 1;1 2 4 2;1 2 4 8];
% Householder reduction of A to tridiagonal form A = [c\d\c]
% Extract c and by d = diag(A) andc = diag(A,1)
>> householder(A)
ans =
4.0000 –2.4495 1.0000 1.0000
4.4495 12.4082 –1.8257 –2.6330
1.0000 2.6422 2.4400 0.3427
1.0000 –1.6330 2.3427 7.1517
>> b = [10;20;30;40];
>> A = [4 2 1 1;2 10 2 1;1 2 4 2;1 2 4 8];
>> [L,U] = lu(A)
L =
1.0000 0 0 0
0.5000 1.0000 0 0
0.2500 0.1667 1.0000 0
0.2500 0.1667 1.0000 1.0000