Chapter 6: Problem Solutions

Problems Home

 P6.1.1 P6.2.1 P6.3.1 P6.4.1 P6.1.2 P6.2.2 P6.3.2 P6.4.2 P6.1.3 P6.2.3 P6.3.3 P6.4.3 P6.1.4 P6.2.4 P6.3.4 P6.4.4 P6.1.5 P6.3.5 P6.1.6 P6.3.6 P6.1.7 P6.3.7 P6.3.8 P6.3.9 P6.3.10 P6.3.11

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