Chapter 6 M-Files

M-File Home

Script Files

 ShowTri Uses LTriSol to get inverse of Forsythe Matrix. ShowTriD Illustrates tridiagonal system solving. ShowHessLU Illustrates Hessenberg LU factorization. ShowGE Illustrates Gaussian Elimination. ShowSparseSolve Examines \$\backslash\$ with sparse arrays. NoPivot Illustrates the dangers of no pivoting. ShowGEPiv Illustrates GEPiv. ShowGE2 Applies GE2 to three different examples. CondEgs Accuracy of some ill-conditioned Pascal linear systems.

Function Files

 LTriSol Solves lower triangular system Lx = b. UTriSol Solves upper triangular system Ux = b. LTriSolM Solves multiple right-hand-side lower triangular systems. TriDiLU LU factorization of a tridiagonal matrix. LBiDiSol Solves lower bidiagonal systems. UBiDiSol Solves upper bidiagonal systems. HessLU Hessenberg LU factorization. GE General LU factorization without pivoting. GEpiv General LU factorization with pivoting. GE2 Illustrates 2-by-2 GE in three-digit arithmetic.

ShowTri

```% Script File: ShowTri
% Inverse of the 5-by-5 Forsythe Matrix.

n = 5;
L = eye(n,n) - tril(ones(n,n),-1)
X = LTriSolM(L,eye(n,n))
```

```% Script File: ShowTriD
% Tridiagonal Solver test

clc
disp('Set L and U as follows:')
L = [1 0 0 0 0; 2 1 0 0 0 ; 0 3 1 0 0; 0 0 4 1 0; 0 0 0 5 1]
U = [2 3 0 0 0; 0 1 -1 0 0; 0 0 2 1 0; 0 0 0 3 1; 0 0 0 0 6]
%
input('Strike Any Key To Continue');
clc
disp('Form A = LU and set b so solution to Ax=b is ones(5,1).')
pause(2)
A = L*U
b = A*ones(5,1)
input('Strike Any Key To Continue');
%
clc
disp('Extract diagonal part of A:')
A = A
d = diag(A)'
input('Strike Any Key To Continue');
clc
disp('Extract subdiagonal part of A:')
A = A
e(2:5) = diag(A,-1)
input('Strike Any Key To Continue');
clc
disp('Extract superdiagonal part of A:')
A = A
f(1:4) = diag(A,1)
input('Strike Any Key To Continue');
%
clc
disp('Solve Ax = b via TriDiLU, LBiDisol, and UBiDiSol')
[l,u] = TriDiLU(d,e,f);
y = LBiDiSol(l,b);
x = UBiDiSol(u,f,y)
```

ShowHessLU

```% Script File: ShowHessLU
% Illustrates computation of a 5-by-5 LU factorization
% of upper Hessenberg system without pivoting.

clc
disp('Steps through the LU factorization of a 5-by-5')
disp('upper Hessenberg matrix.')
input('Strike Any Key to Continue');
clc
n = 5;
A = rand(n,n);
A = triu(A,-1);
[n,n] = size(A);
v = zeros(n,1);
for k=1:n-1
clc
Before = A
v(k+1) = A(k+1,k)/A(k,k);
disp(sprintf('Zero A(%1.0f,%1.0f)',k+1,k))
disp(sprintf('Multiplier = %7.4f / %7.4f = %7.4f',A(k+1,k),A(k,k),v(k+1)))
A(k+1,k:n) = A(k+1,k:n) - v(k+1)*A(k,k:n);
After = A
input('Strike Any Key to Continue')
end
U = A;
```

ShowGE

```% Script File: ShowGE
% Illustrates computation of a 5-by-5 LU factorization
% without pivoting.

clc
format short
disp('Steps through the LU factorization of a 5-by-5 matrix')
input('Strike Any Key to Continue.');
clc
n = 5;
A = magic(5);
ASave = A;
[n,n] = size(A);
v = zeros(n,1);
L = eye(n,n);
for k=1:n-1
v(k+1:n) = A(k+1:n,k)/A(k,k);
for i=k+1:n
clc
Before = A
disp(sprintf('Multiply row %1.0f by %7.4f / %7.4f ',k,A(i,k),A(k,k)))
disp(sprintf('and subtract from row %1.0f:',i))
A(i,k:n) = A(i,k:n) - v(i)*A(k,k:n);
After = A
input('Strike Any Key to Continue.');
end
clc
disp(sprintf('Have computed column %1.0f of L.',k))
disp(' ')
disp('Here is L so far:')
L(k+1:n,k) = v(k+1:n)
input('Strike Any Key to Continue.');
end
clc
disp('Factorization is complete:')
L = L
U = triu(A)
pause
clc
disp('Error check. Look at A - L*U:')
Error = ASave - L*U
```

```% Script File: ShowSparseSolve
% Illustrates how the \ operator can exploit sparsity

clc
disp('        n     Flops Full    Flops Sparse ')
disp('------------------------------------------')
for n=[25 50 100 200 400 800]
T = randn(n,n)+1000*eye(n,n);
T = triu(tril(T,1),-1); T(:,1) = .1; T(1,:) = .1;
b = randn(n,1);
flops(0); x = T\b; fullFlops   = flops;
T_sparse = sparse(T);
flops(0); x = T_sparse\b; sparseFlops = flops;
disp(sprintf('%10d  %10d  %10d ',n,fullFlops,sparseFlops))
end
```

NoPivot

```% Script File: NoPivot
% Examines solution to
%
%        [ delta 1 ; 1 1][x1;x2] = [1+delta;2]
%
% for a sequence of diminishing delta values.

clc
disp(' Delta            x(1)                   x(2)  ' )
disp('-----------------------------------------------------')
for delta = logspace(-2,-18,9)
A = [delta 1; 1 1];
b = [1+delta; 2];
L = [ 1 0; A(2,1)/A(1,1) 1];
U = [ A(1,1) A(1,2) ; 0 A(2,2)-L(2,1)*A(1,2)];
y(1) = b(1);
y(2) = b(2) - L(2,1)*y(1);
x(2) = y(2)/U(2,2);
x(1) = (y(1) - U(1,2)*x(2))/U(1,1);
disp(sprintf(' %5.0e   %20.15f  %20.15f',delta,x(1),x(2)))
end
```

ShowGEPiv

```% Script File: ShowGEpiv
% Illustrates computation of a 5-by-5 LU factorization
% with pivoting.

clc
format short
disp('Steps through Gaussian elimination of a 5-by-5')
disp('matrix showing pivoting.')
input('Strike Any Key to Continue.');
clc
n = 5;
A = magic(5);
[n,n] = size(A);
L = eye(n,n);
piv = 1:n;
for k=1:n-1
clc
[maxv,r] = max(abs(A(k:n,k)));
q = r+k-1;
Before = A
disp(sprintf('piv = [ %1.0f  %1.0f  %1.0f  %1.0f  %1.0f]',piv(1),piv(2),piv(3),piv(4),piv(5)))
disp(' ')
disp(sprintf('Interchange rows k = %1.0f and q = %1.0f',k,q))
piv([k q]) = piv([q k]);
A([k q],:) = A([q k],:);
After = A
disp(sprintf('piv = [ %1.0f  %1.0f  %1.0f  %1.0f  %1.0f]',piv(1),piv(2),piv(3),piv(4),piv(5)))
disp(' ')
disp(sprintf('Zero A(%1.0f:%1.0f,%1.0f):',k,k+1,k))
input('Strike Any Key to Continue.');
clc
if A(k,k) ~= 0
L(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - L(k+1:n,k)*A(k,k+1:n);
A(k+1:n,k) = zeros(n-k,1);
end
end
clc
home
L=L
U = A
piv = piv
```

ShowGE2

```% Script File: ShowGE2
% Illustrates 2-by-2 gaussian elimination in 3-digit floating point
% arithmetic.

A = [.981  .726; .529 .394];
ge2(A)
disp('Strike any key for next example.')
pause

A = [.981  .726; .529 .384];
ge2(A)
disp('Strike any key for next example.')
pause

A = [.981  .726; .529 .294];
ge2(A)
```

CondEgs

```% Script File: CondEgs
% Examines errors for a family of linear equation problems.

for n = [4 8 12 16]
clc
A = pascal(n);
disp(sprintf('cond(pascal(%2.0f)) = %8.4e',n,cond(A)));
disp('True solution is vector of ones. Computed solution =')
xTrue = ones(n,1);
b = A*xTrue;
format long
[L,U,piv] = GEpiv(A);
y = LTriSol(L,b(piv));
x = UTriSol(U,y)
format short
relerr = norm(x - xTrue)/norm(xTrue);
disp(sprintf('Relative error = %8.4e',relerr))
bound = eps*cond(A);
disp(sprintf('Predicted value = EPS*cond(A) = %8.4e',bound))
input('Strike Any Key to Continue.');
end
```

LTriSol

```  function x = LTriSol(L,b)
% x = LTriSol(L,b)
%
% Solves the nonsingular lower triangular system  Lx = b
% where L is n-by-n, b is n-by-1, and x is n-by-1.

n = length(b);
x = zeros(n,1);
for j=1:n-1
x(j) = b(j)/L(j,j);
b(j+1:n) = b(j+1:n) - L(j+1:n,j)*x(j);
end
x(n) = b(n)/L(n,n);
```

UTriSol

```  function x = UTriSol(U,b)
% x = UTriSol(U,b)
%
% Solves the nonsingular upper triangular system  Ux = b.
% where U is n-by-n, b is n-by-1, and X is n-by-1.

n = length(b);
x = zeros(n,1);
for j=n:-1:2
x(j) = b(j)/U(j,j);
b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j);
end
x(1) = b(1)/U(1,1);
```

LTriSolM

```  function X = LTriSolM(L,B)
% X = LTriSolM(L,B)
%
% Solves the nonsingular lower triangular system  LX = B
% where L is n-by-n, B is n-by-r, and X is n-by-r.

[n,r] = size(B);
X = zeros(n,r);
for j=1:n-1
X(j,1:r)     = B(j,1:r)/L(j,j);
B(j+1:n,1:r) = B(j+1:n,1:r) - L(j+1:n,j)*X(j,1:r);
end
X(n,1:r) = B(n,1:r)/L(n,n);
```

```% [l,u] = TriDiLU(d,e,f)
%
% Tridiagonal LU without pivoting. d,e,f are  n-vectors that define a tridiagonal matrix
% A = diag(e(2:n),-1) + diag(d) + diag(f(1:n-1),1). Assume that A has an LU factorization.
% l and u are n-vectors with the property that if L = eye + diag(l(2:n),-1)
% and U = diag(u) + diag(f(1:n-1),1), then A = LU.

n = length(d);
l = zeros(n,1);
u = zeros(n,1);
u(1) = d(1);
for i=2:n
l(i) = e(i)/u(i-1);
u(i) = d(i) - l(i)*f(i-1);
end
```

LBiDiSol

```  function x = LBiDiSol(l,b)
% x = LBiDiSol(l,b)
%
% Solves the n-by-n unit lower bidiagonal system  Lx = b
% where l and b are n-by-1 and L = I + diag(l(2:n),-1).

n = length(b);
x = zeros(n,1);
x(1) = b(1);
for i=2:n
x(i) = b(i) - l(i)*x(i-1);
end
```

UBiDiSol

```  function x = UBiDiSol(u,f,b)
% x = UBiDiSol(u,f,b)
%
% Solves the n-by-n nonsingular upper bidiagonal system  Ux = b
% where u, f, and b are n-by-1 and U = diag(u) + diag(f(1:n-1),1).

n = length(b);
x = zeros(n,1);
x(n) = b(n)/u(n);
for i=n-1:-1:1
x(i) = (b(i) - f(i)*x(i+1))/u(i);
end
```

HessLU

```  function [v,U] = HessLU(A)
% [v,U] = HessLU(A)
%
% Computes the factorization H = LU where H is an n-by-n upper Hessenberg
% and L is an n-by-n lower unit triangular and U is an n-by-n upper triangular
% matrix.
%
% v is a column n-by-1 vector with the property that L = I + diag(v(2:n),-1).

[n,n] = size(A);
v = zeros(n,1);
for k=1:n-1
v(k+1) = A(k+1,k)/A(k,k);
A(k+1,k:n) = A(k+1,k:n) - v(k+1)*A(k,k:n);
end
U = triu(A);
```

GE

```  function [L,U] = GE(A)
% [L,U] = GE(A)
%
% The LU factorization without pivoting. If A is n-by-n and has an
% LU factorization, then L is unit lower triangular and U is upper
% triangular so A = LU.

[n,n] = size(A);
for k=1:n-1
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
end
L = eye(n,n) + tril(A,-1);
U = triu(A);
```

GEpiv

```  function [L,U,piv] = GEpiv(A)
% [L,U,piv] = GE(A)
%
% The LU factorization with partial pivoting. If A is n-by-n, then
% piv is a permutation of the vector 1:n and L is unit lower triangular
% and U is upper triangular so A(piv,:) = LU. |L(i,j)|<=1 for all i and j.

[n,n] = size(A);
piv = 1:n;
for k=1:n-1
[maxv,r] = max(abs(A(k:n,k)));
q = r+k-1;
piv([k q]) = piv([q k]);
A([k q],:) = A([q k],:);
if A(k,k) ~= 0
A(k+1:n,k) = A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
end
end
L = eye(n,n) + tril(A,-1);
U = triu(A);
```

GE2

```  function GE2(A)
% GE2(A)
%
% Displays the results of 2-by-2 Gaussian elimination when it is applied to
% the linear system Ax = [A(1,1)-A(1,2);A(2,1)-A(2,2)] using 3-digit arithmetic.

clc
condA = cond(A);
a11 = represent(A(1,1));
a12 = represent(A(1,2));
a21 = represent(A(2,1));
a22 = represent(A(2,2));
b1 = float(a11,a12,'-');
b2 = float(a21,a22,'-');
disp(['Stored A =   ' pretty(a11)  '   ' pretty(a12) ])
disp(['             ' pretty(a21)  '   ' pretty(a22) ])
L11 = represent(1);
L12 = represent(0);
L21 = float(a21,a11,'/');
L22 = represent(1);
disp(' ');
disp(['Computed L = ' pretty(L11)  '   ' pretty(L12)])
disp(['             ' pretty(L21)  '   ' pretty(L22)])
U11 = a11;
U12 = a12;
U21 = represent(0);
U22 = float(a22,float(L21,a12,'*'),'-');
disp(' ');
disp(['Computed U = ' pretty(U11)  '   ' pretty(U12)])
disp(['             ' pretty(U21)  '   ' pretty(U22)])
y1 = b1;
y2 = float(b2,float(L21,y1,'*'),'-');
x2 = float(y2,U22,'/');
x1 = float(float(y1,float(U12,x2,'*'),'-'),U11,'/');
xe1 = represent(1);
xe2 = represent(-1);
disp(' ');
disp(['Exact b           =  '  pretty(b1) ])
disp(['                     '   pretty(b2)  ])
disp(' ');
disp(['Exact Solution    =  '  pretty(xe1) ])
disp(['                    '   pretty(xe2)  ])
disp(' ');
disp(['Computed Solution =  ' pretty(x1)])
disp(['                    '  pretty(x2)])
disp(sprintf('\ncond(A) = %5.3e',condA))
xtilde = [convert(x1);convert(x2)];
x = [1;-1];
b = A*x;
r = A*xtilde-b;
E = -r*xtilde'/(xtilde'*xtilde);
disp(' ');
disp('Computed solution solves (A+E)*x = b where')
disp(sprintf('\n                   E =  %12.6f %12.6f',E(1,1),E(1,2)));
disp(sprintf('                        %12.6f %12.6f',E(2,1),E(2,2)));
disp(' ')
```