LAB MANUAL

NUMERICAL LINEAR ALGEBRA

(UMA021)
B.E. Second Year
Department of Mathematics
Thapar Institute of Engineering and Technology, Patiala

Name:………………………….
Roll No.:……………………….
Group:…………………………
Instructor:………………………
 Contents
S.No. Experiment Date Signature
of
Instructor
1. (a) Use intermediate value theorem to find the interval of the roots.
(b) Find the root of non-linear equation f(x) = 0 using bisection
method.
2. (a) Find the root of non-linear equation f(x) = 0 using Newton’s
method.
(b) To solve the system of non-linear equations using Newton’s
method.
3. Find the root of non-linear equation f(x) = 0 using fixed-point
iteration method
4. Implement Lagrange interpolating polynomials of degree ≤ n on
n+1 discrete data points.
5. Implement Newton’s divided difference interpolating polynomials
for n+1 discrete data points.
6. Integrate a function numerically using composite trapezoidal and
Simpson’s rules.
7. Solve system of linear equations Ax = b using Gauss elimination
method.
8. Solve system of linear equations Ax = b using Gauss-Seidel method.
9. Fit a curve for given data points by using principle of least squares.
10.(a) Find a dominant eigen-value and associated eigen-vector by Power
method.
(b) Find all the eigenvalues of a matrix using QR algorithm.
 Experiment 1: Bisection Method

1. Algorithm of Intermediate Value Theorem (IVT): To determine all the subintervals [a, b]
of [-N, N] that containing the roots of f(x) = 0.
Input: function f(x), and the values of h, N
for i = -N : h : N
if f(i) * f(i + h) < 0 then a = i and b = i + h
end if
end i
2. Algorithm of Bisection Method: To determine a root of f(x) = 0 that is accurate within a
specified tolerance value ϵ, given values a and b such that f(a) * f(b) < 0.
Define c = (a + b)/2.
if f(a) * f(c) < 0, then set b = c, otherwise a = c.
end if.
Until |a – b| ≤ ϵ (tolerance value).
Print root as c.
Stopping Criteria: Since this is an iterative method, we must determine some stopping
criteria that will allow the iteration to stop. Criteria |f(ck)| very small can be misleading since it
is possible to have |f(ck)| very small, even if ck is not close to the root.
The interval length after N iterations is (b – a)/2N. So, to obtain an accuracy of ϵ, we must
have
log ( b − a ) − log ∈
N≥ .
log 2

3. Students are required to write both the programs (IVT and Bisection) and implement it on the
following examples.
(i) Use bisection method in computing of √ 29 with ϵ = 0.001, N = 10, h = 1.
(ii) Determine the number of iterations necessary to solve f(x) = x3 + 4x2 – 10 = 0 with
accuracy 10-3 using a = 1 and b = 2 and hence find the root with desired accuracy.
4. Thermistors are temperature-measuring devices based on the principle that the thermistor
material exhibits a change in electrical resistance with a change in temperature.

By measuring the resistance of the thermistor material, one can then determine the
temperature. For a 10K3A Betatherm thermistor, the relationship between the resistance R of
the thermistor and the temperature is given by
 1 3
=1.129241 ×10 +2.341077 ×10 ln ( R ) +8.775468 ×10 { ln ( R ) } where T is in
−3 −4 −8
T
Kelvin and R is in ohms. Use the bisection method to find the resistance R at 18.99 0C.
Experiment 2: Newton’s Method for root finding and solution of system of non-linear equations

1. Algorithm for Newton’s method: Find a solution to f(x) = 0, given an initial approximation x0.
Input: Initial approximation x0, tolerance value ϵ, maximum number of iterations N.
Output: Approximate solution or message of failure.
Step 1: Set i = 1.
Step 2: While i ≤ N do Steps 3 to 6.
f ( x0)
Step 3: Set x 1= x 0 − . (Compute xi).
df ( x 0 )
|x1 − x 0|
Step 4: If |x1 − x 0|≤ ϵ or ≤ ϵ then OUTPUT x1; (The procedure is successful)
|x 1|
STOP.
Step 5: Set i = i + 1.
Step 6: Set x0 = x1. (Update x0)
Step 7: Print (‘The method failed after N iterations, N =’, N); (The procedure is unsuccessful)
STOP

2. Students are required to write both the program and implement it on the following examples.
Take tolerance value ϵ = 0.00001
(i) Compute √ 17 .
(ii) The root of exp(– x)(x2 + 5x + 2) + 1 = 0. Take initial guess –1.0.
(iii) Find a non-zero solution of x = 2sin x. (Apply IVT to find an initial guess)

3. An oscillating current in an electric circuit is described by i=9 e −t sin ⁡( 2 πt ) , where t is in

seconds. Determine the lowest value of t such that i = 3.5.

4. Algorithm to solve the system of non-linear equations using Newton’s method.

To approximate the solution of the nonlinear system F(x) = 0 given an initial approximation x:
INPUT number n of equations and unknowns; initial approximation x = (x 1, x2, …, xn)t , tolerance
TOL; maximum number of iterations N.
OUTPUT approximate solution x = (x 1, x2, …, xn)t or a message that the number of iterations was
exceeded.

Step 1 Set k = 1.

Step 2 While (k ≤ N) do Steps 3-7.

Step 3 Calculate F(x) and J(x), where J(x)i,j = (∂fi(x)/∂xj) for 1 ≤ i, j ≤ n.

Step 4 Solve the n x n linear system J(x)y = - F(x).

Step 5 Set x = x + y.

Step 6 If ||y|| < TOL then OUTPUT (x);

(The procedure was successful)
STOP.

Step 7 Set k = k+1.

Step 8 OUTPUT (‘Maximum number of iterations exceeded’);

(The procedure was unsuccessful)
STOP.

5. (a) For the given system of non-linear equations:

1
3 x1 − cos ( x 1 x 2 ) − =0 ,
2
2 2
x 1 − 81 ( x 2 +0.1 ) +sin ( x 3 )+1.06=0 ,
−x x 10 π −3
e 1 2
+20 x 3 + =0 ,
3
find the approximate solution (0.5, 0, -0.52359877) t to this problem with initial guess
x(0)= (0.1, 0.1, -0.1)t.

(b) For the given system of non-linear equations:

x2 + y2 -1 = 0,
x2 + y2 = 0,

find the approximate solution to this problem with initial guess

x(0)= (1, 1)t.
 Experiment 3: Fixed-point Iteration Method

1. Algorithm for Fixed-point iteration method: To find a solution to x = g(x), given an initial
approximation x0.
Input: Initial approximation x0, tolerance value ϵ, maximum number of iterations N.
Output: Approximate solution or message of failure.
Step 1: Set i = 1.
Step 2: While i ≤ N do Steps 3 to 6.
Step 3: Set x 1= g ( x 0 ) . (Compute xi).
|x1 − x 0|
Step 4: If |x1 − x 0|≤ ϵ or ≤ ϵ then OUTPUT x1; (The procedure is successful)
|x 1|
STOP.
Step 5: Set i = i + 1.
Step 6: Set x0 = x1. (Update x0)
Step 7: Print the output and STOP.

2. The equation f(x) = x3 + 4x2 – 10 = 0 has a unique root in [1,2]. There are many ways to
change the equation to the fixed-point form x = g(x) using simple algebraic manipulation. Let
g1, g2, g3, g4 and g5 are iteration functions obtained by the given function, then check which of
the following iteration functions will converge to the fixed point? (Tolerance ϵ = 10-3)
3 2
(a) g1 ( x )= x − x − 4 x +10

(b) g2 ( x )=
√10
x
−4 x

√
(c) g3 ( x )=0.5 10− x 3

(d) g 4 ( x )=
√
10
4+x
3 2
x +4 x − 10
(e) g5 ( x )= x − 2
3 x +8 x

3. Find the smallest and second smallest positive roots of the equation tan(x) = 4x, with an
accuracy of 10-3 using fixed-point iterations.

4. Use a fixed-point iteration method to determine a solution accurate to within 10-2 for
2 sin πx + x=0 on [1,2]. Use initial guess x0 = 1.
 Experiment 4: Lagrange interpolation

1. Algorithm for Lagrange interpolation: Given a set of function values

x x1 x2 …… xn
f(x) f(x1) f(x2) …… f(xn)

To approximate the value of a function f(x) at x = p using Lagrange’s interpolating

polynomial P n−1 ( x ) of degree ≤ n – 1, given by
P n−1 ( x ) = l1(x) f(x1) + l2(x) f(x2) +…………. + ln(x) f(xn)
n
( p − x j)
where li ( x )= ∏ at x = p.
j=1 ; i≠ j ( x i − x j )

We write the following algorithm by taking n points and thus we will obtain a polynomial of
degree ≤ n – 1.

Input:The degree of the polynomial, the values x(i) and f(i), i = 1, 2,….,n and the point of
interpolation p.
Output: Value of P n−1 ( p ).
Algorithm:
Step 1. Calculate the Lagrange’s fundamental polynomials li(x) using the following loop:
for i = 1 to n
l(i) = 1
for j = 1 to n
if j ≠ i
p− x j
l ( i )= l (i)
x (i) − x ( j)
end j
end i

Step 2. Calculate the approximate value of the function at x=p using the following loop:
sum = 0
for i = 1 to n
sum = sum + l(i)*f(i)
end i

Step 3. Print sum.

2. The following data define the sea-level concentration of dissolved oxygen for fresh water as a
function of temperature:
t 0 8 16 24 32 40
O(t) 14.621 11.843 9.870 8.418 7.305 6.413

Use Lagrange’s interpolation formula to approximate the value of O(15) and O(27).
3. Generate eight equally-spaced points from the function f(x) = sin2x from x = 0 to 2π. Use
Lagrange interpolation to approximate f(0.5), f(3.5), f(5.5) and f(6.0).

Experiment 5: Newton’s Divided Difference Interpolation

1. Algorithm for Newton’s divided difference interpolation:

Given n distinct numbers x 1 , x 2 ,… , x n and their corresponding function values
f ( x 1 ) , f ( x 2 ) , … , f ( x n ).
Approximate the value of a function f(x) at x = p using Newton’s divided difference
interpolating polynomial P n−1 ( x ) of degree ≤ n – 1.

Input: Enter n the number of data points; enter n distinct numbers x 1 , x 2 ,… , x n ; enter
corresponding function values f ( x 1 ) , f ( x 2 ) , … , f ( x n ) as F 1,1 , F 1,2 , … , F 1 ,n ; enter an
interpolating point p.
Output: the numbers F 2,2 , F 3,3 , … , F n,n such that
n i −1
P n−1 ( p )=∑ F i ,i ∏ ( p − x j )
i=1 j=1

where Fk, k is the (k-1) divided difference f [ x1 , x 2 ,… , x k ]

th

Step 1: Evaluate F 2,2 , F 3,3 ,… , F n, n

for i=2 :n
for j=i : n
F j ,i −1 − F j −1,i −1
Evaluate F j ,i = .
x j − x j − i+1
end j
end i
i −1
Step 2: Evaluate ∏ ( p − x j ) for each i=1 :n
j=1
for i=1 :n
Set product (i) = 1
for j=1 :i − 1
product (i) = product (i) * ( p − x j )
end j
end i
Step 3: Evaluate P n−1 ( p )
Set Sum = 0
for i=1 :n
Sum = Sum + ( F i , i * product (i))
end i
Step 4: OUTPUT Sum ≡ Pn −1 ( p )
STOP
2. The following data represents the function f ( x )=e x .
x 1 1.5 2.0 2.5
f(x) 2.7183 4.4817 7.3891 12.1825

Estimate the value of f(2.25) using the Newton’s divided difference interpolation. Compare
with the exact value.
3. Approximate f(0.43) by using Newton’s divided difference interpolation, construct the
interpolating polynomials for the following data.
f(0) = 1, f(0.25) = 1.64872, f(0.5) = 2.71828, f(0.75) = 4.4816.
 Experiment 6: Trapezoidal and Simpson’s Rules

1. Algorithm for Composite Trapezoidal rule:

Step 1: Input function f(x); end points a and b; and N number of subintervals.
b −a
Step 2: Set h= .
N
Step 3: Set sum = 0
Step 4: for i = 1 to N-1
Step 5: Set x=a+h ∗ i
Step 6: Set sum= sum +2 ∗ f ( x )
end i
Step 7: Set sum= sum +f ( a )+ f ( b )
h
Step 8: Set ans= sum*
2
STOP

2. Algorithm for Composite Simpson’s rule:

Step 1: Input function f(x); end points a and b; and N number of subintervals (even).
b −a
Step 2: Set h= .
N
Step 3: Set sum = 0
Step 4: for i = 1 to N-1
Step 5: Set x=a+h∗ i
Step 6: if rem(i,2) = = 0
sum= sum +2 ∗ f ( x )
else
sum= sum +4 ∗ f ( x )
end if
end i
Step 7: Set sum= sum +f ( a )+ f ( b )

Step 8: Set ans = sum* ( h3 )

STOP

3. Approximate the following integrals using the composite trapezoidal and Simpson rule by
taking different subintervals (e.g. 4, 6, 10, 20)
0.25
(a) I = ∫ 2
( cosx ) dx
− 0.25
e +1
1
(b) I = ∫ dx
e xlnx
1
(c) I =∫ e
2
−x
cosx dx
−1
4. The length of the curve represented by a function y = f(x) on an interval [a, b] is given by the
integral
b

√ 2
I =∫ 1+[ f ' ( x ) ] dx .
a
Use the trapezoidal rule and Simpson’s rule with 4 and 8 subintervals, compute the length of
the curve
y=tan ( 1+ x ) ,0 ≤ x ≤ 2.
−1 2
 Experiment 7: Gauss Elimination Method

1. Algorithm for Gauss elimination method: Find a solution of system of linear equations.
Input: Number of unknowns and equations n,
Augmented matrix A=[ aij ], where 1 ≤ i ≤ n, and 1 ≤ j ≤ n+1.
Output: Solution ( x1 , x 2 ,⋯ , x n ) or message that the linear system has no unique solution.
Step 1: For i=1,2 , ⋯ , n −1 do Steps 2 – 4. (Elimination process)
Step 2: Let p be the smallest integer with i ≤ p ≤ n and a pi ≠ 0.
If no integer p can be found then
OUTPUT (‘no unique solution exists’); STOP.
Step 3: If p ≠ i then perform ( E p ) ↔ ( Ei ).
Step 4: For j=i+1 ,⋯ , n do Steps 5 and 6.
Step 5: Set m ji=a ji / aii.
Step 6: Perform ( E j − m ji E i ) ↔ ( E j );
Step 7: If a nn=0 then
OUTPUT (‘no unique solution exists’); STOP.
Step 8: Set x n =a n,n+1 / ann. (Start backward substitution)

[ ]
n
Step 9: For i=n −1 , n − 2 ,⋯ ,1 set x i = ai ,n+1 − ∑ aij x j / aii.
j=i+1

Step 10: OUTPUT ( x1 , x 2 ,⋯ , x n ). (Procedure completed successfully)

STOP.

2. Use Gauss elimination method to find the solution of the following linear system of
equations:
10x + 8y – 3z + u = 16
2x + 10y + z – 4u = 9
3x – 4y + 10z + u = 10
2x + 2y – 3z + 10u = 11

3. Solve the following linear system of equations:

π x1 + √ 2 x 2 − x 3 + x 4 =0
e x 1 − x 2 + x 3 +2 x 4 =1
x 1 + x 2 − √ 3 x 3 + x 4=2
− x1 − x 2 + x 3 − √ 5 x 4 =3
4. Kirchhoff’s laws of electrical circuits state that both the net flow of current through each
junction and the net voltage drop around each closed loop of a circuit are zero. Suppose that a
potential of V volts is applied between the points A and G in the circuit and that i1, i2, i3, i4 and
i5 represent current flow as shown in the diagram. Using G as a reference point, Kirchhoff’s
laws imply that the currents satisfy the following system of linear equations:

5 i1 +5 i2 =V
i3 − i4 − i5=0
2 i4 −3 i5 =0
i1 −i2 − i3=0
5 i2 − 7 i3 − 2 i4 =0

Take V = 5.5 and solve the system.

 Experiment 8: Gauss-Seidel Method

1. Algorithm for Gauss Seidel Method: Find a solution of system of linear equations Ax = b.
Input: Number of unknowns n; Coefficient matrix A=[ aij ], where 1 ≤ i ≤ n, and 1 ≤ j ≤ n;
column vector b; Initial solution vector x0; tolerance value tol; maximum number of iterations
N.
Output: Solution ( x1 , x 2 ,⋯ , x n ).
Step 1: For k =1,2 , ⋯ , N do Steps 2 – 4.
Step 2: For i=1,2 , ⋯ , n

[ ]
i− 1 n
1
x i= b −∑ (a x )− ∑ (a x )
aii i j=1 ij j j=i+1 ij 0 j
Step 3: If ‖x − x 0‖<tol then OUTPUT ( x1 , x 2 ,⋯ , x n ).
STOP
Step 4: Set x0 = x. (Update x0)
Step 5: Print OUTPUT ( x1 , x 2 ,⋯ , x n ) (Procedure completed successfully)
STOP.

2. Use Gauss Seidel method to find the solution of the following linear systems with an initial
vector [0,0,0,0] and tolerance value 10-5 in the ‖.‖∞ norm:

(a) 10x + 8y – 3z + u = 16 (b) 4 x1 + x 2 − x3 + x 4 =−2

2x + 10y + z – 4u = 9 x 1 +4 x 2 − x3 − x 4 =− 1
3x – 4y + 10z + u = 10
− x1 − x 2 +5 x 3 + x 4 =0
2x + 2y – 3z + 10u = 11
x 1 − x2 + x 3 +3 x 4 =1

3. Use Gauss Seidel method to solve the following linear system with an initial vector [0,0,0]
and tolerance value 10-3 in the ‖.‖∞ norm:
4.63 x 1 −1.21 x 2 +3.22 x 3=2.22
−3.07 x 1 +5.48 x 2 +2.11 x 3=− 3.17
1.26 x1 +3.11 x 2 +4.57 x 3=5.11
 Experiment 9: Least Square Approximation

1. Write an algorithm for least square approximations to fit any curve of the forms:
B
y=a+bx , y=a+bx+ c x 2, y= A √ x + , …
x

2. Use the method of least squares to fit the linear and quadratic polynomial to the following
data:
x -2 -1 0 1 2
f(x) 15 1 1 3 19

3. By the method of least square fit a curve of the form y=a x b to the following data:
x 2 3 4 5
y 27.8 62.1 110 161

B
4. Use the method of least squares to fit a curve y= A √ x + to the following data:
x
x 0.1 0.2 0.4 0.5 1 2
y 21 11 7 6 5 6
 Experiment 10: QR Algorithm, Gram-Schmidt Process and Power Method
1. Algorithm
This method picks all the eigenvalues. Firstly we write the QR factorization of matrix A using Gram-
Schmidt Process i.e. A = QR.
The QR algorithm is an iterative procedure. To find all the eigenvalues of matrix A, we
generate a sequence of matrices Ak that are orthonormally similar to A and thus have same
eigenvalues.
Basic QR:
We write A = QR and assign A1 = A = Q1R1.
Then for k = 1, 2, … we have Ak = QkRk and Ak+1 = QkRk
end.
The iterates converge to an upper triangular matrix with eigenvalues on main diagonal.

2. Input Matrix (A) of size (m x n) with linearly independent columns.

3.Initialization:
- Get the dimensions of (A) and store in [m,n].
- Initialize matrix (Q) as an (m x n) zero matrix.
- Initialize matrix (R) as an (n x n) zero matrix.

4. Gram-Schmidt Orthogonalization:
- For each column (j) from 1 to (n):
(i). Set (v) to be the (j)-th column of (A).
(ii). For each previous column (i) from 1 to (j-1):
- Compute the inner product of (Q(:, i)) and (v) and store in (R(i, j)).
- Subtract (R(i, j) x Q(:, i)) from (v).

5. Compute the norm of (v) and store in (R(j, j)).

6. Normalize (v) by dividing by (R(j, j)) and store in (Q(:, j)).
7. Output: Orthogonal matrix (Q). Upper triangular matrix (R).
8. Display Results: Print matrices (Q) and (R).

5. Determine the eigenvalues of the following matrix using QR algorithm.

[ ]
1 1 0
A= 1 0 1
0 1 1
6. Algorithm for Power method:
Step 1: START
Step 2: Define matrix A and initial guess x.
Step 3: Calculate y = Ax
Step 4: Find the largest element in magnitude of matrix y and assign to K.
Step 5: Calculate fresh value x = (1/K) * y.
Step 6: If |K ( n ) − K ( n −1 )| > error, go to Setp 3.
Step 7: STOP.

7. Determine the largest eigen-value and the corresponding eigen-vector of the following
matrices using the power method. Use x0 = [1,1,1]T and ϵ = 10-3:

[ ]
4 1 0
(a) 1 20 1 .Use x0 = [1,1,1]T and ϵ = 10-3
0 1 4

[ ]
1 1 0 0
1 2 0 1
(b) . Use x0 = [1,1,0,1]T and ϵ = 10-3
0 0 3 3
0 1 2 3