4.3 Gaussian Elimination#
In the lecture we learned that we can solve a linear equation using the Gaussian elimination method. The idea is to transform the linear equation into an upper (or lower) triangular matrix. We did this by transforming each row of the matrix. We want to solve for \(x\) in the following equation:
As an example let’s look at the following equation:
We can rewrite the matrix by extending it with the vector \(b\) into \([A b]\):
We can now perform row operations to transform the matrix into an upper triangular matrix. The first step is to subtract the first row multiplied by \(4/1\) from the second row. This will result in the following matrix:
We can now solve the second equation for \(x_2\) and substitute it into the first equation. This will result in the following equation:
We can now solve the first equation for \(x_1\) and get the solution \(x_1=1\) and \(x_2=2\).
As a more general way to show our example of a 2x2 matrix \(A\) we have:
To perform the Gaussian elimination of this linear equation we can use the following algorithm:
\(\alpha = a_{21}/a_{11}\)
\(a_{22} = a_{22} - \alpha * a_{12}\)
\(b_2 = b_2 - \alpha * b_1\)
\(x_2 = b_2/a_{22}\)
\(b_1 = b_1 - a_{12} * x_2\)
\(x_1 = b_1/a_{11}\)
return \(x_1\) and \(x_2\)
A = [1 2; 3 4]
b = [5; 6]
A_ = float(copy(A)) # we copy A since changing an element of A will change A as well
b_ = float(copy(b))
α = A_[2,1] / A_[1,1]
A_[2,:] = A_[2,:] - α * A_[1,:]
b_[2] = b_[2] - α * b_[1]
x_2 = b_[2] / A_[2,2]
x_1 = (b_[1] - A_[1,2] * x_2) / A_[1,1]
x = [x_1; x_2]
2-element Vector{Float64}:
-4.0
4.5
We can test our implementation by estimating the error of the solution:
The error should be zero.
A*x - b
2-element Vector{Float64}:
0.0
0.0
Can you think of a way to generalize this algorithm to a square matrix of arbitrary size? One difference is that we have to perform the algorithm for each row. Each row will have a different \(\alpha\) value. We can use a loop to perform the algorithm for each row. In the example above we have applied the subtraction of the first row multiplied by \(4/1\) from the second row. We can generalize this by subtracting the first row multiplied by \(a_{21}/a_{11}\) from the second row. Also note that we can actually substract the whole row at once.
The task is to write a function that implements the Gaussian elimination method for arbitrary square matrices. The function should take the matrix \(A\) and the vector \(b\) as input and return the solution vector \(x\). We will divide this task into two parts. First we will write a function that transforms the linear equation into an upper triangular matrix. Then we will write a function that performs the back substitution to get the solution vector \(x\).
Let’s start by defining a new matrix \(A\) and vector \(b\):
A = [1 2 3; 4 5 6; 7 8 10]
b = [1, 2, 3]
3-element Vector{Int64}:
1
2
3
Now we can write a function that transforms the matrix \(A\) into an upper triangular matrix. The function should take the matrix \(A\) and the vector \(b\) as input and return the transformed matrix \(A\) and the transformed vector \(b\).
function gauss_elimination(A::Matrix{T}, b::Vector{T}) where T <: Number
n = size(A, 1)
U = float(copy(A))
b_ = float(copy(b))
for k = 1:n-1
for i = k+1:n
factor = U[i,k] / U[k,k]
U[i,:] = U[i, :] - (factor * U[k,:])
b_[i] = b_[i] - (factor * b_[k])
end
end
return U, b_
end
gauss_elimination (generic function with 1 method)
Let’s test the function with the matrix \(A\) and the vector \(b\) defined above. We should get an upper triangular matrix.
U, b_ = gauss_elimination(A, b)
([1.0 2.0 3.0; 0.0 -3.0 -6.0; 0.0 0.0 1.0], [1.0, -2.0, 0.0])
Now we can write a function that performs the back substitution to get the solution vector \(x\). The function should take the matrix \(U\) and the transformed vector \(\hat{b}\) as input and return the solution vector \(x\).
function backsub(U,b)
n = size(U,1)
x = zeros(n)
for i = n:-1:1
x[i] = b_[i]
for j = i+1:n
x[i] -= U[i,j] * x[j]
end
x[i] /= U[i,i]
end
return x
end
backsub (generic function with 1 method)
We can now test the function with the matrix \(U\) and the transformed vector \(\hat{b}\) defined above. We should get the solution vector \(x\).
x = backsub(U,b_)
3-element Vector{Float64}:
-0.33333333333333326
0.6666666666666666
0.0
A*x-b
3-element Vector{Float64}:
0.0
0.0
0.0