Solving a System of Linear Equations

Currently, we have only been making use of graphs and plotting the point of intersection to get the solution, but there are ways to do it mathematically. We will be looking at how to obtain these solutions using Substitution and Gaussian Elimination.

Solving a System Using Substitution Method

First, lets take the example of this system of equations.

1. 3x+2y-z=1

2. 2x-2y+4z=-2

3. -x+1/2 y-z=0

We can write x in terms of other variables to substitute it into other equations. Taking from the third equation,


We then substitute this value into the second equation,



We will get an equation of y in terms of z.


Then we proceed to substitute the value of y and x into the first equation.






With this value of z, we can then solve for y and then x by substituting back the values of z.





Therefore, the solution to this system is




Solving a System Using Gaussian Elimination Method

Gaussian elimination is an algorithm for solving systems of linear equations. It transforms matrices, which represent a linear system, into Row Elchleon form, where it can be easily operated on and used to solve a system. To represent a system with a matrix, we follow the format of

\begin{bmatrix}a_{1} & b_{1}& c_{1}&x \\a_{2} & b_{2}& c_{2}&y \\a_{3} & b_{3}& c_{3}&z\end{bmatrix}

where the system is

1. a_{1}x+b_{1}y+c_{1}z=x

2. a_{2}x+b_{2}y+c_{2}z=y

3. a_{3}x+b_{3} y+c_{3}z=z

The columns correspond to the individual variables (e.g. the first column is the coefficients for x, second for y, etc.) and the rows represent the equation(first row represents the first equation, etc.).

Row Elchleon Form

However, a matrix by itself is not usually useful to us. We can simplify the matrix into row elchleon form to make it easier to read and compute the solution. A matrix is in row elchleon form if it fulfills these conditions.

1. All rows consisting of only zeroes are at the bottom.

2. The leading coefficient (also called the pivot) of a nonzero row is always
to the right of the leading coefficient of the row above.

In summary, a matrix in row elchleon form will look something like this:

\begin{bmatrix}a_{1} & b_{1}& c_{1}&x \\0&b_{2}& c_{2}&y \\0&0&c_{3}&z\end{bmatrix}

Reduced Row Elchleon Form

You can further simplfy the matrix into a reduced row elchleon form when the matrix the leading entry in each nonzero row is a 1 (called a leading 1) and the values above a leading entry is 0. The matrix will therefore look like this:

\begin{bmatrix}1 & 0& 0&x \\0&1& 0&y \\0&0&1&z\end{bmatrix}

which will allow us to derive the values of a, b and c by just looking at the equation. However, in some cases, simplfying the matrix into reduced row elchleon form takes a longer amount time to compute as compared to deriving the answer from a matrix that is in row elchleon form.

Transformation Operations

There are mainly three ways to simplfy a matrix into Row Elchleon form, which are

1. Row switching: switch one row with another row
2. Row multiplication: multiply a row by a nonzero constant

3. Row addition/subtraction: add the multiple of another row to a row.

Row Switching would simply switch the position of one row with another. In this case,


We can switch Row 2 and Row 3 to turn the matrix into Row Elchleon Form, where it will look like


Row multiplication can be used to scale the rows by a factor of a non-zero constant,n.

\begin{bmatrix}a_{1} & b_{1}& c_{1}&x \\2a_{2} & 2b_{2}& 2c_{2}&2y \\a_{3} & b_{3}& c_{3}&z\end{bmatrix}

In this example, we scaled the second row by a factor of 2.

Row addition/subtraction can be used to add or subtract values from one row with that of another row. A matrix


can be converted into Row Elchleon form by subtracting Row 3 with 1/2 of Row 2, resulting in

\begin{bmatrix}1&1&1&6\\0&2&1&7\\0&1-1/2*2&2-1/2*1&8-1/2*7\end{bmatrix}= \begin{bmatrix}1&1&1&6\\0&2&1&7\\0&0&1.5&4.5\end{bmatrix}


Lets use an example to showcase turning a system of equations into matrix form, then simplifying to row elchleon form and reduced row elchleon form. The system,

1. 2x+y-z=7

2. 2x-y+2z=6

3. x-2y+z=0

can be represented in the matrix as


This system is now represented as a Matrix. To simplfy it to Row Elchleon Form, we will first need to remove the leading coefficient in Row 2, which is a 2. We will simply perform the operation of Row 2 – Row 1 to turn the leading coefficient in row 2 to a 0.


Next, we will have to turn the leading coefficient in Row 3 to a 0. This can be done by Row 3 – 1/2*Row 1

\begin{bmatrix}2&1&-1&7\\0&-2&3&-1\\1-1&-2-1/2&1-(-1/2)&0-7/2\end{bmatrix} = \begin{bmatrix}2&1&-1&7\\0&-2&3&-1\\0&-2.5&1.5&-3.5\end{bmatrix}

Next, we will remove the leading coefficient from Row 3, which is now -2.5. This is done by Row 3 – 5/4*Row 2

\begin{bmatrix}2&1&-1&7\\0&-2&3&-1\\0&-2.5-(-2*5/4)&1.5-(3*5/4)&-3.5-(-1*5/4)\end{bmatrix} =\begin{bmatrix}2&1&-1&7\\0&-2&3&-1\\0&0&-9/4&-9/4\end{bmatrix}

The matrix is now in Row Elchleon Form. We can also do the same by using the NM Dev class GaussianElimination, which performs elementary row operations to reduce a matrix to the row echelon form.

					%use s2
//var a is the matrix for the linear system
var a = DenseMatrix(arrayOf(
    doubleArrayOf(2.0, 1.0, -1.0, 7.0),
    doubleArrayOf(2.0, -1.0, 2.0 ,6.0),
    doubleArrayOf(1.0, -2.0, 1.0,0.0)
// GaussianElimination(A,usePivoting,epsilon) takes in 3 parameters:
//  A - a matrix
//  usePivoting - true if to use partial pivoting, e.g., for numerical stability.
//  epsilon - a precision parameter: when a number |x| ≤ ε, it is considered 0
val ops:GaussianElimination = GaussianElimination(a, true, 0.0)
//GaussianElimination has 4 methods, but we will only be looking at 1,
//which is U, as it returns to us a upper triangular matrix, 
//which is the matrix in row elchleon form
val ans:Matrix = ops.U()

This will return a different matrix that is also in Row Elchleon Form, which looks like


Now, we can begin deriving a solution from this simplified matrix. From the first matrix, we can already tell that

1.8z = 1.8

z = 1

Using the value z = 1, we can then derive the value of y using Row 2,




Lastly, we can use z = 1 and y = 2 to derive the value of x in Row 1,





Therefore, using a matrix in Row Elchleon Form, we have solved for the solution, where




We can use NM Dev class BackwardSubstitution to implement a backward substitution algorithm that will compute the answer in the same way as our manual working. It works on an upper triangular matrix, which in our case means a matrix in Row Elchleon Form, and starts first by computing X𝑛, then substitutes that backward into the next equation to solve for X𝑛−1, and repeats until X1. An example of the code is shown below.

					val tri:UpperTriangularMatrix = UpperTriangularMatrix(
    doubleArrayOf(2.0, 1.0, -1.0),
    doubleArrayOf(-2.5, 1.5),
)) // matrix of the variables in upper triangular shape
val b:Vector = ans.getColumn(4) //column for constants
val solver:BackwardSubstitution  = BackwardSubstitution()
output:[3.000000, 2.000000, 1.000000]

This gives us the same answer of




We are also able to further reduce the matrix into Reduced Row Elchleon Form. Starting from this matrix,


we will try to manipulate the leading entry to a 1 while ensuring that each column containing a leading 1 has zeros in all its other entries. Firstly, we can divide Row 3 by 1.8,

\begin{bmatrix}2&1&-1&7\\0&-2.5&1.5&-3.5\\0&0&1.8/1.8&1.8/1.8\end{bmatrix} =  \begin{bmatrix}2&1&-1&7\\0&-2.5&1.5&-3.5\\0&0&1&1\end{bmatrix}

then we can divide Row 2 by -2.5,

\begin{bmatrix}2&1&-1&7\\0&-2.5/-2.5&1.5/-2.5&-3.5/-2.5\\0&0&1&1\end{bmatrix} = \begin{bmatrix}2&1&-1&7\\0&1&-0.6&1.4\\0&0&1&1\end{bmatrix}

then we can subtract Row 1 with Row 2 to give us

\begin{bmatrix}2&1-1&-1+0.6&7-1.4\\0&1&-0.6&1.4\\0&0&1&1\end{bmatrix} = \begin{bmatrix}2&0&-0.4&5.6\\0&1&-0.6&1.4\\0&0&1&1\end{bmatrix}

From here, we can divide Row 1 by 2

\begin{bmatrix}2/2&0&-0.4/2&5.6/2\\0&1&-0.6&1.4\\0&0&1&1\end{bmatrix} = \begin{bmatrix}1&0&-0.2&2.8\\0&1&-0.6&1.4\\0&0&1&1\end{bmatrix}

Now, we will have to subtract Row 2 with -0.6*Row 3

\begin{bmatrix}1&0&-0.2&2.8\\0&1&-0.6+0.6&1.4+0.6\\0&0&1&1\end{bmatrix} = \begin{bmatrix}1&0&-0.2&2.8\\0&1&0&2\\0&0&1&1\end{bmatrix}

And add 0.2*Row 3 to Row 1

\begin{bmatrix}1&0&-0.2+0.2&2.8+0.2\\0&1&0&2\\0&0&1&1\end{bmatrix} = \begin{bmatrix}1&0&0&3\\0&1&0&2\\0&0&1&1\end{bmatrix}

From here, the Matrix is finally now in Reduced Row Elchleon Form, and we can easily deduce the answer from looking at the matrix




Though reducing the matrix into Reduced Row Elchleon Form takes a longer time it is easier to compute or read out the solution when the matrix is in such a form. We can use the NM Dev class GaussJordanElimination to simplfy our matrix into Reduced Row Elchleon Form.

					%use s2
//var a the matrix for the linear system
var a = DenseMatrix(arrayOf(
    doubleArrayOf(2.0, 1.0, -1.0, 7.0),
    doubleArrayOf(2.0, -1.0, 2.0 ,6.0),
    doubleArrayOf(1.0, -2.0, 1.0,0.0)
// GaussJordanElimination(A,usePivoting,epsilon) takes in 3 parameters:
//  A - a matrix
//  usePivoting - true if to use partial pivoting, e.g., for numerical stability.
//  epsilon - a precision parameter: when a number |x| ≤ ε, it is considered 0
val ops:GaussJordanElimination = GaussJordanElimination(a, true, 0.0)
val ans:Matrix = ops.U()
	[,1] [,2] [,3] [,4] 
[1,] 1.000000, 0.000000, 0.000000, 3.000000, 
[2,] 0.000000, 1.000000, 0.000000, 2.000000, 
[3,] 0.000000, 0.000000, 1.000000, 1.000000,