MECH20042/AERO20542 Numerical Methods and Computing
Laboratory exercise 1: Direct methods for the solution of
tridiagonal systems of linear equations
Solution of systems of linear equations is one of the most frequently encountered problems in
numerical modelling and simulation. Efficient numerical methods, both in terms of the execution time
and memory storage are essential to complete this task. Sparse systems of linear equations arise in
many applications, such as finite element or finite volume solution of differential equations. Sparse
linear systems have coefficient matrices that are sparse, i.e., a large proportion of the elements are
equal to zero. Banded matrices are a special class of sparse matrices in which the non-zero coefficients
are concentrated about the main diagonal.
Storing sparse matrices in computer memory as two-dimensional arrays is inefficient, as many zero
elements are kept needlessly in computer memory. Banded matrices can be stored by their diagonals,
where each diagonal is stored as a one-dimensional array (a vector). With this setup a tridiagonal
matrix 𝑇 of size 𝑛 × 𝑛
can be stored using three vectors as follows:
𝐴 = [𝑎11 𝑎22 ⋯ 𝑎𝑛𝑛]
𝑇 ∈ 𝑅
𝑛
,
w**; = [𝑎21 𝑎** ⋯ 𝑎𝑛,𝑛−1]
𝑇 ∈ 𝑅
𝑛−1
,
𝐶 = [𝑎12 𝑎23 ⋯ 𝑎𝑛−1,𝑛]
𝑇 ∈ 𝑅
𝑛−1
.
The Gaussian elimination technique applied to a tridiagonal system 𝑇𝒙 = 𝒇 is particularly simple,
because only the non-zero elements in the sub-diagonal held in vector w**; need to be eliminated. This
algorithm, known as the Thomas algorithm, proceeds as follows:
FORWARD ELIMINATION BACKSUBSTITUTION
𝑎𝑖𝑖 = 𝑎𝑖𝑖 −
𝑎𝑖,𝑖−1
𝑎𝑖−1,𝑖−1
𝑎𝑖−1,𝑖 w**9;𝑛 =
𝑓𝑛
𝑎𝑛𝑛
𝑓𝑖 = 𝑓𝑖 −
𝑎𝑖,𝑖−1
𝑎𝑖−1,𝑖−1
𝑓𝑖−1 w**9;𝑖 =
1
𝑎𝑖𝑖
(𝑓𝑖 − 𝑎𝑖,𝑖+1 w**9;𝑖+1)
𝑖 = 2, … , 𝑛 𝑖 = 𝑛 − 1, … ,1
TASK 1. Calculate the number of arithmetic operations that are required to solve a tridiagonal system
𝑇𝒙 = 𝒇 of size 𝑛 using the Thomas algorithm. Based on this result, determine the asymptotic
complexity of the Thomas algorithm, and compare it to the asymptotic complexity of the standard
Gaussian elimination.
TASK 2. Rewrite the Thomas algorithm in terms of the arrays 𝐴,w**;, and 𝐶 introduced to store the matrix
𝑇 efficiently.
TASK 3. Implement the Thomas algorithm from TASK 2 as a Python function. The input parameters to
the function should be the coefficient matrix 𝑇 (stored as three arrays 𝐴,w**;, and 𝐶) and the right-hand
side vector 𝒇. The output should be the solution vector 𝒙. The coefficient matrix and the right-hand
side should be defined in the main script and passed to the function that solves the system.
TASK 4. Test your code by solving the linear system of size 𝑛 = 10 with the values 𝐴 = 2, and w**; = 𝐶 =
−1. Set the right-hand side to 𝒇 = 𝟏. To verify the correctness of your code, compare the solution
vector obtained from the Thomas algorithm to that obtained by applying the direct solver
numpy.linalg.solve(). For the latter, the coefficient matrix should be assembled.
TASK 5. Solve five linear systems 𝑇𝒙 = 𝒇 with 𝐴 = 2, w**; = 𝐶 = −1 and 𝒇 = 𝟏 varying the problem size
𝑛 between 106
and 108
. Record the execution times in seconds for each case. To accomplish this task,
explore the Python function timer() from the package timeit (refer to the code for matrix
multiplication covered in lectures). Plot a graph where the obtained execution times are represented
as the function of the problem size 𝑛. What are your conclusions about the cost of the Thomas
請加QQ：99515681 郵箱：99515681@qq.com WX：codehelp
