Numerical calculation method chapter6 Iterative method for solving linear equations

• Numerical calculation method Chapter6. Iterative methods for solving linear equations
  • 0. Problem description
  • 1. Jacobi iteration
    • 1. Jacobi Iterative method
    • 2. Jacobi Iterative matrix
    • 3. Jacobi Iterative convergence condition
    • 4. python Pseudo code implementation
  • 2. Gauss-Seidel iteration
    • 1. Gauss-Seidel Iterative method
    • 2. Gauss-Seidel Iterative matrix
    • 3. Gauss-Seidel Iterative convergence condition
    • 4. Pseudo code implementation
  • 3. Relaxation iteration
    • 1. Relaxation iteration method
    • 2. Relaxation iteration matrix
    • 3. Convergence condition of relaxation iteration
    • 4. Pseudo code implementation
  • 4. Inverse matrix calculation

0. Problem description

The problems to be solved in this chapter are the same as those in the previous chapter , Still $n$ The problem of solving a system of linear equations .

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

however , In the previous chapter, the linear transformation of the matrix is transformed into a triangular matrix or a diagonal matrix that can be solved quickly , The calculation result is accurate .

The idea of this chapter is to solve the problem $Ax=y$ convert to $x=A^{\prime }x$ The iterative form of , thus , Then we can give a group of iterated algebras $x_{k+1}=A^{\prime }{x}_k$.

here , If $x_k$ Satisfy the convergence condition , that $x_k$ Will converge to $x=A^{\prime }x$ In a set of solutions of , The above questions can also be answered .

1. Jacobi iteration

1. Jacobi Iterative method

Linear transformation of the original equation , We obviously have ：

$$\{\begin{array}{rlrlrl}{x}_1=& & -\frac{a_{12}}{a_{11}}{x}_2& -...& -\frac{a_{1n}}{a_{11}}{x}_n& +\frac{y_1}{a_{11}}\\ x_2=& -\frac{a_{21}}{a_{22}}{x}_2& & -...& -\frac{a_{2n}}{a_{22}}{x}_n& +\frac{y_2}{a_{22}}\\ ...& & & & & \\ x_n=& -\frac{a_{n1}}{a_{nn}}{x}_1& -\frac{a_{n2}}{a_{nn}}{x}_2& -...& & +\frac{y_n}{a_{nn}}\end{array}$$

We write it in matrix form, that is ：

$$\left(\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_n\end{array}\right)=\left(\begin{array}{cccc}0& -\frac{a_{12}}{a_{11}}& ...& -\frac{a_{1n}}{a_{11}}\\ -\frac{a_{21}}{a_{22}}& 0& ...& -\frac{a_{2n}}{a_{22}}\\ ...& ...& ...& ...\\ -\frac{a_{n1}}{a_{nn}}& -\frac{a_{n2}}{a_{nn}}& ...& 0\end{array}\right)\left(\begin{array}{c}{x}_1\\ x_2\\ ...\\ x_n\end{array}\right)+\left(\begin{array}{c}\frac{y_1}{a_{11}}\\ \frac{y_2}{a_{22}}\\ ...\\ \frac{y_n}{a_{nn}}\end{array}\right)$$

We simplify it into a symbolic expression that is ：

$$x=Bx+g$$

Based on this , We can give it Jacobi The iterative formula is as follows ：

$$x_{k+1}=B{x}_k+g$$

2. Jacobi Iterative matrix

We give a more rigorous data expression as follows , Make ：

$$D=\left(\begin{array}{cccc}{a}_{11}& & & \\ & a_{22}& & \\ & & ...& \\ & & & a_{nn}\end{array}\right)$$

Then there are ：

$$Ax=(D+A-D)x=y$$

Then we can export ：

$$Dx=(D-A)x+y$$

Double left $D^{-1}$ That is to say ：

$$x=D^{-1}(D-A)x+D^{-1}y$$

remember $B=D^{-1}(D-A)$,$g=D^{-1}y$, You can get exactly the same content as in the previous section .

Iterative matrix $B$ That is to say Jacobi Iterative matrix .

3. Jacobi Iterative convergence condition

Jacobi The necessary and sufficient condition for iterative convergence is ：

Iterative matrix $B$ The spectral radius of $\rho (B)<1$.

or ：

• For iterative matrices $B$ All eigenvalues of ${\lambda }_i$, Both $|{\lambda }_i|<1$, Or equivalent $max(|{\lambda }_i|)<1$.

however , In some specific cases , The above necessary and sufficient conditions can be appropriately relaxed .

To be specific , There are the following theorems ：

Theorem 6.1
If the equations $Ax=y$ The coefficient matrix of $A$, Meet one of the following conditions , Then Jacobi Iterative convergence ：
(1) A Is a row diagonally dominant matrix , namely $|a_{ii}|>{\sum }_{j\ne i}|a_{ij}|$,$i=1,2,...,n$;
(2) A Is a column diagonally dominant matrix , namely $|a_{jj}|>{\sum }_{i\ne j}|a_{ij}|$,$j=1,2,...,n$;

4. python Pseudo code implementation

Last , We give Jacobi Iterative python The pseudocode is as follows ：

def jacobi_iter(A, y, epsilon=1e-6):
    n = len(A)
    B = [[0 for _ in range(n)] for _ in range(n)]
    g = [0 for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if j == i:
                continue
            B[i][j] = -A[i][j] / A[i][i]
        g[i] = y[i] / A[i][i]
    
    x = [0 for _ in range(n)]
    for _ in range(10**6):
        z = [0 for _ in range(n)]
        for i in range(n):
            z[i] += g[i]
            for j in range(n):
                z[i] += B[i][j] * x[j]
        if max(abs(z[i] - x[i]) for i in range(n)) < epsilon:
            break
        x = z
    return z

2. Gauss-Seidel iteration

1. Gauss-Seidel Iterative method

Gauss-Seidel Iterative equations and the above Jacobi Iterations are actually very similar , The only difference is that Jacobi Iteration is based on $x^{(k)}$ Update iteratively for the whole each time , and Guass-Seidel Iteration is to calculate each $x_i^{(k)}$ Is to use all that have been iteratively calculated $x_j^{(k)}$ Result .

That is , Each element is updated iteratively .

Expressed by formula, it is ：

$$\{\begin{array}{rlrlrl}{x}_1^{(k+1)}=& & -\frac{a_{12}}{a_{11}}{x}_2^{(k)}& -...& -\frac{a_{1n}}{a_{11}}{x}_n^{(k)}& +\frac{y_1}{a_{11}}\\ x_2^{(k+1)}=& -\frac{a_{21}}{a_{22}}{x}_1^{(k+1)}& & -...& -\frac{a_{2n}}{a_{22}}{x}_n^{(k)}& +\frac{y_2}{a_{22}}\\ ...& & & & & \\ x_n^{(k+1)}=& -\frac{a_{n1}}{a_{nn}}{x}_1^{(k+1)}& -\frac{a_{n2}}{a_{nn}}{x}_2^{(k+1)}& -...& & +\frac{y_n}{a_{nn}}\end{array}$$

2. Gauss-Seidel Iterative matrix

alike , We give a more rigorous mathematical derivation as follows ：

$$\begin{array}{rl}A& =D+L+U\\ & =\left(\begin{array}{cccc}{a}_{11}& & & \\ & a_{22}& & \\ & & ...& \\ & & & a_{nn}\end{array}\right)+\left(\begin{array}{cccc}0& & & \\ a_{21}& 0& & \\ ...& & & \\ a_{n1}& a_{n2}& ...& 0\end{array}\right)\left(\begin{array}{cccc}0& a_{12}& ...& a_{1n}\\ & 0& ...& a_{2n}\\ & & ...& \\ & & & 0\end{array}\right)\end{array}$$

thus , We have ：

$$Ax=(D+L+U)x=y$$

or ：

$$(D+L)x=-Ux+y$$

Multiply both sides by $(D+l{)}^{-1}$ That is to say ：

$$x=-(D+L{)}^{-1}Ux+(D+L{)}^{-1}y$$

Make $S=-(D+L{)}^{-1}U$,$f=(D+L{)}^{-1}y$, We can write Gauss-Seidel Iterative formula ：

$$x^{(k+1)}=S{x}^{(k)}+f$$

It is called iterative matrix $S$ That is to say Gauss-Seidel Iterative matrix .

Of course , As mentioned in the previous section , In fact, it is not necessary to calculate $S$ and $f$, Just according to the previous $Jacobi$ The iteration matrix can be updated in real time .

3. Gauss-Seidel Iterative convergence condition

alike , We give the book about Gauss-Seidel The convergence conditions of the iteration are as follows ：

Theorem 6.2
If the coefficient matrix of the system of equations is diagonally dominant in rows or columns , be Gauss-Seidel Iterative convergence .
Theorem 6.3
If the coefficient matrix of the system of equations $A$ Is a symmetric positive definite matrix , be Gauss-Seidel Iterative convergence .

4. Pseudo code implementation

alike , We use it python The pseudo code given is as follows ：

def gauss_seidel_iter(A, y, epsilon=1e-6):
    n = len(A)
    B = [[0 for _ in range(n)] for _ in range(n)]
    g = [0 for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if j == i:
                continue
            B[i][j] = -A[i][j] / A[i][i]
        g[i] = y[i] / A[i][i]
    
    x = [0 for _ in range(n)]
    cnt = 0
    for _ in range(10**6):
        for i in range(n):
            z = g[i]
            for j in range(n):
                z += B[i][j] * x[j]
            if abs(z-x[i]) < epsilon:
                cnt += 1
            else:
                cnt = 0
            x[i] = z

            if cnt >= n:
                break
        if cnt >= n:
            break
    return x

3. Relaxation iteration

1. Relaxation iteration method

The prototype of relaxation iteration is still the same as before Jacobi iteration , however , and Gauss-Seidel The real-time parameter update of the iteration is different , Relaxation iteration is right here Jacobi Iterative batch update and Gauss-Seidel Iterative real-time updates take a compromise , Through a super parameter $w$ To control the parameter update speed .

To be specific , That is to say ：

$$\{\begin{array}{rl}{x}_1^{(k+1)}& =(1-w)x_1^{(k)}+w(b_{12}{x}_2^{(k)}+b_{13}{x}_3^{(k)}+...+b_{1n}{x}_n^{(k)}+g_1)\\ x_2^{(k+1)}& =(1-w)x_2^{(k)}+w(b_{21}{x}_1^{(k+1)}+b_{23}{x}_3^{(k)}+...+b_{2n}{x}_n^{(k)}+g_2)\\ ...& \\ x_n^{(k+1)}& =(1-w)x_n^{(k)}+w(b_{n1}{x}_1^{(k+1)}+b_{n2}{x}_2^{(k+1)}+...+b_{n,n-1}{x}_{n-1}^{(k+1)}+g_n)\end{array}$$

2. Relaxation iteration matrix

alike , We give the mathematical expression of relaxation iteration as follows ：

$$x^{(k+1)}=(1-w)x^{(k)}+w(\tilde{L}{x}^{(k+1)}+\tilde{U}{x}^{(k)}+g)$$

Modelled on the Gauss-Seidel iteration , We can also give the iteration matrix format as follows ：

$$x^{(k+1)}=S_w{x}^{(k)}+f$$

among ,

$$\{\begin{array}{rl}{S}_w& =(I+w{D}^{-1}L{)}^{-1}[(1-w)I-w{D}^{-1}U]\\ f& =w(I+w{D}^{-1}L{)}^{-1}g\end{array}$$

3. Convergence condition of relaxation iteration

There are two theorems to judge the convergence of relaxation iteration ：

Theorem 6.4
The necessary condition for the convergence of relaxation iteration is $0<w<2$;
Theorem 6.5
if $A$ For symmetric positive definite matrix , Then when $0<w<2$ when , Relaxation iteration has constant convergence ;

Special , When $0<w<1$ when , Call it Sub relaxation iteration ; When $1<w<2$ when , Call it Over relaxation iteration ; When $w=1$ when , Iterations degenerate into Gauss-Seidel iteration .

4. Pseudo code implementation

Last , We also give the pseudo code implementation of relaxation iteration as follows ：

def loose_iter(A, y, w=0.5, epsilon=1e-6):
    n = len(A)
    B = [[0 for _ in range(n)] for _ in range(n)]
    g = [0 for _ in range(n)]
    for i in range(n):
        for j in range(n):
            if j == i:
                continue
            B[i][j] = -A[i][j] / A[i][i]
        g[i] = y[i] / A[i][i]
    
    x = [0 for _ in range(n)]
    cnt = 0
    for _ in range(10**6):
        for i in range(n):
            z = (1-w) * x[i] + w * g[i]
            for j in range(n):
                z += w * B[i][j] * x[j]
            if abs(z-x[i]) < epsilon:
                cnt += 1
            else:
                cnt = 0
            x[i] = z

            if cnt >= n:
                break
        if cnt >= n:
            break
    return x

4. Inverse matrix calculation

Last , Let's look at the inverse matrix calculation .

In principle, the calculation of inverse matrix is actually a special application of solving linear equations mentioned above , As a matter of fact $n$ Then we can get our inverse matrix by splicing the solutions of the unit vectors .

We won't go into details here , Just give python The pseudo code is implemented as follows ：

def cal_inverse_matrix(A):
    n = len(A)
    res = [[0 for _ in range(n)] for _ in range(n)]
    for j in range(n):
        y = [0 for _ in range(n)]
        y[j] = 1
        x = gauss_seidel_iter(A, y)
        for i in range(n):
            res[i][j] = x[i]
    return res