Solving Linear Systems (SymPy)
Most linear algebra courses start by considering how to solve a system of linear equations.
In matrix notation, this linear system of \(n\) equations in \(n\) unknowns (\(x_0, \ldots, x_{n-1}\)), can be written as
where
We'll look at how to use SymPy to solve linear systems, including those with free variables. We'll also look at how to do row operations (Gaussian elimination).
First we import the necessary libraries.
.solve
and .linsolve
Let's start with a 3x3 linear system:
Verify (by substitution) that \(\vec{x} = [0,2,-1]\) is a solution.
System of Equations
To solve a linear system represented by a system of linear equations we use the sympy.solve
funtion.
- \(A\) may can be a rectangular matrix.
Matrix Form
For a linear system expressed in matrix form \(A\vec{x} = \vec{b}\) we can use linsolve
.
Augmented Matrix Form
For a linear system expressed as an augmented matrix \([A|\vec{b}]\) we can use linsolve
.
M = Matrix([[1,1,-2],[0,2,1],[2,1,-5]])
b = Matrix([4,3,7])
M_aug = M.row_join(b)
linsolve(M_aug, (x, y, z))
System with Free Variables
The three methods above can be used to solve a system with free variables.
Consider the system
which has solution space \(\{ (-t-1, t, 2) : t \in \mathbb{R} \}\).
Example
Consider another example: $$ 2x_1 + 3x_2 = 11 $$ $$ 3x_1 + 6x_2 = 7 $$
\(\left\{ \left( 15, -\frac{19}{3} \right) \right\}\)
We can verify this answer is correct by observing \(A\vec{x} = \vec{b}\).
pprint
is pretty print, which displays matrices in the proper format.
\(\left[ \begin{matrix} 11 \\ 7 \\ \end{matrix} \right]\)
Gaussian Elimination
The general procedure learned to solve a system of linear equations is Gaussian elimination. The goal is to apply row operations to "eliminate" all the variables except for one in each equation (if possible), and this resulting linear system has the same set of solutions as the original, but is easier to solve. We generally do this on the matrix form of the system, and call this the reduced row echelon form (rref) of the matrix. The elementary row operations we use are:
Elementary Row Operations
- Replace row \(j\) with \(k\) times row \(i\) added to row \(j\). \(\quad R_j \leftarrow kR_i +R_j\)
- Multiply row \(i\) by a scalar \(k\). \(\quad R_i \leftarrow kR_i\)
- Swap rows \(i\) and \(j\). \(\quad R_i \leftrightarrow R_j\)
For example, consider the matrix
Since we are working in python we will index rows starting at \(0\). In other words, an \(n\times m\) matrix has rows \(R_0\) to \(R_{n-1}\).
Here are the row operations taking it to rref:
To do row operations in SymPy us the elemetary_row_op
method on a Matrix
object.
elementary_row_op(op='n->kn', row=None, k=None, row1=None, row2=None)
Performs the elementary row operation
"n->kn"
(row n goes to k*n)"n<->m"
(swap row n and row m)"n->n+km
" (row n goes to row n + k*row m)row
is the index of the row to be replaced.row1
is the row that will be added torow
row1
androw2
are the two rows in the swap.
Example
We'll redo the example above, this time using SymPy.
A = Matrix([[2,1,0],[1,1,-1],[0,-1,2]])
pprint(A)
print("Swap rows 0 and 1:")
A1 = A.elementary_row_op("n<->m",0,1)
pprint(A1)
print("Add -2 times row 0 to row 1:")
A2 = A1.elementary_row_op("n->n+km",row=1,k=-2,row1=0)
pprint(A2)
print("Add -1 times row 1 to row 2:")
A3 = A2.elementary_row_op("n->n+km",row=2,k=-1,row1=1)
pprint(A3)
print("Add 1 times row 1 to row 0:")
A4 = A3.elementary_row_op("n->n+km",row=0,k=1,row1=1)
pprint(A4)
print("Multiply row 1 by -1:")
A5 = A4.elementary_row_op("n->kn",row=1,k=-1)
pprint(A5)
[ 2 1 0]
[ 1 1 -1]
[ 0 -1 2]
Swap rows 0 and 1:
[ 1 1 -1 ]
[ 2 1 0 ]
[ 0 -1 2 ]
Add -2 times row 0 to row 1:
[ 1 1 -1 ]
[ 0 -1 2 ]
[ 0 -1 2 ]
Add -1 times row 1 to row 2:
[ 1 1 -1 ]
[ 0 -1 2 ]
[ 0 0 0 ]
Add 1 times row 1 to row 0:
[ 1 0 1 ]
[ 0 -1 2 ]
[ 0 0 0 ]
Multiply row 1 by -1:
[ 1 0 1 ]
[ 0 1 -2 ]
[ 0 0 0 ]
Solving a linear system with Gaussian Elimination
We'll consider the following linear system which has more unknowns than equations. This means the coefficient matrix will not be square.
Write it as an augmented matrix:
Use row operations to take it to rref.
A = Matrix([[2,0,4,4],[-1,0,-1,0],[-3,0,0,6]])
b = Matrix([14,-3,3])
aug_mat = A.row_join(b) # augmented matrix
aug_mat
\(\left[ \begin{matrix} 2 & 0 & 4 & 4 & 14 \\ -1 & 0 & -1 & 0 & -3 \\ -3 & 0 & 0 & 6 & 3 \\ \end{matrix} \right]\)
\(\left[ \begin{matrix} -1 & 0 & -1 & 0 & -3 \\ 2 & 0 & 4 & 4 & 14 \\ -3 & 0 & 0 & 6 & 3 \\ \end{matrix} \right]\)
A2 = A1.elementary_row_op("n->n+km",row=1,k=2,row1=0)
A3 = A2.elementary_row_op("n->n+km",row=2,k=-3,row1=0)
A4 = A3.elementary_row_op("n->kn",row=0,k=-1)
A4
\(\left[ \begin{matrix} 1 & 0 & 1 & 0 & 3 \\ 0 & 0 & 2 & 4 & 8 \\ 0 & 0 & 3 & 6 & 12 \\ \end{matrix} \right]\)
- Use
k=Rational(1,3)
to avoid using floating point numbers.
\(\left[ \begin{matrix} 1 & 0 & 1 & 0 & 3 \\ 0 & 0 & 2 & 4 & 8 \\ 0 & 0 & 1 & 2 & 4 \\ \end{matrix} \right]\)
A6 = A5.elementary_row_op("n->n+km",row=1,k=-2,row1=2)
A7 = A6.elementary_row_op("n->n+km",row=0,k=-1,row1=2)
A7
\(\left[ \begin{matrix} 1 & 0 & 0 & -2 & -1 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 2 & 4 \\ \end{matrix} \right]\)
\(\left[ \begin{matrix} 1 & 0 & 0 & -2 & -1 \\ 0 & 0 & 1 & 2 & 4 \\ 0 & 0 & 0 & 0 & 0 \\ \end{matrix} \right]\)
At this point we have reduced the matrix to rref. We see there are two free variables (\(x_2\) and \(x_4\)). The general solution is:
.rref
SymPy has a built-in function to put a matrix in reduced row echelon form (rref). To put a matrix into reduced row echelon form use the method .rref
. .rref
returns a tuple of two elements. The first is the reduced row echelon form, and the second is a tuple of indices of the pivot columns. The length of this tuple is the rank of the matrix.
Example
Compute the rref of the matrix
\(\left(\left[ \begin{matrix} 1 & 0 & 1 \\ 0 & 1 & -2 \\ 0 & 0 & 0 \\ \end{matrix} \right], (0,1) \right)\)
Example
Compute the rref of the matrix
\(\left(\left[ \begin{matrix} 1 & 0 & 0 & -2 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 0 & 0 \\ \end{matrix} \right], (0,2) \right)\)
Example
For the linear system, find the rref of the augmented matrix.
A = Matrix([[2,0,4,4],[-1,0,-1,0],[-3,0,0,6]])
b = Matrix([14,-3,3])
aug_mat = A.row_join(b) # augmented matrix
pprint(aug_mat.rref())
\(\left(\left[ \begin{matrix} 1 & 0 & 0 & -2 & -1\\ 0 & 0 & 1 & 2 & 4\\ 0 & 0 & 0 & 0 & 0\\ \end{matrix} \right], (0,2) \right)\)
Computing an Inverse with Gaussian Elimination
Recall, if we have an invertible matrix \(A\), then we can compute the inverse using Gaussian elimination as follows:
Matrix \(B\) is the inverse of \(A\).
Note
If you are wondering why the inverse matrix can be computed in this way, it comes back to elementary matrices. Suppose \(n\) row operations are applied to take \(A\) to the identity \(I\). Let \(E_1, E_2, \ldots E_n\) be the \(n\) elementary matrices corresponding to these row operations. Then, since row operations could be done by multiplying by elementary matrices, it follows that
Therefore, multiplying \(A^{-1}\) to the right or each side:
This means \(A^{-1}\) is obtained by doing the row operations on \(I\).
Example
Let's find the inverse of
\(\left[ \begin{matrix} 1 & 2 & 1 & 0\\ 3 & 4 & 0 & 1\\ \end{matrix} \right]\)
\(\left[ \begin{matrix} 1 & 2 & 1 & 0\\ 0 & -2 & -3 & 1\\ \end{matrix} \right]\)
A2 = A1.elementary_row_op("n->n+km",row=0,k=1,row1=1)
A3 = A2.elementary_row_op("n->kn",row=1,k=-1/2)
A3
\(\left[ \begin{matrix} 1 & 0 & -2 & 1\\ 0 & 1 & 1.5 & -0.5\\ \end{matrix} \right]\)
Therefore, we have
We can verify we found the inverse by checking that \(A A^{-1} = I\).
\(\left[ \begin{matrix} 1.0 & 0 \\ 0 & 1.0 \\ \end{matrix} \right]\)
Nullspace
The nullspace of a matrix can be computed with .nullspace
. This returns an orthogonal basis for the nullspace.
\(\left[ \left[ \begin{matrix} 0 \\ 1 \\ 0 \\ 0 \end{matrix} \right], \left[ \begin{matrix} -1 \\ 0 \\ 1 \\ 1 \end{matrix} \right] \right]\)