41.2: Solve Linear Systems

Last updated
Save as PDF

Page ID: 70549

\( \newcommand{\vecs}[1]{\overset { \scriptstyle \rightharpoonup} {\mathbf{#1}} } \) \( \newcommand{\vecd}[1]{\overset{-\!-\!\rightharpoonup}{\vphantom{a}\smash {#1}}} \)\(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\) \(\newcommand{\id}{\mathrm{id}}\) \( \newcommand{\Span}{\mathrm{span}}\) \( \newcommand{\kernel}{\mathrm{null}\,}\) \( \newcommand{\range}{\mathrm{range}\,}\) \( \newcommand{\RealPart}{\mathrm{Re}}\) \( \newcommand{\ImaginaryPart}{\mathrm{Im}}\) \( \newcommand{\Argument}{\mathrm{Arg}}\) \( \newcommand{\norm}[1]{\| #1 \|}\) \( \newcommand{\inner}[2]{\langle #1, #2 \rangle}\) \( \newcommand{\Span}{\mathrm{span}}\)\(\newcommand{\AA}{\unicode[.8,0]{x212B}}\)

import numpy as np
import sympy as sym

When we have the same number of equations as unknowns, we have the following system \(Ax=b\) with a square matrix \(A\). In this section, we assume that the matrix \(A\) has full rank. That is the system has an unique solution. We talked about many ways to solve this system of equations. Some examples are:

Jacobian iteration/ Gauss-Seidel iteration
\(x=A^{−1}b\)
Gaussian elimination
LU decomposition

In this assignment, we will show that we can solve it by QR decomposion.

Consider the following system of equations:

\[\begin{split}\begin{bmatrix}5&-2&2 \\ 4 & -3 &4 \\ 4& -6 &7 \end{bmatrix}\begin{bmatrix}x_1\\x_2\\x_3\end{bmatrix}=\begin{bmatrix}1\\2\\3\end{bmatrix}\end{split} \nonumber \]

A = np.matrix([[5, -2, 2], [4, -3, 4], [4,-6,7]])
b = np.matrix([[1],[2],[3]])
display(sym.Matrix(A))
display(sym.Matrix(b))

Back substitution

Let’s first implement the back substitution in Python. The back substitution function back_subst solves the system \(Rx=b\) where \(R\) is an upper-triangular matrix.

When we solve for \(x\), we start with \(x_n\): \(x_n = {b_n\over R_{n,n}}\). Then we solve for \(x_{n-1}\) as \(x_{n-1} = {b_{n-1}-R_{n-1,n}x_n\over R_{n-1,n-1}}\). Then we can find \(x_{n-2},\cdots,x_1\): \(x_{n-2} = {b_{n-2}-R_{n-2,n-1}x_{n-1}-R_{n-2,n}x_n\over R_{n-2,n-2}}\). We can solve for all components of \(x\) in the reversed order. So this is called back substitution.

Do This

Complete the following code for back substitution.

def back_subst(R,b):
    n = R.shape[0]; x = np.zeros(n);
    for i in reversed(range(n)):
        x[i] = b[i]
        for j in range(i+1,n):            
## Your code starts ##            
            x[i] =   # Complet this line of code.
## Your code ends ##            
        x[i] = x[i]/R[i,i]
    return np.matrix(x).T

Do This

When we solve for \(Ax=b\) with QR decomposition. We have the following steps:

Find the QR decomposition of \(A\)
From \(QRx=b\), we obtain \(Rx=Q^{\top}b\)
Solve for \(x\) using back substitution.

Use these steps to solve \(Ax=b\) with the given \(A\) and \(b\). Compare the result with np.linalg.solve.

## Your code starts ##
x =  
## Your code ends ##
print(type(x))   # x is a column vector in the np.matrix type
np.allclose(x, np.linalg.solve(A,b))