41.2: Solve Linear Systems

Login with LibreOne to run this code cell interactively.

If you have already signed in, please refresh the page.

Login

import numpy as np
import sympy as sym

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$$

Login with LibreOne to run this code cell interactively.

If you have already signed in, please refresh the page.

Login

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))

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.

Login with LibreOne to run this code cell interactively.

If you have already signed in, please refresh the page.

Login

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

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.

Login with LibreOne to run this code cell interactively.

If you have already signed in, please refresh the page.

Login

## 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))

## 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))