共轭梯度法求解线性方程组

import numpy as np


def conjugate_gradient(A, b, x0=None, tol=1e-10, max_iter=1000):
    """共轭梯度法求解线性方程组 Ax = b

    参数:
        A: 系数矩阵(n×n)
        b: 右侧向量(n)
        x0: 初始猜测解(可选)
        tol: 收敛容差(默认1e-10)
        max_iter: 最大迭代次数(默认1000)

    返回:
        x: 方程组的近似解
    """
    n = len(b)
    # 初始化解向量
    if x0 is None:
        x = np.zeros(n)  # 如果没有初始猜测，使用零向量
    else:
        x = x0.copy()  # 复制初始猜测以避免修改原数组

    # 计算初始残差 r = b - Ax
    r = b - np.dot(A, x)
    p = r.copy()  # 初始搜索方向
    r_norm_sq = np.dot(r, r)  # 残差的平方范数

    # 开始迭代
    for k in range(max_iter):
        Ap = np.dot(A, p)  # 矩阵A与搜索方向p的乘积
        alpha = r_norm_sq / np.dot(p, Ap)  # 计算步长
        x = x + alpha * p  # 更新解
        r_new = r - alpha * Ap  # 计算新残差
        r_new_norm_sq = np.dot(r_new, r_new)  # 新残差的平方范数

        # 检查收敛性
        if np.sqrt(r_new_norm_sq) < tol:
            break

        # 计算共轭系数
        beta = r_new_norm_sq / r_norm_sq
        # 更新搜索方向
        p = r_new + beta * p
        # 更新残差及其范数
        r = r_new
        r_norm_sq = r_new_norm_sq

    return x


# 示例用法
if __name__ == "__main__":
    # 定义对称正定矩阵A和向量b
    A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
    b = np.array([1, 0, 1])
    # 调用共轭梯度法求解
    x = conjugate_gradient(A, b)
    print("近似解:", x)

过程详细版

import numpy as np

def conjugate_gradient(A, b, x0=None, tol=1e-10, max_iter=2):
    """共轭梯度法求解线性方程组 Ax = b

    参数:
        A: 系数矩阵(n×n，要求对称正定)
        b: 右侧向量(n)
        x0: 初始猜测解(可选)
        tol: 收敛容差(默认1e-10)
        max_iter: 最大迭代次数(默认1000)

    返回:
        x: 方程组的近似解
    """
    n = len(b)
    # 初始化解向量
    if x0 is None:
        x = np.zeros(n)  # 初始解设为[0,0,...,0]
        print(f"初始解 x⁽⁰⁾ = {x}")  # 上标(0)表示第0次迭代
    else:
        x = x0.copy()
        print(f"初始解 x⁽⁰⁾ = {x}")

    # 计算初始残差 r⁽⁰⁾ = b - Ax⁽⁰⁾
    r = b - np.dot(A, x)
    print(f"\n迭代0:")
    print(f"r⁽⁰⁾ = b - A·x⁽⁰⁾ = {b} - {np.dot(A, x)} = {r}")

    p = r.copy()  # 初始搜索方向 p⁽⁰⁾ = r⁽⁰⁾
    r_norm_sq = np.dot(r, r)  # ‖r⁽⁰⁾‖²
    print(f"p⁽⁰⁾ = r⁽⁰⁾ = {p}")
    print(f"‖r⁽⁰⁾‖² = {r_norm_sq:.6f}")

    # 开始迭代
    for k in range(max_iter):
        print(f"\n迭代{k+1}:")

        # 计算 A·p⁽ᵏ⁾
        Ap = np.dot(A, p)
        print(f"A·p⁽{k}⁾ = {A}·{p} = {Ap}")

        # 计算步长 αₖ = ‖r⁽ᵏ⁾‖² / (p⁽ᵏ⁾ᵀ·A·p⁽ᵏ⁾)
        alpha = r_norm_sq / np.dot(p, Ap)
        print(f"α_{k} = ‖r⁽{k}⁾‖²/(p⁽{k}⁾ᵀ·A·p⁽{k}⁾) = {r_norm_sq:.6f}/{np.dot(p, Ap):.6f} = {alpha:.6f}")

        # 更新解 x⁽ᵏ⁺¹⁾ = x⁽ᵏ⁾ + αₖ·p⁽ᵏ⁾
        x = x + alpha * p
        print(f"x⁽{k+1}⁾ = x⁽{k}⁾ + α_{k}·p⁽{k}⁾ = {x-alpha*p} + {alpha:.6f}·{p} = {x}")

        # 计算新残差 r⁽ᵏ⁺¹⁾ = r⁽ᵏ⁾ - αₖ·A·p⁽ᵏ⁾
        r_new = r - alpha * Ap
        print(f"r⁽{k+1}⁾ = r⁽{k}⁾ - α_{k}·A·p⁽{k}⁾ = {r} - {alpha:.6f}·{Ap} = {r_new}")

        # 计算新残差的范数平方 ‖r⁽ᵏ⁺¹⁾‖²
        r_new_norm_sq = np.dot(r_new, r_new)
        print(f"‖r⁽{k+1}⁾‖² = {r_new_norm_sq:.6f}")

        # 检查收敛性
        if np.sqrt(r_new_norm_sq) < tol:
            print(f"\n收敛于第{k+1}次迭代，残差范数√‖r⁽{k+1}⁾‖² = {np.sqrt(r_new_norm_sq):.6f} < {tol}")
            break

        # 计算共轭系数 βₖ = ‖r⁽ᵏ⁺¹⁾‖² / ‖r⁽ᵏ⁾‖²
        beta = r_new_norm_sq / r_norm_sq
        print(f"β_{k} = ‖r⁽{k+1}⁾‖²/‖r⁽{k}⁾‖² = {r_new_norm_sq:.6f}/{r_norm_sq:.6f} = {beta:.6f}")

        # 更新搜索方向 p⁽ᵏ⁺¹⁾ = r⁽ᵏ⁺¹⁾ + βₖ·p⁽ᵏ⁾
        p = r_new + beta * p
        print(f"p⁽{k+1}⁾ = r⁽{k+1}⁾ + β_{k}·p⁽{k}⁾ = {r_new} + {beta:.6f}·{p-beta*p} = {p}")

        # 更新残差及其范数
        r = r_new
        r_norm_sq = r_new_norm_sq

    return x

# 示例用法
if __name__ == "__main__":
    # 定义对称正定矩阵A和向量b

    # example 01
    # A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
    # b = np.array([1, 0, 1])

    # example 02
    A = np.array([[4, 3, 0], [3, 4, -1], [0, -1, 4]])
    b = np.array([ 24, 30,-24])

    # example 03
    # A = np.array([[6, 2], [3,2]])
    # b = np.array([0,-1])
    print("求解线性方程组 Ax = b")
    print(f"其中 A = \n{A}")
    print(f"b = {b}\n")

    # 调用共轭梯度法求解
    x = conjugate_gradient(A, b)
    print("\n最终近似解:", x)

    # 验证解的正确性
    print("\n验证 A·x = ", np.dot(A, x))
    print("应与 b = ", b, "非常接近")