Solve linear equations using Python to implement iterative methods (Jacobi iteration, Gauss-Seidel iteration, relaxation iteration)

1. Jacobi iteration

The Jacobi iteration method is an iterative algorithm for solving systems of linear equations. It belongs to the direct iteration method among the iterative methods, and it approximates the solution of the linear equation system by continuously updating the solution vector iteratively.

The basic concepts of Jacobi iteration method are as follows:

Given the coefficient matrix A of the system of linear equations and the constant vector b on the right side.
Diagonally decompose the coefficient matrix A to obtain three matrices D, L and U, where D is the diagonal matrix of A, L is the strict lower triangular matrix of A (that is, the elements below the main diagonal are 0), and U is A strictly upper triangular matrix (that is, the elements above the main diagonal are 0).
Initialize the solution vector x to any initial value (can be an all-0 vector).
Iteratively update the solution vector x until the convergence condition is met:
- Calculate the residual vector r = b – Ax, where A is the coefficient matrix and x is the current solution vector.
- Use the update formula: x = D^(-1) * (b – (L + U)x), where D^(-1) represents the inverse matrix of D.
- Repeat the above two steps until the norm of the residual vector r is less than the set convergence threshold or the maximum number of iterations is reached.
Returns the approximate solution vector x.

The convergence and speed of the Jacobi iteration method are related to the characteristics of the coefficient matrix A. When the coefficient matrix A is diagonally dominant or strictly diagonally dominant, the Jacobi iteration method has convergence. However, in some cases, the convergence speed of Jacobi iteration method is slow, so it can be combined with other iteration methods, such as Gauss-Seidel iteration method or successive super-relaxation method, to improve the iteration effect.

"""
@Time: 2023/10/19 0019 17:11
@Auth:yeqc

        Jacobi iteration method
"""
import numpy as np


def Jacobi(A, b, max_iter=100, epsilon=1e-10):
    n = len(b)
    x = np.zeros(n) # initial solution
    x_new = np.zeros(n) # Store the updated solution vector

    for iter in range(max_iter):
        for i in range(n):
            sum_ax = np.dot(A[i], x) - A[i][i] * x[i] # Calculate the sum excluding the diagonal elements
            x_new[i] = (b[i] - sum_ax) / A[i][i] # Update the i-th component of the solution vector

        # Determine the termination condition
        if np.linalg.norm(x_new - x) < epsilon:
            break

        x = np.copy(x_new) # Update solution vector
    return x


# Example
A = np.array([[9, -1, -1], [-1, 8, 0], [-1, 0, 9]])

b = np.array([7, 7, 8])

x = Jacobi(A, b)
print("Solution:", x)

2. Gauss-Seidel iteration

The Gauss-Seidel iteration method is an iterative algorithm for solving linear equations. It is an improved version of the Jacobi iteration method. This method approximates the solution of a system of linear equations by continuously iteratively updating each component of the solution vector.

The basic concepts of the Gauss-Seidel iteration method are as follows:

Given the coefficient matrix A of the system of linear equations and the constant vector b on the right side.
Diagonally decompose the coefficient matrix A to obtain three matrices D, L and U, where D is the diagonal matrix of A, L is the strict lower triangular matrix of A (that is, the elements below the main diagonal are 0), and U is A strictly upper triangular matrix (that is, the elements above the main diagonal are 0).
Initialize the solution vector x to any initial value (can be an all-0 vector).
Iteratively update the solution vector x until the convergence condition is met:
- For the i-th unknown number, calculate its new value: x(i) = (1 / a(i, i)) * (b(i) – Σ(a(i, j) * x(j), j=1 to i-1) – Σ(a(i, j) * x(j), j=i + 1 to n)), where a(i, j) represents the i-th row and j-th column element of the coefficient matrix A.
- By continuously updating each component in the solution vector x, a new solution vector is obtained.
- Repeat the above two steps until the convergence condition is met: the norm of the residual vector is less than the set convergence threshold or the maximum number of iterations is reached.
Returns the approximate solution vector x.

The improvement of the Gauss-Seidel iteration method compared to the Jacobi iteration method is that in each iteration, it uses the new components of the updated solution vector to calculate the new value of the next unknown, thus speeding up the convergence speed. However, similar to the Jacobi iteration method, the Gauss-Seidel iteration method also needs to meet certain convergence conditions to obtain effective results.

"""
@Time: 2023/10/19 0019 17:11
@Auth:yeqc

        Gauss-Seidel Gauss-Seidel iteration
"""
import numpy as np


def Jacobi(A, b, max_iter=100, epsilon=1e-10):
    n = len(b)
    x = np.zeros(n) # initial solution
    x_new = np.zeros(n) # Store the updated solution vector

    for iter in range(max_iter):
        for i in range(n):
            sum_ax = np.dot(A[i], x_new) - A[i][i] * x[i] # Calculate the sum excluding the diagonal elements, note: Gauss-Seidel, use x_new instead of x here
            x_new[i] = (b[i] - sum_ax) / A[i][i] # Update the i-th component of the solution vector

        # Determine the termination condition
        if np.linalg.norm(x_new - x) < epsilon:
            break

        x = np.copy(x_new) # Update solution vector
    return x


# Example
A = np.array([[9, -1, -1], [-1, 8, 0], [-1, 0, 9]])

b = np.array([7, 7, 8])

x = Jacobi(A, b)
print("Solution:", x)

3. Relaxation iteration

Relaxation Iteration is an iterative algorithm used to solve linear equations. It is often used to solve sparse matrices and large-scale linear equations. It introduces a relaxation factor based on the Gauss-Seidel iteration method to accelerate the convergence process.

The basic concept of relaxed iteration is as follows:

Given the coefficient matrix A of the system of linear equations and the constant vector b on the right side.
Diagonally decompose the coefficient matrix A to obtain three matrices D, L and U, where D is the diagonal matrix of A, L is the strict lower triangular matrix of A (that is, the elements below the main diagonal are 0), and U is A strictly upper triangular matrix (that is, the elements above the main diagonal are 0).
Initialize the solution vector x to any initial value (can be an all-0 vector).
Iteratively update the solution vector x until the convergence condition is met:
- For the i-th unknown number, calculate its new value: x(i) = (1 – w) * x(i) + (w / a(i, i)) * (b(i) – Σ(a(i, j) * x(j), j=1 to i-1) – Σ(a(i, j) * x(j), j=i + 1 to n)), where a(i, j) represents the coefficient The i-th row and j-th column element of matrix A, w is the relaxation factor (the value range is 0 < w < 2).
- By continuously updating each component in the solution vector x, a new solution vector is obtained.
- Repeat the above two steps until the convergence condition is met: the norm of the residual vector is less than the set convergence threshold or the maximum number of iterations is reached.
Returns the approximate solution vector x.

The relaxation iteration method introduces a relaxation factor w, which is used to adjust the amplitude of each iteration update. When the relaxation factor is close to 1, the algorithm converges slowly; when the relaxation factor is close to 0 or 2, the algorithm converges quickly. Appropriate selection of the relaxation factor can significantly improve the convergence speed of the algorithm, but too large or too small a relaxation factor may cause the algorithm to diverge.

It should be noted that in the relaxation iteration method, the selection of relaxation factor is an important issue and needs to be debugged and optimized according to specific problems. Typically, appropriate relaxation factors can be determined through experimentation and experience.

"""
@Time: 2023/10/19 0019 17:11
@Auth:yeqc

        Relaxation iteration
"""
import numpy as np


def relaxation_iteration(A, b, x0, max_iter=100, tol=1e-6, omega=1.0):
    n = len(A)
    x = x0.copy()

    for k in range(max_iter):
        for i in range(n):
            x[i] = (1 - omega) * x[i] + (omega / A[i, i]) * (
                        b[i] - np.dot(A[i, :i], x[:i]) - np.dot(A[i, i + 1:], x0[i + 1:]))

        if np.linalg.norm(x - x0) < tol:
            break

        x0 = x.copy()

    return x


# Example
A = np.array([[9, -1, -1], [-1, 8, 0], [-1, 0, 9]])

b = np.array([7, 7, 8])
x0 = np.zeros(3) #Initial solution vector
omega = 1.2 # relaxation factor
x = relaxation_iteration(A, b, x0, omega=omega)

print("Solution:", x)