Tutorial 6

Interactive Notebooks

Click the badges above to run this tutorial interactively in your browser without installing anything!

• Binder: Free cloud-based Jupyter environment
• Colab: Google's free Jupyter notebook environment

Topics Covered

This tutorial covers the following topics using Python and SymPy:

• Matrix algebra: transpose, trace, determinant, inverse
• Orthogonal matrices
• Conditioning of a system and the determinant test
• Cramer's rule and the noise sensitivity of ill-conditioned systems
• Gauss elimination with partial pivoting (GEwPP)
• REF, RREF and rank

Introduction

import sympy as sy
from sympy.abc import x, y, z, t
from sympy import Matrix, Rational, eye, zeros, sqrt, N, simplify, nsimplify
sy.init_printing(use_latex=True)

def cramer(M, b):
    """Solve M x = b with Cramer's rule. Raises if the system is singular."""
    D = M.det()
    if D == 0:
        raise ZeroDivisionError(f"det(A) = {D}: Cramer's rule does not apply (singular system)")
    cols = []
    for i in range(M.cols):
        Mi = M.copy()
        Mi[:, i] = b            # replace column i with the right-hand side
        cols.append(Mi.det() / D)
    return Matrix(cols)

caution

In Python 1/4 is the float 0.25, so 1/4*(...) silently converts an exact SymPy result into a floating-point one. Write (...)/4, or use sy.Rational(1, 4).

Question 1

Solve for $\boldsymbol{F}$ in

$\frac{1}{4}\left[\frac{220}{\operatorname{trace}(C)}A^{T}B-C\right]=\frac{1}{3}EF-\frac{20}{\det(D)}D$

Rearranging, $F=\left(\frac{1}{3}E\right)^{-1}\left(\frac{1}{4}\left[\frac{220}{\operatorname{trace}(C)}A^{T}B-C\right]+\frac{20}{\det(D)}D\right)$.

A = Matrix([[1, 0, 2], [3, -1, -3]])
B = Matrix([[4, 6, 8], [-2, 5, 7]])
C = Matrix([[10, 11, 3], [-10, -27, -25], [70, -29, -27]])
D = Matrix([[1, 2, 0], [-1, 3, 0], [0, 0, 4]])
E = Matrix([[2, -2, 1], [1, 2, 2], [2, 1, -2]])

display(C.trace(), D.det())          # -44, 20
display(A.T * B)

First verify that $\frac{1}{3}E$ really is orthogonal, so that $\left(\frac{1}{3}E\right)^{-1}=\left(\frac{1}{3}E\right)^{T}$.

E3 = E/3
display(simplify(E3 * E3.T))                       # the identity matrix
display(simplify(E3.inv() - E3.T) == zeros(3, 3))  # True -> orthogonal

Now evaluate the bracket and solve for $F$.

inner = ((220/C.trace()) * (A.T * B) - C)/4 + (20/D.det()) * D
display(inner)

F = E3.T * inner                  # = (E/3)^-1 * inner
display(F)

$$F=\begin{bmatrix}-23 & -\frac{16}{3} & -\frac{25}{3}\\[4pt] -13 & \frac{97}{3} & \frac{121}{3}\\[4pt] 23 & -\frac{17}{3} & -\frac{41}{3}\end{bmatrix}$$

Question 2

For each case, comment on the condition of the coefficient matrix in terms of its determinant, without solving.

case1 = Matrix([[3, 7], [1, -2]])                      # 3x2 + 7x1 = -5 ; x2 - 2x1 = 1/2
case2 = Matrix([[1, 2], [2, Rational(3999, 1000)]])    # x2 + 2x1 = 4 ; 2x2 + 3.999x1 = 7.999
case3 = Matrix([[5, 20], [1, 4]])                      # 5x2 + 20x1 = -1 ; x2 + 4x1 = 2

display(case1.det())                 # -13        -> well-conditioned, unique solution
display(case2.det())                 # -1/1000    -> ill-conditioned, det is almost 0
display(case3.det())                 # 0          -> singular

Case	$\det(A)$	Condition	Geometry	Solution
1	$-13$	well-conditioned	one intersection point (very different slopes)	unique
2	$-1/1000$	ill-conditioned	nearly parallel lines (almost equal slopes)	exists but extremely noise-sensitive
3	$0$	singular	parallel lines (identical slopes)	none

Question 3

Continue Q2 with Cramer's rule.

display(cramer(case1, Matrix([-5, Rational(1, 2)])))                 # <-1/2, -1/2>
display(cramer(case2, Matrix([4, Rational(7999, 1000)])))            # <2, 1>

Case 1 gives $x_{1}=x_{2}=-0.5$ and Case 2 gives $x_{1}=2$, $x_{2}=1$; both check out when substituted back. Case 3 has $\det(A)=0$, so Cramer's rule fails:

try:
    cramer(case3, Matrix([-1, 2]))
except ZeroDivisionError as err:
    display(err)                     # det(A) = 0 -> Cramer's rule does not apply

Because the two equations are inconsistent (they describe parallel lines), the system has no solution.

Question 4

Repeat the Cramer calculation with the measured (noisy) right-hand side.

b1_noise = Matrix([Rational(-5001, 1000), Rational(4999, 10000)])
display([N(v, 4) for v in cramer(case1, b1_noise)])     # [-0.5002, -0.5001]

b2_noise = Matrix([Rational(4001, 1000), Rational(7998, 1000)])
display([N(v, 4) for v in cramer(case2, b2_noise)])     # [-3.999, 4.000]

Case 1 is well-conditioned, so the tiny perturbation is not amplified and the answer stays close to $(-0.5,-0.5)$. Case 2 is ill-conditioned: a change of only $0.001$ in the right-hand side moves the solution from $(2,1)$ to about $(-3.999,4.000)$. Small measurement noise is hugely magnified.

Question 5

Solve the following system by GEwPP (scaling, partial pivoting, forward elimination, backward substitution).

$y+z+2t=0,\quad x+2y+3z+t=5,\quad x+2y+z=4,\quad 2x+y+z+t=3$

G = Matrix([[0, 1, 1, 2],
            [1, 2, 3, 1],
            [1, 2, 1, 0],
            [2, 1, 1, 1]])
b = Matrix([0, 5, 4, 3])

display(G.LUsolve(b))                            # LU decomposition with partial pivoting
display(G.row_join(b).rref()[0])                 # augmented matrix -> RREF

$\begin{bmatrix}x\\y\\z\\t\end{bmatrix}=\begin{bmatrix}1\\1\\1\\-1\end{bmatrix}$

The LUsolve method performs partial-pivoting elimination (SymPy does not apply the tutorial's row-scaling step, but the answer is the same); the rref of the augmented matrix shows the same answer in one step.

Question 6

Obtain the REF, the RREF, the number of linearly independent vectors, and the rank.

M = Matrix([[1, 3, 4],
            [3, 9, 12],
            [2, 6, 9]])

display(M.echelon_form())        # REF
display(M.rref())                # RREF + pivot columns
display(M.rank())                # 2

i. REF: $\begin{bmatrix}1&3&4\\0&0&1\\0&0&0\end{bmatrix}$

ii. RREF: $\begin{bmatrix}1&3&0\\0&0&1\\0&0&0\end{bmatrix}$ — note this is not the same as the REF, because the entry above the second pivot must also be cleared.

iii. There are 2 linearly independent vectors (2 pivot columns).

iv. Rank = 2, so the matrix is rank deficient (it is $3\times3$ but not full rank).

MA = Matrix([[1, Rational(1, 2), Rational(1, 2), Rational(1, 2)],
             [0, Rational(3, 4), Rational(1, 4), Rational(-1, 4)],
             [0, 0, Rational(2, 3), Rational(1, 3)],
             [0, 0, 0, Rational(6, 7)]])
display(MA.rank())               # 4 -> full rank

Matrix A from the question table is upper triangular with no zero on the diagonal, so its rank is 4 and it is full rank, matching the printed solution.

note

The REF and the RREF differ here: the REF still has a non-zero entry above the second pivot, and the RREF clears it.

Question 7

An electronics company produces transistors $T$, resistors $R$ and chips $C$.

M = Matrix([[4, 3, 2],      # copper
            [1, 3, 1],      # zinc
            [2, 1, 3]])     # glass
b = Matrix([960, 510, 610])

weekly = M.LUsolve(b)
display(weekly)                  # <120, 100, 90> per week
display(weekly / 5)              # <24, 20, 18> per day

So the company produces $120$ transistors, $100$ resistors and $90$ chips per week, i.e. $24$, $20$ and $18$ per day on average.

Question 8

Gravel pits — use Cramer's rule.

M = Matrix([[Rational(52, 100), Rational(20, 100), Rational(25, 100)],
            [Rational(30, 100), Rational(50, 100), Rational(20, 100)],
            [Rational(18, 100), Rational(30, 100), Rational(55, 100)]])
b = Matrix([4800, 5800, 5700])

display(M.det())                                       # 43/500 = 0.086
display([N(v, 8) for v in cramer(M, b)])               # [4005.814, 7131.395, 5162.791]

$$\begin{bmatrix}P_{1}\\P_{2}\\P_{3}\end{bmatrix} =\begin{bmatrix}4005.8\\7131.4\\5162.8\end{bmatrix}\text{m}^{3}$$

note

All percentages enter as decimals: $0.30$ fine gravel and $0.55$ coarse gravel in the first numerator, so the $30$ and $55$ of the raw table are not used directly.

Question 9

Three reactors linked by pipes. The mass-balance equations are $120c_{1}-20c_{2}=400$, $-80c_{1}+80c_{2}=0$ and $-40c_{1}-60c_{2}+120c_{3}=200$.

M = Matrix([[120, -20, 0],
            [-80, 80, 0],
            [-40, -60, 120]])
b = Matrix([400, 0, 200])

display(M.LUsolve(b))            # naive GE (no pivoting needed here)

$$\begin{bmatrix}c_{1}\\c_{2}\\c_{3}\end{bmatrix} =\begin{bmatrix}4\\4\\5\end{bmatrix}\text{mg/m}^{3}$$

Question 10

Three planes: $3x+6y-3z=-6$, $4x-2y+6z=2$, $3x+2y+z=-2$.

M = Matrix([[3, 6, -3],
            [4, -2, 6],
            [3, 2, 1]])
b = Matrix([-6, 2, -2])

display(M.det())                                        # 0 -> singular
display(M.rank(), M.row_join(b).rank())                 # 2, 2 -> consistent, infinite solutions
display(sy.linsolve((M, b), (x, y, z)))                 # {(-z, z - 1, z)}

Since $\det(M)=0$ but the rank of the coefficient matrix equals the rank of the augmented matrix, the system is consistent with infinitely many solutions. Setting $z=t$:

$x=-t,\qquad y=-1+t,\qquad z=t,\qquad t\in\mathbb{R}$

which is the line along which all three planes intersect. Note that matrix inversion and Cramer's rule are both unusable here — GEwPP (or linsolve) is the right tool.

---

Extra Learning Resources

Key Concepts to Master

1. Determinant test: $\det(A)\neq0$ → unique solution; $\det(A)=0$ → no solution or infinitely many
2. Well- vs ill-conditioned: a determinant close to zero means small input errors produce large output errors
3. Cramer's rule: $x_{i}=\det(A_{i})/\det(A)$, where $A_{i}$ has column $i$ replaced by $b$
4. Rank: the number of pivots; rank deficient means the rank is less than the matrix size
5. GEwPP: scaling to compare rows fairly, partial pivoting to avoid dividing by a tiny pivot, then forward elimination and back substitution

SymPy Resources

• SymPy Matrices — det, inv, rank, rref, echelon_form, LUsolve, row_join
• SymPy linsolve — solve singular and underdetermined systems
• SymPy Matrix expressions — symbolic matrix algebra

Common Pitfalls

• Python integer division: use /4 or sy.Rational(1, 4), never 1/4, when you need an exact result
• rref() returns a tuple (matrix, pivot_columns) — index [0] to get the matrix
• rref() and echelon_form() are not the same thing
• A matrix with a zero row can still be consistent — compare the ranks, don't just look at the determinant
• Floats such as 3.999 are exact decimal literals in SymPy only if you use Rational(3999, 1000); 3.999 alone becomes a binary float