Page 245 - MATLAB an introduction with applications

230 ——— MATLAB: An Introduction with Applications

U =
4.0000 2.0000 1.0000 1.0000
0 9.0000 1.5000 0.5000
0 0 3.5000 1.6667
0 0 0 6.0000
>> L*U
ans =
4 2 1 1
2 10 2 1
1 2 4 2
1 2 4 8
>> d = L\b
d =
10
15
25
10
>> x = U\d
x =
0.2381
0.5159
6.3492
1.6667
function A = householder(A)
% Householder reduction of A to tridiagonal form A = [c\d\c]
% Extract c and d by d = diag(A), c = diag(A,1)
% Usage: A = householder(A)
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);