12.1 — Using SciPy, System of Linear Equations

A system of linear equations is a collection of one or more linear equations:

\small \begin{aligned} a_{0,0}x_0 + a_{0,1}x_2 + \cdots + a_{0,n}x_n & = b_0 \\ a_{1,0}x_0 + a_{1,1}x_2 + \cdots + a_{1,n}x_n & = b_1 \\ & \vdots \\ a_{m,0}x_0 + a_{m,1}x_2 + \cdots + a_{m,n}x_n & = b_m \\ \end{aligned}

In matrix notation, the linear system is given by $\mathbf{Ax = b}$ , where

\small A = \begin{bmatrix} a_{0,0} & a_{0,1} & \cdots & a_{0,n} \\ a_{1,0} & a_{1,1} & \cdots & a_{1,n} \\ \vdots & & & \vdots \\ a_{m,0} & a_{m,1} & \cdots & a_{m,n} \\ \end{bmatrix} \ \ , \ \ \mathbf{x} = \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_n \end{bmatrix} \ \ , \ \ \mathbf{b} = \begin{bmatrix} b_0 \\ b_1 \\ \vdots \\ b_m \end{bmatrix}

Gauss-Jordan elimination to solve linear system

High-level overview of Gauss-Jordan elimination:

Given the matrices $\mathbf{A}$ and $\mathbf{b}$ of the linear system $\mathbf{Ax = b}$ we obtain augment matrix $\mathbf{M = [A | b]}$ by stacking $\mathbf{b}$ as the last column to the matrix $\mathbf{A}$ . $\small M = \left[ \begin{array}{cccc|c} a_{0,0} & a_{0,1} & \cdots & a_{0,n} & b_0 \\ a_{1,0} & a_{1,1} & \cdots & a_{1,n} & b_1 \\ \vdots & & & \vdots & \vdots \\ a_{m,0} & a_{m,1} & \cdots & a_{m,n} & b_m \\ \end{array} \right]$

Using elementary row operations, convert the augmented matrix $M$ into reduced row-echelon form (RREF).
- In this form, the solution is simply the last column of the matrix $M$ (if there is a unique solution). For example,
$\small \left[ \begin{array}{ccc|c} 0 & 8 & 4 & 2 \\ 2 & 5 & 1 & 5 \\ 4 & 10 & -1 & 1 \\ \end{array} \right]$
$\implies$
$\small \left[ \begin{array}{ccc|c} 1 & 0 & 0 & 4.125 \\ 0 & 1 & 0 & -1.25 \\ 0 & 0 & 1 & 3. \\ \end{array} \right]$

Reduced row-echelon form

A matrix is in reduced row-echelon form (RREF) if it satisfies all of the following conditions.

In every row the leftmost non-zero entry is $1$ (called the leading $1$ ).
If a column contains a leading $1$ (i.e. first non-zero entry in the column is $1$ ), then all other entries in that column is $0$ .
The leading $1$ of any given row is always to the right of the leading $1$ of the row above it.
Any rows consisting entirely of zeroes are placed at the bottom of the matrix.

Examples:

\small A = \begin{bmatrix} 1 & 0 & 1 & 0\\ 0 & 1 & 0 & 2 \end{bmatrix} \quad , \quad B = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}

\small C = \begin{bmatrix} 1 &0 \\ 0 & 1\\ 0 & 0\\ 0 & 0\\ 0 & 0 \end{bmatrix} \quad , \quad D = \begin{bmatrix} 1 & 0 & 0 & 4 & 0 & 0 \\ 0 & 1 & 0 & 3 & 0 & 0\\ 0 & 0 & 1 & 3 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}

Elementary row operations

We use the following operations to convert a matrix M into RREF:

Row Swap. Exchange any two rows: $M[i] \leftrightarrow M[j]$
Scalar Multiplication. Multiply any row by a non-zero constant:
$M[i] \gets k \cdot M[i] , k \neq 0$
Row Addition. Add multiple of one row to another row:
$M[i] \gets M[i] + k \cdot M[j] , i \neq j$

Download file elementary_row_ops.py from Ed.

Gauss-Jordan Elimination algorithm using NumPy

Check this link to understand how the algorithm works.

Check the program gaussjordan.py for an implementation of the algorithm. (Just for understanding; not for exam).

Based roughly on the pseudocode for ToReducedRowEchelonForm in https://rosettacode.org/wiki/Reduced_row_echelon_form

Using scipy to solve linear systems

1import numpy as np
2import scipy.linalg as la
3
4A = np.array([[0, 8, 4], [2, 5, 1], [4, 10, -1]], dtype=float)
5b = np.array([2, 5, 1], dtype=float)
6
7result = la.solve(A, b)
8print(result)
9# [ 4.125 -1.25   3.   ]