↜ Back to index Introduction to Numerical Analysis 2
Lecture a2
Solutions of exercises from the previous lecture
Gaussian elimination (ガウスの消去法)
Our goal in this class is to solve linear systems (1次連立方程式) of the type \begin{aligned} a_{11} x_1 + a_{12} x_2 + \cdots + a_{1n} x_n &= b_1,\\ a_{21} x_1 + a_{22} x_2 + \cdots + a_{2n} x_n &= b_2,\\ &\ \ \vdots\\ a_{n1} x_1 + a_{n2} x_2 + \cdots + a_{nn} x_n &= b_n,\\ \end{aligned} where n \in {\mathbb{N}} is the number of equations, a_{11}, a_{12}, …, a_{nn} and b_1, b_2, …, b_n are given numbers and x_1, x_2, …, x_n are the unknowns.
This can be written more compactly as Ax = b, where A is the n\times n matrix \begin{pmatrix} a_{11} & a_{12} &&\cdots & a_{1n}\\ a_{21} & a_{22} & && \vdots\\ a_{31} & a_{32} & a_{33} & \\ \vdots& & & \ddots &\\ a_{n1} & a_{n2} & \cdots & & a_{nn}\\ \end{pmatrix} and x, b \in {\mathbb{R}}^n are the vectors \begin{aligned} b &= (b_1, b_2, \ldots, b_n)\\ x &= (x_1, x_2, \ldots, x_n) \end{aligned}
It is easy to do this if A is upper- or lower triangular.
We therefore need to convert the system Ax=b into a form with a triangular matrix. We can do this rather easily by performing row operations on A and b. Indeed, if we add a multiple of one equation to another equation in the system, we do not change the solution. This is equivalent to adding a multiple of a row of the matrix A to another row, and adding the multiple of the component of b in the same row to the component of b in the other row.
For the following A and b, \tag{1} \begin{pmatrix} 1 & 3 & 1 \\ 1 & 1 & -1\\ 3 & 11 & 6 \end{pmatrix}, \begin{pmatrix} 9 \\ 1 \\ 34 \end{pmatrix}, we can add (-1) \times {\rm first\ row} to the second row, and (-3) \times {\rm first\ row} to the third row and we get \tag{2} \begin{pmatrix} 1 & 3 & 1 \\ 0 & -2 & -2\\ 0 & 2 & 3 \end{pmatrix}, \begin{pmatrix} 9 \\ -8 \\ 7 \end{pmatrix} We can then add the second row to the third row and obtain \tag{3} \begin{pmatrix} 1 & 3 & 1 \\ 0 & -2 & -2\\ 0 & 0 & 1 \end{pmatrix}, \begin{pmatrix} 9 \\ -8 \\ -1 \end{pmatrix} System with this matrix and right hand side is now easy to solve.
This algorithm is called Gaussian elimination (ガウスの消去法). The following is a basic implementation in Fortran.
implicit none
integer, parameter :: ms = 20 ! maximal allowed n (number of equations)
real a(ms,ms), b(ms), x(ms) ! arrays for the matrix A, vectors b, x
integer i, j, k, n
real c, s
! read data
! number of equations
read *, n
! read a in column-major order
! the order of loops is important
do j = 1, n
do i = 1, n
read *, a(i, j)
end do
end do
! read b
do i = 1, n
read *, b(i)
end do
! now we are ready to proceed solving Ax = b
! eliminate columns 1, ..., n-1
do k = 1, n - 1
if (a(k,k) == 0) then
print *, "error: a(",k , ", ", k, ") is zero"
stop
end if
do i = k + 1, n
= -a(i,k) / a(k,k)
c do j = k, n
= a(i, j) + c * a(k, j)
a(i,j) end do
= b(i) + c * b(k)
b(i) end do
end do
! add this point time the matrix stored in `a` is a upper triagonal matrix
! compute the solution x of Ax = b
do k = 1, n
! we have to go over rows in reverse
= n + 1 - k
i = 0
s do j = i+1,n
= s + a(i, j) * x(j)
s end do
= (b(i) - s) / a(i, i)
x(i) end do
! print the vector x
do i = 1, n
print *, x(i)
end do
end
You can find versions of this code online on many websites. Our goal is to understand it and develop the skills to write it ourselves based on the explanation of the algorithm and modify it if necessary.
The above code for example fails for a simple matrix A = \begin{pmatrix} 0 & 1 \\ 1 & 1 \end{pmatrix}, even though this matrix is invertible and Ax = b has a solution. We will learn how to fix the program later.