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

ShowTriD 

% 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

ShowSparseSolve 

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

 

TriDiLU 

% [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(' ')