APMA1170 - Homework 4

Homework 4
Problem 2.31
Let A ∈ Rn×n and let x, y, and z be n-vectors such that Ax = b and Ay = b + z.then prove that
||z||2 −1 −1 exists).
||A||2 ≤ ||x − y||2 ≤ ||A ||2 ||x||2 (assuming that A

Solution:
From the two linear systems, we get A(y − x) = z, and hence ||z||2 ≤ ||A||2 ||x − y||2 . Equaly, we
||z||2
have y − x = A−1 z, so ||x − y||2 ≤ ||A−1 ||2 ||x||2 . Then we get the conclution ||A|| 2
≤ ||x − y||2 ≤
−1
||A ||2 ||x||2 .

Problem 4.17
(a) Let A be an orthogonal matrix. then show that Cond2 (A) = 1.
(b) Show that Cond2 (A) = 1 if and only if a is a scalar multiple of an orthogonal matrix.

Solutions:
(a)
From AT A = I, and Cond2 (AT A) = (Cond2 (A))2 (P67, Property (IV)), we get (Cond2 (A))2 =
Cond2 (AT A) = Cond2 (I) = 1, so Cond2 (A) = 1.
(b)
”If”: For A = aO(a is a scalar and O is an orthogonal matix), Cond2 (A) = ||A||2 ||A−1 ||2 =
a||O||2 a1 ||O−1 ||2 = ||O||2 ||O−1 ||2 = 1 .
”Only if”: Without loss of general, we assume ||A||2 = 1, and now we only need to prove A to be
orthogonal. ||y||2 = ||AA−1 y||2 ≤ ||A||2 ||A−1 y||2 ≤ ||A||2 ||A−1 ||2 ||y||2 = Cond2 (A)||y||2 = ||y||2 , so
all the inequalities must be equanlities. If we take x = A−1 y, then ||Ax||2 = ||A||2 ||x||2 = ||x||2 .
And this is true for any x. So A is norm-preserving and hence orthogonal. (To prove the last argue-
ment, we treat the L2 norm like inner product: (x, x) = ||x||2 = ||Ax||2 = (Ax, Ax) = (AT Ax, x).)

Problem 6.18
||δx||
Prove ||x+δx|| ≤ Cond(A) ||∆A||
||A|| , where Ax = b and (A + ∆A)(x + δx) = b.
Verify the above inequality for the system
    
1 1/2 1/3 x1 1
 1/2 1/3 1/4   x2  =  1 
1/3 1/4 1/5 x3 1
APMA1170 - Homework 4

with  
0 0 0.00003
∆A =  0 0 0 .
0 0 0

Solution:
Substracting Ax = b and (A + ∆A)(x + δx) = b, we get δx = −A−1 ∆A(x + δx). Then we take the
||A||
norm on both hand sides.||δx|| = ||A−1 ∆A(x+δx)|| ≤ ||A−1 ||||∆A||||x+δx|| =||A−1 || ||A|| ||∆A||||x+
δx||=Cond(A) ||∆A|| ||δx|| ||∆A||
||A|| ||x + δx||. And then ||x+δx|| ≤ Cond(A) ||A|| .
For the matrix above, the solution is x=
   
x1 3
 x2  =  −24  .
x3 30

While for the perturbed problem, the solution is x + δx=

 
2.991907283444917
 −23.967629133779667  .
29.973024278149719
||δx||
||x+δx|| = 0.001114657219952 ≤ 0.011163453834773 = Cond(A) ||∆A||
||A|| .

Problem M4.8
the purpose of this exercise is to test that the Hilbert matrix is ill-conditioned with respect to
solving the linear system problem.
(i) creat A=hilb(10). Perturb the (10,1) entry of A by 10−5 . call the perturbed matrix B. Let
b = rand(10, 1). compute x = A \ b, y = B \ b. compute ||x − y|| and ||x−y||
||x|| . what conclusion can
you draw from this?
(ii) compute the condition numbers of both A and B: Cond(A), Cond(B).
(iii) compute the condition number of a using the MATLAB command Cond(A) and then use
it to compute the theoretical upper bound given in theorem 4.23.( ||δx|| ||∆A||
||x|| ≤ Cond(A) ||A|| /(1 −
Cond(A) ||∆A||
||A|| )) Compare this bound with the catural relative error.

Solution:
(i) Here is the computation:

A=hilb(10);
B=A;
B(10,1)=B(10,1)+.00001;
APMA1170 - Homework 4

b=rand(10,1)
b =
0.814723686393179
0.905791937075619
0.126986816293506
0.913375856139019
0.632359246225410
0.097540404999410
0.278498218867048
0.546881519204984
0.957506835434298
0.964888535199277
x=A\b;
y=B\b;
norm(x-y)
ans =
1.248179182526498e+012
ans/norm(x)
ans =
1.326285038669488

||x − y|| = 1.248179182526498e + 012 and ||x−y||

||x|| = 1.326285038669488, which is big for the relative
error. So, we can say Hilbert matrix is ill-conditioned with respect to solving the linear system
problem.
(ii)

Cond(A)
ans =
1.602466539329798e+013
Cond(B)
ans =
4.754875807268704e+012

(iii)

Cond(A)*0.00001/norm(A);
ans/(1-ans)
ans =
-1.000000010932644

Because κ(A) ∗ ||∆A||k|A|| > 1, this negative result can happen. However, the assumption of the
result is that κ(A) ∗ ||∆A||k|A|| < 1. So if it is negative, you are violating the assumptions of the
theorem and all bets are off.
APMA1170 - Homework 4

Problem M6.11
(a) write MATLAB programs called linsyswp and linsyspp to solve Ax=b and to compute the
growth factor(gf) using Caussian elimination with no and partial pivotings, respectively, as follows:
(x̂,gf) = linsyswp(A,b), (x̂,gf) = linsyspp(A,b).
(b) Using the computed solutions and the growth factors obtained in (a) make the following table
for each of the given data set. The function linsyswf is available in MATCOM.
Test data for Problem M6.1: Each of the following matrices of order 10 and 20: Hilbert, Per,
vandermonde and a randomly generated matrix(using rand(n)). For the Pei matrix, take α close
to 1.
Create the vector b in each case such that the solution vector x is a vector with all components
equal to 1. Present your results using Table 6.2.

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp
linsyspp
linsyswf
A−1 b

Solution:
(a) In the function linsyswp and linsyspp, lugewp and lugepp are called. Here are the codes:

function [ x,f ] = linsyswp( A )

%Get Gaussian Elimination without pivoting
%with fucntion linswp:A=LU and its growth factor
%Solve the linear system Ax=b with this factorization
%Input: A - matrix;
%(Create the vector b s.t x is a vector with all 1s)
%Outputs: x-computed solution; f(1)-infinity norm of x;
%f(2)-relative error;f(3)-residual;f(4)-Posterior error bound;
%f(5) - growth factor

[n,m]=size(A); e=ones(n,1); b=A*e; [L,U,f(5)]=lugewp(A);

y=zeros(n,1); x=zeros(n,1);

%---solve forward Ly=b-------------------------------

for i=1:n
sum=0;
if (i~=1)
for j=1:(i-1)
sum=sum+L(i,j)*y(j);
APMA1170 - Homework 4

end
end
y(i)=(b(i)-sum)/L(i,i);
end
%---solve backward Ux=y------------------------------
for i=n:-1:1
sum=0;
if(i~=n)
for j=(i+1):n
sum=sum+U(i,j)*x(j);
end
end
x(i)=(y(i)-sum)/U(i,i);
end
f(1)=norm(x,inf); f(2)=norm(x-e)/norm(e); f(3)=norm(b-A*x);
f(4)=Cond(A)*f(3)/norm(b); end

function [ x,f ] = linsyspp( A )

%Get Gaussian Elimination with partial pivoting
%with fucntion linswp:PA=LU and its growth factor
%Solve the linear system Ax=b with this factorization
%Input: A - matrix;
%(Create the vector b s.t x is a vector with all 1s)
%Outputs: x-computed solution; f(1)-infinity norm of x;
%f(2)-relative error;f(3)-residual;f(4)-Posterior error bound;
%f(5) - growth factor

[n,m]=size(A); e=ones(n,1); b=A*e; [L, U, f(5), P, s, r] = lugepp( A

); b=P’*b; y=zeros(n,1); x=zeros(n,1);

%---solve forward Ly=b-------------------------------

for i=1:n
sum=0;
if (i~=1)
for j=1:(i-1)
sum=sum+L(i,j)*y(j);
end
end
y(i)=(b(i)-sum)/L(i,i);
end
%---solve backward Ux=y------------------------------
APMA1170 - Homework 4

for i=n:-1:1
sum=0;
if(i~=n)
for j=(i+1):n
sum=sum+U(i,j)*x(j);
end
end
x(i)=(y(i)-sum)/U(i,i);
end
f(1)=norm(x,inf); f(2)=norm(x-e)/norm(e); f(3)=norm(b-A*x);
f(4)=Cond(A)*f(3)/norm(b);
end

(b)
(1) For 10 × 10 Hilbert Matrix, we have a table:

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp 1.001191 0.000633 7.85E-16 0.002671 1
linsyspp 9.4E+12 4.81E+12 3.733266 1.27E+13 1
linsyswf 1.000307 0.000163 5.87E-16 0.001999
A−1 b 1.000517 0.000274 2.48E-16 0.000845

(2) For 20 × 20 Hilbert Matrix, we have a table:

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp 145.2 56.42734 2.79E-15 6.66E+02 1
linsyspp 6.81E+17 2.94E+17 22.9138 5.46E+18 1
linsyswf 29.78243 12.51735 1.78E-15 425.1852
A−1 b 100.4214 44.98435 4.07E-15 971.5208

(3) For 10 × 10 Pei Matrix(α = 0.9999), we have a table:

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp 1 7.00E-16 5.62E-15 1.78E-15 1
linsyspp 1 7.00E-16 5.62E-15 1.78E-15 1
linsyswf 1 1.76E-16 0 0
A−1 b 1 1.16E-15 5.33E-15 1.69E-15
APMA1170 - Homework 4

(4) For 20 × 20 Pei Matrix(α = 0.9999), we have a table:

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp 1 1.22E-15 1.70E-14 3.81E-15 1
linsyspp 1 1.22E-15 1.70E-14 3.81E-15 1
linsyswf 1 2.87E-16 5.02E-15 1.12E-15
A−1 b 1 9.56E-16 8.70E-15 1.95E-15

(5) For 10 × 10 Vandermonde Matrix(v = linspace(1, 1.9, 10)), we have a table:

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp 1.000017 9.18E-06 3.75E-13 0.00016 2.003853
linsyspp 1.71E+11 9.16E+10 1600.34 6.82E+11 1
linsyswf 1.000002 9.78E-07 1.42E-13 6.07E-05
A−1 b 1.000001 4.06E-07 1.31E-13 5.59E-05

(6) For 20 × 20 Vandermonde Matrix(v = linspace(1, 2.9, 20)), we have a table:

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp 55.64989 24.59787 0.000155 3.75E+11 15.87062
linsyspp 1.25E+22 5.34E+21 1.22E+11 2.96E+26 1
linsyswf 39562.56 19344.97 7.78E-07 1.88E+09
A−1 b 41745.77 19026.02 2.09E-06 5.04E+09

(7) For 10 × 10 Random Matrix, we have a table:

APMA1170 - Homework 4

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp 1 6.01E-15 1.78E-15 1.33E-14 3.861603
linsyspp 13.63895 7.050191 4.272297 32.10087 1
linsyswf 1 2.75E-15 2.81E-15 2.11E-14
A−1 b 1 7.77E-15 2.55E-15 1.92E-14

(8) For 20 × 20 Random Matrix, we have a table:

Method Norm of the computed Relative error Residual Posterior error Bound Growth
||x−x̂||2
solution ||x̂||∞ ||x||2 ||b − Ax̂||2 Cond2 (A) × ||b−Ax̂||
||b||2
2
Factor
linsyswp 1 2.09E-14 7.51E-14 2.00E-13 36.17695
linsyspp 92.85163 43.1708 58.43972 155.4415 1
linsyswf 1 6.84E-15 8.14E-15 2.17E-14
A−1 b 1 7.85E-15 7.99E-15 2.13E-14