Chapter 6: Problem Solutions

Problems Home

P6.1.1

% Problem P6_1_1
%
% Solving certain singular lower triangular systems.

clc 
n = 6;
L = tril(randn(n,n));
L(1,1) = 0
x = LTriSol1(L)
disp('L*x =')
disp(L*x)

  function x = LTriSol1(L)
% x = LTriSol1(L)
%
% L is an n-by-n  lower triangular matrix. L(2:n,2:n) is
% nonsingular and L(1,1) = 0.
%
% x satisfies Lx = 0 with x(1) = -1. 

[n,n] = size(L);
x = zeros(n,1);
x(1) = -1;
x(2:n) = LTriSol(L(2:n,2:n),L(2:n,1));

P6.1.2

% Problem P6_1_2
%
% Solving multiple righthand side upper triangular systems

clc
n = 6;
r = 4;
U = triu(randn(n,n))
B = randn(n,r)
X = UTriSolM(U,B)
disp('U*X - B')
disp(U*X-B)


  function X = UTriSolM(U,B)
% X = UTriSolM(U,B)
%
% U is an n-by-n nonsingular upper triangular matrix
% B is am n-by-r matrix.
%
% X satisfies  UX = B.

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

P6.1.3

% Problem P6_1_3.
%
% Solving the upper triangular Sylvester equation SX-XT=B.

clc 
m = 6;
n = 4;
S = triu(randn(m,m))
T = triu(randn(n,n))
B = randn(m,n)
X = Sylvester(S,T,B)
disp(sprintf('norm(SX-XT-B)/norm(X) = %8.3e',norm(S*X-X*T-B)/norm(X)))


  function X = Sylvester(S,T,B)
% X = Sylvester(S,T,B)
%
% S is an m-by-m upper triangular 
% T is an n-by-n upper triangular with no diagonal entries in
%         common with S.
% B is an m-by-n matrix.
% 
% X is an m-by-n so SX-XT=B.

[m,n] = size(B);
X = zeros(m,n);
for k=1:n
   if k==1
      v = B(:,k);
   else
      v = B(:,k) + X(:,1:k-1)*T(1:k-1,k);
   end
   X(:,k) = UTriSol(S-T(k,k)*eye(m,m),v);
end

P6.1.4

% Problem P6_1_4
%
% Solving the lower triangular Sylvester equation SX-XT=B.

clc 
m = 6;
n = 4;
S = tril(randn(m,m))
T = tril(randn(n,n))
B = randn(m,n)
X = Sylvester1(S,T,B)
disp(sprintf('norm(SX-XT-B)/norm(X) = %8.3e',norm(S*X-X*T-B)/norm(X)))


 function X = Sylvester1(S,T,B)
% X = Sylvester1(S,T,B)
%
% S is an m-by-m lower triangular. 
% T is an n-by-n lower triangular with no diagonal entries in
%         common with S.
% B is an m-by-n matrix.
% 
% X is an m-by-n so SX-XT=B.
%

[m,n] = size(B);
X = zeros(m,n);
for k=n:-1:1
	if k==n
	   v = B(:,n);
	else
	   v = B(:,k) + X(:,k+1:n)*T(k+1:n,k);
	end
   X(:,k) = LTriSol(S-T(k,k)*eye(m,m),v);
end

P6.1.5

% Problem P6_1_5
%
% The inverse of upper triangular matrix.

clc 
n = 6;
U = triu(randn(n,n))
X = UTriInv(U)
disp(sprintf('norm(U*X - eye(n,n))/norm(X) = %8.3e',norm(U*X - eye(n,n))/norm(X) ));

  function X = UTriInv(U)
% X = UTriInv(U)
%
% U is an n-by-n nonsingular upper triangular matrix.
%
% X is the inverse of U.

[n,n] = size(U);
X = zeros(n,n);
for k=1:n
   v = zeros(k,1);
   v(k) = 1;
   X(1:k,k) = UTriSol(U(1:k,1:k),v);
end

P6.1.6

% Problem P6_1_6
%
% Inverse of Forsythe matrix.

n = input('Enter n: ');
A = tril(-ones(n,n) + 2*eye(n,n))
A_inv = eye(n,n);
for j = 1:n
   for i=j+1:n
      A_inv(i,j) = 2^(i-j-1);	
   end
end
A_inv = A_inv

P6.1.7

Not available.

P6.2.1

% Problem P6_2_1
%
% Solving lower Hessenberg systems

n = 6;
A = triu(randn(n,n),-1)
b = A'*ones(n,1)
x = HessTrans(A,b)



  function x = HessTrans(A,b)
% x = HessTrans(A,b)
%
% A is an n-by-n upper Hessenberg matrix.
% b iS an n-by-1 vector.
% 
% x   solves A'*x = b. 

n = length(b); 
[v,U] = HessLU(A);
% (U'L')x = b
y = LTriSol(U',b);
x = UBidiSol(ones(n,1),v(2:n),y);

P6.2.2

% Problem P6_2_2
%
% Compare CubicSpline and CubicSpline1 where the latter exploits the
% tridiagonal nature of the linear system.

clc 
disp('Complete Spline')
disp(' ')
disp('           CubicSpline          CubicSpline1')
disp('  n      flops      time      flops       time')
disp('-----------------------------------------------')
for n = [20 40 80 160 320]
   x = linspace(0,pi/2,n)';
   y = sin(x);
   
   flops(0)
   tic;
   [a,b,c,d] = CubicSpline(x,y,1,1,0);
   t1 = toc;
   f1 = flops;
   
   flops(0)
   tic;
   [a1,b1,c1,d1] = CubicSpline1(x,y,1,1,0);
   t2 = toc;
   f2 = flops;
   disp(sprintf('%4.0f    %6.0f  %8.3f     %6.0f   %8.3f',n,f1,t1,f2,t2))
end
disp(' ')
disp(' ')	
disp('Second Derivative Spline')
disp(' ')
disp('           CubicSpline          CubicSpline1')
disp('  n      flops      time      flops       time')
disp('-----------------------------------------------')
for n = [20 40 80 160 320]
   x = linspace(0,pi/2,n)';
   y = sin(x);	
   flops(0)
   tic;    
   [a,b,c,d] = CubicSpline(x,y,2,0,1);
   t1 = toc;
   f1 = flops;
   flops(0)
   tic;    
   [a1,b1,c1,d1] = CubicSpline1(x,y,2,0,1);
   t2 = toc;
   f2 = flops;
   
   disp(sprintf('%4.0f    %6.0f  %8.3f     %6.0f   %8.3f',n,f1,t1,f2,t2))
end

disp(' ')
disp(' ')	
disp('Not-a-Knot Spline')
disp(' ')
disp('           CubicSpline          CubicSpline1')
disp('  n      flops      time      flops       time')
disp('-----------------------------------------------')
for n = [20 40 80 160 320]
   x = linspace(0,pi/2,n)';
   y = sin(x);		
   flops(0)
   tic;   
   [a,b,c,d] = CubicSpline(x,y); 
   t1 = toc;
   f1 = flops;  
   flops(0)
   tic;   
   [a1,b1,c1,d1] = CubicSpline1(x,y); 
   t2 = toc;
   f2 = flops; 
   disp(sprintf('%4.0f    %6.0f  %8.3f     %6.0f   %8.3f',n,f1,t1,f2,t2))
end



  function [a,b,c,d] = CubicSpline(x,y,derivative,muL,muR)
% [a,b,c,d] = CubicSpline(x,y,derivative,muL,muR)
% Cubic spline interpolation with prescribed end conditions.
% 
% x,y are column n-vectors. It is assumed that n >= 4 and x(1) < ... x(n).
% derivative is an integer (1 or 2) that specifies the order of the endpoint derivatives.
% muL and muR are the endpoint values of this derivative.
%
% a,b,c, and d are column (n-1)-vectors that define the spline S(z). On [x(i),x(i+1)], 
%  
%          S(z) =  a(i) + b(i)(z-x(i)) + c(i)(z-x(i))^2 + d(i)(z-x(i))^2(z-x(i+1).
%
% Usage:
%   [a,b,c,d] = CubicSpline(x,y,1,muL,muR)   S'(x(1))  = muL, S'(x(n))  = muR
%   [a,b,c,d] = CubicSpline(x,y,2,muL,muR)   S''(x(1)) = muL, S''(x(n)) = muR
%   [a,b,c,d] = CubicSpline(x,y)             S'''(z) continuous at x(2) and x(n-1)
       

% First, set up all but the first and last equations that
% define the vector of interior knot slopes s(2:n-1).

n = length(x); 
Dx = diff(x);
yp = diff(y) ./ Dx;
T = zeros(n-2,n-2);
r = zeros(n-2,1);
for i=2:n-3
   T(i,i)   = 2*(Dx(i) + Dx(i+1));
   T(i,i-1) = Dx(i+1);
   T(i,i+1) = Dx(i);
   r(i)     = 3*(Dx(i+1)*yp(i) + Dx(i)*yp(i+1));
end

% For each of the 3 cases, finish setting up the linear system,
% solve the system, and set s(1:n) to be the vector of slopes.

if nargin==5
   %Derivative information available.
   if derivative==1
      % End values for S'(z) specified.         
      T(1,1) = 2*(Dx(1) + Dx(2));
      T(1,2) = Dx(1);
      r(1) = 3*(Dx(2)*yp(1)+Dx(1)*yp(2)) - Dx(2)*muL;
      T(n-2,n-2) = 2*(Dx(n-2)+Dx(n-1));
      T(n-2,n-3) = Dx(n-1);
      r(n-2) = 3*(Dx(n-1)*yp(n-2) + Dx(n-2)*yp(n-1)) -Dx(n-2)*muR;
      s = [muL; T\r; muR];
   end
   
   if derivative==2
      % End values for S''(z) specified.  
      T(1,1) = 2*Dx(1) + 1.5*Dx(2);
      T(1,2) = Dx(1);
      r(1) = 1.5*Dx(2)*yp(1) + 3*Dx(1)*yp(2) + Dx(1)*Dx(2)*muL/4;
      T(n-2,n-2) = 1.5*Dx(n-2)+2*Dx(n-1);
      T(n-2,n-3) = Dx(n-1);
      r(n-2) = 3*Dx(n-1)*yp(n-2) + 1.5*Dx(n-2)*yp(n-1)-Dx(n-1)*Dx(n-2)*muR/4;
      stilde = T\r;
      s1 = (3*yp(1) - stilde(1) - muL*Dx(1)/2)/2;
      sn = (3*yp(n-1) - stilde(n-2) + muR*Dx(n-1)/2)/2;
      s = [s1;stilde;sn];
   end;
else
   % No derivative information. Compute the not-a-knot spline.
   q = Dx(1)*Dx(1)/Dx(2);
   T(1,1) = 2*Dx(1) +Dx(2) + q;
   T(1,2) = Dx(1) + q;
   r(1) = Dx(2)*yp(1) + Dx(1)*yp(2)+2*yp(2)*(q+Dx(1));
   q = Dx(n-1)*Dx(n-1)/Dx(n-2);
   T(n-2,n-2) = 2*Dx(n-1) + Dx(n-2)+q;
   T(n-2,n-3) = Dx(n-1)+q;
   r(n-2) = Dx(n-1)*yp(n-2) + Dx(n-2)*yp(n-1) +2*yp(n-2)*(Dx(n-1)+q);
   
   [l,u] = TriDiLU(d,e,f);		 
   stilde = UBiDiSol(u,f,LBiDiSol(l,r));
  
   s1 = -stilde(1)+2*yp(1);
   s1 = s1 + ((Dx(1)/Dx(2))^2)*(stilde(1)+stilde(2)-2*yp(2));
   sn = -stilde(n-2) +2*yp(n-1);
   sn = sn+((Dx(n-1)/Dx(n-2))^2)*(stilde(n-3)+stilde(n-2)-2*yp(n-2));
   s = [s1;stilde;sn];
end

% Compute the a,b,c,d vectors.
   
a = y(1:n-1);
b = s(1:n-1);
c = (yp - s(1:n-1)) ./ Dx;
d = (s(2:n) + s(1:n-1) - 2*yp) ./ (Dx.* Dx);

P6.2.3

% Problem P6_2_3
%
% Solve HX-XT=B where H is upper Hessenberg and T is upper triangular.

clc  
m = 6;
n = 4;
H = triu(randn(m,m),-1)
T = triu(randn(n,n))
B = randn(m,n)
X = SylvesterH(H,T,B)
disp(sprintf('norm(H*X-X*T-B)/norm(X) = %8.3e ',norm(H*X-X*T-B)/norm(X)))


  function X = SylvesterH(H,T,B)
% X = SylvesterH(H,T,B)
%
% S is an m-by-m upper Hessenberg matrix.
% T is an n-by-n upper triangular matrix with the property that HX-XT=B
%     is a nonsingular linear system.
% B is an m-by-n matrix.
% 
% X is an m-by-n matrix that satisfies HX-XT=B.
%
    [m,n] = size(B);
	X = zeros(m,n);
	for k=1:n
	   if k==1
	      r = B(:,k);
	   else
	      r = B(:,k) + X(:,1:k-1)*T(1:k-1,k);
	   end
	   [v,U] = HessLU(H-T(k,k)*eye(m,m));
	   y = LBiDiSol(v,r);
	   X(:,k) = UTriSol(U,y);
	end

P6.2.4

 % Clear command window to start.
clc;

% Create W.
d=[1 2 3 4 5 6];
e=[3 5 0 5 0 5];
f=[2 7 2 7 2 8];

% Create b.
b=[7 7 7 7 7 7]';

% Create alpha, beta.
alpha = 2;
beta = 3;

% Solve Tx=b, where T = W+alpha*e_n*e_1'+beta*e_1*e_n';
x = TriCorner(d,e,f,alpha,beta,b);



function x = TriCorner(d,e,f,alpha,beta,b)

% W is a tridiagonal matrix with d encoding the main diagonal, 
%     f encoding the super diagonal, and e encoding the subdiagonal.
%
% T = tridiagonal + corners.
% T = W+alpha*e_n*e_1'+beta*e_1*e_n', 
%     where e_1, e_n = first, last cols. of identity matrix.
%
% Goal:  Find x such that Tx = b.  Do this efficiently.
% Note:  x = z-Y*[1+Y(1,1) Y(1,2); Y(n,1) 1+Y(n,2)]^(-1)*[z1; zn]
%     where z = W^(-1)*b, Y = W^(-1)*[alpha*e_n beta*e_1].
%
%
% Solution:

% Initialization:
n = length(b);
e_1 = zeros(n,1); e_1(1) = 1;
e_n = zeros(n,1); e_n(n) = 1;

% Compute LU factorization of W:
[l,u] = TriDiLU(d,e,f);    % See p. 217 for TriDiLU.m.

% Solve Wz = b by noting LUz = b. 
% First step is to solve Ly=b.  Second step is to solve Uz=y.
% Note L = lower bidiagonal, U = upper bidiagonal, since W = tridiagonal.
y = LBiDiSol(l,b);   % See p.217 for LBiDiSol.m.
z = UBiDiSol(u,f,y);   % See p.218 for UBiDiSol.m.

% Solve for y1, y2 (columns of Y).  First we're solving Wy1 = alpha*e_n.
% Second we're solving Wy2 = beta*e_1.

% Use the same approach as above.  Repeat linear system solves twice.
t = LBiDiSol(l,alpha*e_n);
y1 = UBiDiSol(u,f,t);

t = LBiDiSol(l,beta*e_1);
y2 = UBiDiSol(u,f,t);

% Assemble Y.
Y = [y1 y2];

% Let V = [1+Y(1,1) Y(1,2); Y(n,1) 1+Y(n,2)]^(-1).
% Compute V efficiently.
% Compute determinant.
delta = 1/((1+Y(1,1))*(1+Y(n,2))-(Y(1,2))*(Y(n,1)));
V(1,1) = 1+Y(n,2);
V(1,2) = -Y(1,2);
V(2,1) = -Y(n,1);
V(2,2) = 1+Y(1,1);
V = delta*V;

% Let s = [z(1); z(n)].
s=[z(1); z(n)];

% Final computations.
x=z-Y*V*s;

% Check!
W = diag(e(2:n),-1) + diag(d) + diag(f(1:n-1),1);
T = W+alpha*e_n*e_1'+beta*e_1*e_n';
difference = norm(T*x-b);
disp(['The norm of the difference is ' num2str(difference) '.']);

P6.3.1

% Problem P6_3_1
%
% Multiple righthand side solver.

clc
n = 6;
r = 4;
A = randn(n,n)
B = randn(n,r)
X = MultRhs(A,B)
disp(sprintf('norm(AX-B)/(norm(A)*norm(X)) = %8.3e',norm(A*X-B)/(norm(A)*norm(X))))


  function X = MultRhs(A,B)
% X = MultRhs(A,B)
%
% A is an n-by-n nonsingular matrix.
% B is an n-by-r matrix.
%
% X is an n-by-r matrix that satisfies AX = B.
[n,r] = size(B);
X = zeros(n,r);
[L,U,piv] = GEpiv(A);
for k=1:r
   y = LTriSol(L,B(piv,k));
   X(:,k) = UTriSol(U,y);
end

P6.3.2

% Problem P6_3_2
%
% Illustrates computation of a 5-by-5 LU factorization
% with parameterized pivoting.

clc
format short
disp('Steps through Gaussian elimination of a 5-by-5')
disp('matrix showing pivoting.')
alfa = input('Enter alfa (0 < alfa <= 1): ');
input('Strike Any Key to Continue.');
clc
n = 7;
A = magic(n);
[n,n] = size(A); 
L = eye(n,n);
piv = 1:n;
for k=1:n-1
   clc
   [maxv,r] = max(abs(A(k+1:n,k)));
   if abs(A(k,k)) >= alfa*maxv
      % No Interchange
      q = k;
   else
      q = r+k;
   end
   Before = A
   disp(sprintf('piv = [ %1.0f  %1.0f  %1.0f  %1.0f  %1.0f  %1.0f  %1.0f  %1.0f ]',piv))
   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  %1.0f  %1.0f  %1.0f ]',piv))
   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
L=L
U = A
piv = piv

P6.3.3

% Problem P6_3_3
%
% Computing a specified entry of A^-1.

clc
n = 6;
A = hilb(n)
X = zeros(n,n);
for i=1:n
   for j=1:n
      X(i,j) = Ainv1(A,i,j);
   end
end
disp('Computed inverse using Ainv:')
disp(' ' )
for i=1:n
   disp(sprintf('  %10.0f  %10.0f  %10.0f  %10.0f  %10.0f  %10.0f',X(i,:)))
end
Xexact = invhilb(n)


  function z = Ainv(A,i,j)
% z = Ainv(A,i,j)
%
% A is an n-by-n nonsingular matrix.
% i,j  integers with 1<=i<=n and 1<=j<=n.
%
% z is the (i,j) entry of A^-1.

[n,n] = size(A);
[L,U,piv] = GEpiv(A);
b = zeros(n,1);
b(j) = 1;
y = LTriSol(L,b(piv));
% Just need the i-th component of the solution to Ux = y.
xtilde = UTriSol(U(i:n,i:n),y(i:n));
z = xtilde(1);

P6.3.4

% Problem  P6_3_4
%
% Block back substitution, the 2-by-2 case.

clc 
n = 6;
A = randn(n,n);
B = randn(n,n);
C = randn(n,n);
g = randn(n,1);
h = randn(n,1);
[y,z] = BlockSolve(A,B,C,g,h);
xtilde = [y;z];
x = [A B; zeros(n,n) C]\[g;h];
disp('           BlockSolve            Standard Solve')
disp('---------------------------------------------------')
for k=1:2*n
   disp(sprintf(' %23.15f  %23.15f',xtilde(k),x(k)))
end


  function [y,z] = BlockSolve(A,B,C,g,h)
% [y,z] = BlockSolve(A,B,C,g,h)
%
% A,B,C are n-by-n matrices with A and C nonsingular
% g,h are column n-vectors
%
% y,z are column n-vectors so Ay+Bz = g and Cz = h.

z = C\h;
y = A\(g-B*z);

P6.3.5

% Problem P6_3_5
%
% Examines the solution to ABCx = f.

clc 
disp('f1 = ProdSolve flops')
disp('f2 = Multiply and Solve flops')
disp(' ')
disp(' ')
disp('   n      f2/f1  ')
disp('-----------------')
for n = [10 20 40 80 160]
   A = randn(n,n);
   B = randn(n,n);
   C = randn(n,n);
   f = randn(n,1);
   flops(0);
   x = ProdSolve(A,B,C,f);
   f1 = flops;
   flops(0)
   M = A*B*C;
   [LM,UM,pivM] = GEpiv(M);
   y = LTriSol(LM,f(pivM));
   xtilde = UTriSol(UM,y);
   f2 = flops;
   disp(sprintf('%4.0f   %8.3f',n,f2/f1))
end


  function x = ProdSolve(A,B,C,f)
% x = ProdSolve(A,B,C,f)
%
% A,B,C are n-by-n nonsingular matrices
% f is a column n-vector
%
% x is a columnn-vector so ABCx = f
%
[LA,UA,pivA] = GEpiv(A);
yA = Ltrisol(LA,f(pivA));
xA = UtriSol(UA,yA);
[LB,UB,pivB] = GEpiv(B);
yB = Ltrisol(LA,xA(pivB));
xB = UtriSol(UA,yB);
[LC,UC,pivC] = GEpiv(C);
yC = Ltrisol(LA,xB(pivC));
x  = UtriSol(UC,yC);

P6.3.6

Not Available

P6.3.7

Not Available

P6.3.8

% Problem P6_3_8
%
% Intersection in 3-space.

P = [1;0;0];
Q = [0;1;1];
R = [1;1;1];
S = [1;0;1];

t = Intersect1(P,Q,R,S)




  function t = Intersect1(P,Q,R,S)
% t = Intersect(P,Q,R,S)
%
% P, Q, R, and S are column 3-vectors. 
% t is a scalar so that P + t*Q is a linear combination
% or R and S.

% This means P + t*Q = x(1)*R + x(2)*S for some pair of scalars
% x(1) and x(2). In other words,
%
%          P = [R S -Q]*[x(1);x(2);t]

A = [R S -Q]
x = A\P;
t = x(3);

P6.3.9

Not Available

P6.3.10

Not Available

P6.3.11

Not Available

P6.4.1

% Problem P6_4_1
%
% Explore relative error in computed x.

clc 

del = .0000002;
A = [ 981 726 ; 529 del+726*529/981]
disp(sprintf('cond(A) = %10.3e',cond(A)))
disp(' ')
xexact = [ 1.234567890123456; .0000123456789012];
b = A*xexact;
x = A\b;
disp(sprintf('xexact = %20.16f  %20.16f',xexact))
disp(sprintf('x      = %20.16f  %20.16f',x))

P6.4.2

% Problem P6_4_2
%
% Explore computed relative error.

clc 
xexact = [.1234567890123456 ; 1234.567890123456];
for OurEps = logspace(-5,-16,12)
   % Simulate a computed x so norm(x - xexact)/norm(xexact) is about OurEps*cond
   x = xexact + rand(2,1)*10^4*OurEps*norm(xexact);
   disp(' ')
   disp(sprintf('OurEps = %8.3e',OurEps))
   disp(' ')
   disp(sprintf('xexact = %20.16f  %20.16f',xexact))
   disp(sprintf('x      = %20.16f  %20.16f',x))
end

P6.4.3

Not Available

P6.4.4

Not Available