Numerical calculation method chapter5 Direct method for solving linear equations

• Numerical calculation method Chapter5. A direct method for solving linear equations
  • 0. Problem description
  • 1. Elimination method
    • 1. Trigonometric equations
      • 1. Diagonal equations
      • 2. Lower trigonometric equations
      • 3. Upper triangular equations
    • 2. Gauss Elimination method
    • 3. Gauss-Jordan Elimination method
  • 2. Direct decomposition method
    • 1. Dolittle decompose
    • 2. Courant decompose
    • 3. Catch up method
    • 4. Symmetric positive definite matrix $LD{L}^T$ decompose

0. Problem description

This chapter examines how to solve linear equations ：

$$\{\begin{array}{rl}{a}_{11}{x}_1+a_{12}{x}_2+...+a_{1n}{x}_n& =b_1\\ a_{21}{x}_1+a_{22}{x}_2+...+a_{2n}{x}_n& =b_2\\ ...& \\ a_{n1}{x}_1+a_{n2}{x}_2+...+a_{nn}{x}_n& =b_n\end{array}$$

Or it can be expressed as a matrix ：

$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1n}\\ a_{21}& a_{22}& ...& a_{2n}\\ ...& & & \\ a_{n1}& a_{n2}& ...& a_{nn}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_n\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_n\end{array}\right)$$

1. Elimination method

1. Trigonometric equations

First , Let's look at some special forms of equations ：

1. Diagonal equations

The functional form of diagonal equations is as follows ：

$$\left(\begin{array}{cccc}{a}_{11}& & & \\ & a_{22}& & \\ & & ...& \\ & & & a_{nn}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_n\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_n\end{array}\right)$$

obviously , You can write directly $x_i=b_i/a_{ii}$, among ,$i\in [1,n]$.

2. Lower trigonometric equations

below , Let's look at a slightly more complicated situation , That is, the case of lower triangular matrix ：

$$\left(\begin{array}{cccc}{a}_{11}& & & \\ a_{21}& a_{22}& & \\ & & ...& \\ a_{n1}& a_{n2}& ...& a_{nn}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_n\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_n\end{array}\right)$$

here , The difficulty of the problem will be a little more complicated , But it is also possible to give a direct answer as follows ：

$$x_i=\frac{b_i-{\sum }_{j<i}{a}_{ij}{x}_j}{a_{i}i}$$

We can give python The pseudocode is as follows ：

def down_triangle(A, b):
    n = len(A)
    x = [0 for _ in range(n)]
    for i in range(n):
        y = b[i]
        for j in range(i):
            y -= A[i][j] * x[j]
        x[i] = y / A[i][i]
    return x

3. Upper triangular equations

alike , For the case of upper trigonometric function , We also have ：

$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1n}\\ & a_{22}& ...& a_{2n}\\ & & ...& \\ & & & a_{nn}\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_n\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ ...\\ b_n\end{array}\right)$$

It can be solved that ：

$$x_i=\frac{b_i-{\sum }_{j>i}{a}_{ij}{x}_j}{a_{i}i}$$

alike , We can give python The pseudocode is as follows ：

def up_triangle(A, b):
    n = len(A)
    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        y = b[i]
        for j in range(i+1, n):
            y -= A[i][j] * x[j]
        x[i] = y / A[i][i]
    return x

2. Gauss Elimination method

Now? , Let's look at the solution of a general form of multivariate linear equation .

The core idea is to convert it into a triangular matrix and then solve it .

We also give python The pseudocode is as follows ：

def gauss_elimination(A, b):
    n = len(A)
    for i in range(n-1):
        for k in range(i+1, n):
            c = A[k][i] / A[i][i]
            for j in range(i, n):
                A[k][j] -= c * A[i][j]
            b[k] -= c * b[i]
    return up_triangle(A, b)

3. Gauss-Jordan Elimination method

Gauss-Jordan Elimination method and the above Gauss The elimination method is essentially the same , however Gauss The elimination method is to convert a general matrix into a triangular matrix , and Gauss-Jordan The elimination method is to convert a general matrix into a diagonal matrix .

alike , We give python The pseudocode is as follows ：

def gauss_jordan_elimination(A, b):
    n = len(A)
    for i in range(n-1):
        for k in range(i+1, n):
            c = A[k][i] / A[i][i]
            for j in range(i, n):
                A[k][j] -= c * A[i][j]
            b[k] -= c * b[i]
    
    for j in range(n-1, 0, -1):
        for i in range(j):
            c = A[i][j] / A[j][j]
            A[i][j] = 0
            b[i] -= c * b[j]

    res = [b[i] / A[i][i] for i in range(n)]
    return res

2. Direct decomposition method

There is no essential difference between the direct decomposition method and the above elimination method , But the difference is , The core idea of the direct decomposition method is to constantly try to convert the general matrix into the form of triangular matrix and then solve it based on the specificity of triangular matrix .

1. Dolittle decompose

Dolittle The idea of decomposition is to say that the general $n$ Order matrix $A$ convert to $A=LU$, among ,$L$ Is a diagonal element for 1 The lower triangular matrix of , and $U$ It's an upper triangular matrix .

namely ：

$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1n}\\ a_{21}& a_{22}& ...& a_{2n}\\ & & ...& \\ a_{n1}& a_{n2}& ...& a_{nn}\end{array}\right)=\left(\begin{array}{cccc}1& & & \\ l_{21}& 1& & \\ & ...& & \\ l_{n1}& l_{n2}& ...& 1\end{array}\right)\left(\begin{array}{cccc}{u}_{11}& u_{12}& ...& u_{1n}\\ & u_{22}& ...& u_{2n}\\ & & ...& \\ & & & u_{nn}\end{array}\right)$$

And two trigonometric matrices $L$ and $U$ Are relatively easy to solve .

We are easy to have ：

$$\{\begin{array}{rl}{u}_{ij}& =a_{ij}-\sum _{k=1}^{i-1}{l}_{ik}{u}_{kj}\\ l_{ij}& =(a_{ij}-\sum _{k=1}^{j-1}{l}_{ik}{u}_{kj})/u_{jj}\end{array}$$

here , For the equation to be solved $Ax=b$ Can become ：

$$LUx=b$$

further , We can split it into two equations ：

$$\{\begin{array}{rl}Ly& =b\\ Ux& =y\end{array}$$

Solve the above two equations respectively , The final solution can be obtained $x$：

$$\{\begin{array}{rl}{y}_i& =b_i-\sum _{j=1}^{i-1}{l}_{ij}{y}_j\\ x_i& =(y_i-\sum _{j=i+1}^n{u}_{ij}{x}_j)/u_{ii}\end{array}$$

alike , We can give python The pseudocode is as follows ：

def dolittle_solve(A, b):
    n = len(A)
    L = [[0 for _ in range(n)] for _ in range(n)]
    U = [[0 for _ in range(n)] for _ in range(n)]

    for k in range(n):
        L[k][k] = 1
        for j in range(k, n):
            U[k][j] = A[k][j]
            for r in range(k):
                U[k][j] -= L[k][r] * U[r][j]
        for i in range(k+1, n):
            L[i][k] = A[i][k]
            for r in range(k):
                L[i][k] -= L[i][r] * U[r][k]
            L[i][k] /= U[k][k]
    
    y = [0 for _ in range(n)]
    for i in range(n):
        y[i] = b[i]
        for j in range(i):
            y[i] -= L[i][j] * y[j]

    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = y[i]
        for j in range(i+1, n):
            x[i] -= U[i][j] * x[j]
        x[i] /= U[i][i]

    return x

2. Courant decompose

Courant Break down and Dolittle Decomposition is essentially indistinguishable , It is still the general parameter matrix $A$ Split into two three matrices , But for the Courant In terms of decomposition ,$L$ Is a general lower triangular matrix , and $U$ Is that the diagonal elements are all 1 The upper triangular matrix of .

namely ：

$$\left(\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1n}\\ a_{21}& a_{22}& ...& a_{2n}\\ & & ...& \\ a_{n1}& a_{n2}& ...& a_{nn}\end{array}\right)=\left(\begin{array}{cccc}{l}_{11}& & & \\ l_{21}& l_{22}& & \\ & ...& & \\ l_{n1}& l_{n2}& ...& l_{nn}\end{array}\right)\left(\begin{array}{cccc}1& u_{12}& ...& u_{1n}\\ & 1& ...& u_{2n}\\ & & ...& \\ & & & 1\end{array}\right)$$

here , In imitation, we will soon get $L$ and $U$ The expression of is as follows ：

$$\{\begin{array}{rl}{l}_{ij}& =a_{ik}-\sum _{r=1}^{k-1}{l}_{ir}{u}_{rk}\\ u_{ij}& =(a_{kj}-\sum _{r=1}^{k-1}{l}_{kr}{u}_{rj})/l_{kk}\end{array}$$

alike , We can determine the final solution by the following function $x$ To solve the ：

$$\{\begin{array}{rl}Ly& =b\\ Ux& =y\end{array}$$

The final result expression is as follows ：

$$\{\begin{array}{rl}{y}_i& =(b_i-\sum _{j=1}^{i-1}{l}_{ij}{y}_j)/l_{ii}\\ x_i& =y_i-\sum _{j=i+1}^n{u}_{ij}{x}_j\end{array}$$

Also given python The pseudo code is implemented as follows ：

def courant_solve(A, b):
    n = len(A)
    L = [[0 for _ in range(n)] for _ in range(n)]
    U = [[0 for _ in range(n)] for _ in range(n)]

    for k in range(n):
        for i in range(k, n):
            L[i][k] = A[i][k]
            for r in range(k):
                L[i][k] -= L[i][r] * U[r][k]
        U[k][k] = 1
        for j in range(k+1, n):
            U[k][j] = A[k][j]
            for r in range(k):
                U[k][j] -= L[k][r] * U[r][j]
            U[k][j] /= L[k][k]
        
    
    y = [0 for _ in range(n)]
    for i in range(n):
        y[i] = b[i]
        for j in range(i):
            y[i] -= L[i][j] * y[j]
        y[i] /= L[i][i]

    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = y[i]
        for j in range(i+1, n):
            x[i] -= U[i][j] * x[j]

    return x 

3. Catch up method

The catch-up method is essentially just Courant A special case of decomposition .

Its solution and Courant Decomposition is completely consistent , The only difference is , there $A$ Matrix is a special ternary matrix , namely ：

$$A=\left(\begin{array}{ccccc}{a}_1& b_1& & & \\ c_2& a_2& b_2& & \\ & & ...& & \\ & & c_{n-1}& a_{n-1}& b_{n-1}\\ & & & c_n& a_n\end{array}\right)$$

thus , Expand the matrix $L$ and $U$ It can also be synchronously modified to ：

$$L=\left(\begin{array}{cccc}{u}_1& & & \\ w_2& u_2& & \\ & ...& ...& \\ & & w_n& u_n\end{array}\right),U=\left(\begin{array}{cccc}1& v_1& & \\ & 1& v_2& \\ & & ...& \\ & & & 1\end{array}\right)$$

here , Because of the matrix $L$ and $U$ The particularity of , We can directly write the final solution as follows ：

$$\{\begin{array}{rl}{u}_i& =a_i-c_i{v}_{i-1}\\ v_i& =b_i/u_i\\ y_i& =(f_i-c_i{y}_{i-1})/u_i\\ x_i& =(y_i-v_i{x}_{i+1})\end{array}$$

4. Symmetric positive definite matrix $LD{L}^T$ decompose

alike ,$LD{L}^T$ Decomposition and the above Dolittle Break down and Courant Decomposition is essentially the same , But here's the difference , The problem here is Symmetric positive definite matrices .

therefore , For symmetric positive definite matrices , We can always have ：

$$\begin{array}{rl}A& =LD{L}^T\\ & =\left(\begin{array}{cccc}1& & & \\ l_{21}& 1& & \\ & ...& & \\ l_{n1}& l_{n2}& ...& 1\end{array}\right)\left(\begin{array}{cccc}{d}_1& & & \\ & d_2& & \\ & & ...& \\ & & & d_n\end{array}\right)\left(\begin{array}{cccc}1& l_{12}& ...& l_{1n}\\ & 1& ...& l_{2n}\\ & & ...& \\ & & & 1\end{array}\right)\end{array}$$

here , We can solve ：

$$\{\begin{array}{rl}{d}_k& =a_{kk}-\sum _{r=1}^{k-1}{d}_r{l}_{kr}^2\\ l_{ik}& =(a_{ik}-\sum _{r=1}^{k-1}{d}_r{l}_{ir}{l}_{kr})/d_k\end{array}$$

, in turn, , We can put the equation $Ax=LD{L}^{T}x=b$ Disassemble into three steps to complete ：

$$\{\begin{array}{rl}Lz& =b\\ Dy& =z\\ L^{T}x& =y\end{array}$$

Solution ：

$$\{\begin{array}{rl}{z}_i& =b_i-\sum _{j=1}^{i-1}{l}_{ij}{z}_j\\ y_i& =z_i/d_i\\ x_i& =y_i-\sum _{j=i+1}^n{l}_{ji}{x}_j\end{array}$$

Again , We give python The pseudocode is as follows ：

def ldl_solve(A, b):
    
    n = len(A)
    l = [[0 for _ in range(n)] for _ in range(n)]
    d = [0 for _ in range(n)]
    for k in range(n):
        d[k] = A[k][k]
        for r in range(k):
            d[k] -= d[r] * l[k][r] * l[k][r]
        
        l[k][k] = 1
        for i in range(k+1, n):
            l[i][k] = A[i][k]
            for r in range(k):
                l[i][k] -= d[r] * l[i][r] * l[k][r]
            l[i][k] /= d[k]

    z = [0 for _ in range(n)]
    for i in range(n):
        z[i] = b[i]
        for j in range(i):
            z[i] -= l[i][j] * z[j]
    y = [z[i]/d[i] for i in range(n)]
    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = y[i]
        for j in range(i+1, n):
            x[i] -= l[j][i] * x[j]
    return x