41.2: Solve Linear Systems

Matrix Algebra with Computational Applications · LibreTexts

0%

Contents

    Text

    Size
    100%
    Theme
    Typeface
    Line length
    Spacing

    From LibreTexts, under CC BY-NC 4.0, as the source states it; its equations are written in TeX, which the PDF shows typeset. This page at LibreTexts.

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

    This page titled 41.2: Solve Linear Systems was last modified on Sat, 17 Sep 2022 15:14:42 GMT and is shared under a CC BY-NC 4.0 license and was authored, remixed, and/or curated by Dirk Colbry via source content that was edited to the style and standards of the LibreTexts platform.