解线性代数迭代法代码实现

jacobi 迭代法

import numpy as np
from tabulate import tabulate


def jacobi_method(A, b, x0, num_iterations=20):
    # 获取系数矩阵对角矩阵
    D = np.diag(np.diag(A))
    #     计算下三角和上三角矩阵的和
    L_plus_U = A - D
    #     求对角矩阵的逆
    D_inv = np.linalg.inv(D)
    # 初始化解向量
    x = x0

    # 用于保存每次迭代的结果
    steps = []

    for i in range(num_iterations):
        x_new = D_inv@(b - L_plus_U@x)
        steps.append([f"迭代{i + 1}次"] + [f"{val:4f}" for val in x_new])
        x = x_new
    return x, steps


A = np.array([
    [5, -1, 1],
    [2, 8, -1],
    [-1, 1, 4]
], dtype=float)
b = np.array([10, 11, 3], dtype=float)
x0 = [0, 0, 0]


solution, steps = jacobi_method(A, b, x0)
headers = ["迭代", "x", "y", "z"]
print("jacobi 迭代法详细结果:")
print(tabulate(steps, headers=headers))
# print("迭代",num_iterations,"次后的解:")

print("\n最终解为:\t")
print("x=", solution[0])
print("y=", solution[1])
print("z=", solution[2])


# 验证：使用NumPy的直接解法求解线性方程组Ax=b
exact_solution = np.linalg.solve(A, b)
print("NumPy 直接解:", exact_solution)

gauss-seidel迭代法

import numpy as np
from tabulate import tabulate


# Gauss-Seidel 迭代法（矩阵形式）
def gauss_seidel_matrix(A, b, x0, max_iterations=20, tolerance=1e-6):
    """
    使用Gauss-Seidel迭代法求解线性方程组Ax=b。

    参数:
    A (numpy.ndarray): 系数矩阵，形状为(n, n)。
    b (numpy.ndarray): 常数项向量，形状为(n,)。
    x0 (numpy.ndarray): 初始猜测解向量，形状为(n,)。
    max_iter (int): 最大迭代次数。
    tol (float): 收敛容差。

    返回:
    tuple: 包含最终解向量和迭代过程记录的元组。
    """
    # 提取上三角矩阵U（不包括对角线）
    U = np.triu(A, k=1)  # 上三角矩阵

    # 计算D + L矩阵（A - U）
    D_plus_L = A - U  # D + L 矩阵

    # 初始化解向量x为初始猜测值x0
    x = x0

    # 初始化用于记录每一步迭代结果的列表
    steps = []

    # 开始迭代过程
    for k in range(max_iterations):
        # 计算新的解向量x_new
        x_new = np.linalg.inv(D_plus_L) @ (b - U @ x.copy())

        # 记录当前迭代结果，保留四位小数
        steps.append([k + 1] + [f"{val:.4f}" for val in x_new])

        # 检查是否满足收敛条件
        if np.linalg.norm(x_new - x) < tolerance:
            print(f"\n在第 {k + 1} 次迭代后收敛")
            return x_new, steps

        # 更新解向量x为新计算的x_new
        x = x_new

    # 如果达到最大迭代次数仍未收敛，输出提示信息并返回当前解和迭代过程
    print("\n达到最大迭代次数，未完全收敛")
    return x, steps


# 定义系数矩阵A
A = np.array([
    [5, -1, 1],
    [2, 8, -1],
    [-1, 1, 4]
], dtype=float)

# 定义常数项向量b
b = np.array([10, 11, 3], dtype=float)

# 初始猜测解向量x0
x0 = np.array([0, 0, 0], dtype=float)


# 运行Gauss-Seidel迭代法
solution, steps = gauss_seidel_matrix(A, b, x0)

# 定义表格头部
headers = ["迭代次数", "x", "y", "z"]

# 打印Gauss-Seidel迭代法的详细结果
print("Gauss-Seidel 迭代法详细结果:")
print(tabulate(steps, headers=headers))

# 输出最终解
print("\n最终解:", solution)

# 验证：使用NumPy的直接解法求解线性方程组Ax=b
exact_solution = np.linalg.solve(A, b)
print("NumPy 直接解:", exact_solution)

SOR迭代法

import numpy as np


# SOR 迭代法（矩阵形式）
def sor_method(A, b, x0, omega, max_iter, tol):
    """
    使用SOR迭代法求解线性方程组Ax=b。

    参数:
    A (numpy.ndarray): 系数矩阵，形状为(n, n)。
    b (numpy.ndarray): 常数项向量，形状为(n,)。
    x0 (numpy.ndarray): 初始猜测解向量，形状为(n,)。
    omega (float): 松弛因子，通常在 (1, 2) 之间选择。
    max_iter (int): 最大迭代次数。
    tol (float): 收敛容差。

    返回:
    numpy.ndarray: 最终解向量。
    """
    n = A.shape[0]
    D = np.diag(np.diag(A))  # 对角矩阵
    L = np.tril(A, k=-1)  # 严格下三角矩阵
    U = np.triu(A, k=1)  # 严格上三角矩阵

    for k in range(max_iter):
        # 计算新的解向量 x_new
        x_new = np.zeros_like(x0)
        for i in range(n):
            sigma1 = L[i, :i] @ x_new[:i]  # 已更新的部分
            sigma2 = U[i, i:] @ x0[i:]  # 未更新的部分
            x_new[i] = (1 - omega) * x0[i] + omega / A[i, i] * (b[i] - sigma1 - sigma2)

        # 输出当前迭代结果
        print(f"迭代 {k + 1}: {x_new}")

        # 检查收敛性
        if np.linalg.norm(x_new - x0) < tol:
            print(f"\n在第 {k + 1} 次迭代后收敛")
            return x_new

        x0 = x_new  # 更新解向量

    print("\n达到最大迭代次数，未完全收敛")
    return x0


A = np.array([
    [10, 3, 1],
    [2, -10, 3],
    [1, 3, 10]
], dtype=float)
b = np.array([14, -5, 14], dtype=float)

# 初始猜测
x0 = np.array([0, 0, 0], dtype=float)

# 参数
omega = 1.25  # 松弛因子，通常在 (1, 2) 之间选择
max_iterations = 100
tolerance = 1e-8

# 运行 SOR
solution = sor_method(A, b, x0, omega, max_iterations, tolerance)

# 输出最终结果
print("\n最终解:", solution)

# 验证：用 NumPy 的直接解法
exact_solution = np.linalg.solve(A, b)
print("NumPy 直接解:", exact_solution)