编制 Jacobi 迭代法、 Gauss-Seidel 迭代法和超松弛(SOR)迭代法解方程组 Ax=b 的程序

Jacobi 迭代法、Gauss-Seidel 迭代法和SOR(Successive Over-Relaxation)迭代法都是数值线性代数中用于求解矩阵方程 \(Ax = b\) 的迭代算法，它们通常应用于大型稀疏系统，因为直接求解往往效率低下。

以下是简化的 Python 示例，使用 NumPy 库来演示这些方法：

import numpy as np

def jacobi(A, b, tolerance=1e-6, max_iterations=100):
    n = len(b)
    x = np.zeros(n)
    D = np.diagflat(A.diagonal())
    L = D - A
    U = D + A.T
    for _ in range(max_iterations):
        old_x = x.copy()
        x = (b / D).T + (L @ x)
        if np.allclose(x, old_x, atol=tolerance):
            break
    return x

def gauss_seidel(A, b, tolerance=1e-6, max_iterations=100):
    n = len(b)
    x = np.zeros(n)
    for _ in range(max_iterations):
        for i in range(n):
            sum_row = sum(A[i, :i] * x[:i])
            sum_col = sum(A[i, i+1:] * x[i+1:])
            x[i] = (b[i] - sum_row - sum_col) / A[i, i]
        if np.allclose(x, x_old, atol=tolerance):
            break
    return x

def sor(A, b, omega=1.5, tolerance=1e-6, max_iterations=100):
    n = len(b)
    x = np.zeros(n)
    D = np.diagflat(A.diagonal())
    L = omega * (D - A)
    U = (1 - omega) * (D + A.T)
    for _ in range(max_iterations):
        old_x = x.copy()
        x = ((omega * b) + (U @ x)) / (omega * D + U.sum(axis=0))
        if np.allclose(x, old_x, atol=tolerance):
            break
    return x

# 示例用法
A = np.array([[4, 1], [1, 3]])
b = np.array([7, 11])
solution_jacobi = jacobi(A, b)
solution_gauss_seidel = gauss_seidel(A, b)
solution_sor = sor(A, b)

print(f"Jacobi method solution: {solution_jacobi}")
print(f"Gauss-Seidel method solution: {solution_gauss_seidel}")
print(f"SOR method solution: {solution_sor}")

在这个例子中，我们首先创建了矩阵 A 和向量 b，然后分别计算三种方法的解。`max_iterations` 和 `tolerance` 参数可以用来控制迭代次数和收敛精度。每个方法内部都包含了一个循环，在满足停止条件（通常是相邻迭代步长足够小）时跳出循环。