Lab Session 4
MA571: Numerical Linear Algebra January-May, 2026 S. Bora
Important instructions:
(i) Switch to format long e for all experiments.
(ii) Submit a single livescript program that contains all comments, answers and codes
necessary to produce the required outputs. Ensure that the answers are correctly
numbered and the file does not include any irrelevant material. The livescript
program should be saved as [Link]
(iii) Read the note before attempting the questions 2 to 4.
Note : The Rule-of-thumb of ill-conditioning says that if
(i) κ(A) ≈ 10t ;
(ii) the computed solution of Ax = b has been obtained from a backward stable algorithm;
(iii) the entries of A and b are accurate to about s decimal places where s > t;
then one should expect an agreement of at least s − t decimal places between the corresponding
entries of the exact solution and the computed solution of the system.
1. The purpose of this experiment is to note the extreme ill-conditioning of the Hilbert matrix.
Convince yourself that the condition number of H grows quickly with n. Try
>> C=[];
>> N= 2:2:16;
for n=N
H=hilb(n); C=[C; cond(H)];
end
>> semilogy(N,C)
Based on this graph formulate a conjecture as to the relationship between cond(H) and n.
[Note: The Matlab command cond(H) computes the 2-norm condition number of H. Type
help cond for details.]
Modify the above code and compute the condition number in 1-norm and ∞-norm.
2. Given a Hilbert matrix H, the aim of this‘ exercise is to verify whether the rule of thumb for
solutions of the system Hx = b obtained via a number of different solution methods holds. In
practice, the exact solution is unavailable. But for the experiments we can use the following
trick to gain access to it:
Choose x first and set b := Hx and then solve Hx = b by setting x1 = H\b. Then x is the
exact solution of Hx = b and x1 is the computed solution!
The system Hx = b may be solved by using the invhilb command that generates the inverse
of a Hilbert matrix. Moreover you can also use your gepp function formulated in the first
problem to solve Hx = b. You may have to use format long e to see more digits. Note
1
that H \ b will NOT find a solution via GEPP in this case as H is positive definite. Instead
it will use Cholesky method to find the solution. Now execute the following steps:
>> n=8;
>> H=hilb(n); HI = invhilb(n);
>> x= rand(n,1);
>> b =H*x;
>> x1 = H\ b; x2 = HI*b;
% obtain x3 by solving Hx=b using GEPP.
>> [x x1 x2 x3]
>> [cond(H) norm(x-x1)/norm(x) norm(x-x2)/norm(x) norm(x-x3)/norm(x)]
Repeat for n = 10 and n = 11 and list the results corresponding to n = 8, 10, 12, and
determine correct digits in x1, x2, x3. Now answer the following questions:
How many digits are lost in computing x1, x2 and x3? Does the loss of accuracy in each
case agree with the value predicted by the Rule-of-thumb mentioned earlier? Note that since
the experiments are performed on randomly generated data and their is no error other than
rounding error, you should take s ≈ 16 (wait for a few more classes to know why!) in the
Rule-of-thumb analysis.
Which is better among x1, x2 and x3 or isn’t there much of a difference?
3. The purpose of this experiment is to illustrate that a small value of the norm of the residual
alone is not enough to guarantee an accurate answer.
If x̂ is the computed solution of Ax = b then r := Ax̂ − b is called the residual. Of course
r = 0 if and only if x = x̂. But usually r ̸= 0. Does a small ∥r∥/∥b∥ imply ∥x − x̂∥/∥x∥ small?
Try the following:
>> n=10;
>> H=hilb(n); x = randn(n,1);
>> b = H*x;
>> x1= H \ b;
>> r = H*xt-b;
>> disp( [norm(r)/norm(b) norm(x-xt)/norm(x)])
What is your conclusion?
4. The purpose of this experiment is two-fold. Firstly, the aim is to show that only the coefficient
matrix A having a small condition number is not enough to ensure accuracy in the computed
solution of the system Ax = b. Secondly, the aim is to understand the role of the assumption
of backward stability of the algorithm for the predictions from Rule-of-thumb analysis to
hold.
Recall the Wilkinson’s matrix that you had generated in one of your previous classes. For
n = 32, pick a random x and then compute b := W ∗ x. Store the solution obtained by using
GEPP in x̂ and compute the (forward) error ∥x − x̂∥∞ /∥x∥∞ . Also compute cond(A) in the
infinity norm and ∥r∥∞ /∥b∥∞ . Repeat the test for n = 64.
Further repeat the above experiment using QR decomposition (More about this later but find-
ing it is easy in MATLAB. Just typing [Q, R] = qr(A) gives unitary Q and upper triangular R
such that A = QR.) Solve W x = b using QR decomposition (using x = colbackward(R,Q’*b)).
Tabulate all the quantities for the different experiments. Noting that s ≈ 16 in the rule-of-
thumb analysis, answer the following.
2
(a) Which of the two methods appear to give a lower error in the computed solution?
(b) For which of the two methods is the Rule-of-thumb predicting the correct answer rea-
sonably well?
(c) Which of the two methods give rise to a lower value of ∥r∥∞ /∥b∥∞ ?
(d) What can you say about the backward stability of GEPP and QR decomposition methods
from the experiments?