超松弛迭代法求解线性方程组

import numpy as np
import matplotlib.pyplot as plt


def sor_solver(A, b, omega, x0, tol=1e-5, max_iterations=200):
    """
    使用SOR方法求解线性方程组 Ax = b.

    参数:
        A (np.array): 系数矩阵.
        b (np.array): 常数向量.
        omega (float): 松弛因子 (0 < omega < 2).
        x0 (np.array): 初始猜测向量.
        tol (float): 收敛容差.
        max_iterations (int): 最大迭代次数.

    返回:
        x (np.array): 方程组的解.
        iterations (int): 达到收敛所用的迭代次数.
        errors (list): 每次迭代的误差历史 (无穷范数).
    """
    n = len(b)
    x = x0.copy()
    # 文档中例子的精确解，用于计算真实误差
    x_exact = np.array([1.1, 1.2, 1.3])
    errors = []

    for k in range(max_iterations):
        x_old = x.copy()

        # 应用SOR迭代公式
        # x_i^(k+1) = (1-ω)x_i^(k) + (ω/a_ii)*(b_i - Σ(j<i) a_ij*x_j^(k+1) - Σ(j>i) a_ij*x_j^(k))
        for i in range(n):
            sigma1 = np.dot(A[i, :i], x[:i])
            sigma2 = np.dot(A[i, i + 1:], x_old[i + 1:])
            x[i] = (1 - omega) * x_old[i] + (omega / A[i, i]) * (b[i] - sigma1 - sigma2)

        # 记录与精确解的误差（使用无穷范数）
        error = np.linalg.norm(x - x_exact, np.inf)
        errors.append(error)

        # 检查迭代是否收敛（基于相邻两次迭代结果的差异）
        if np.linalg.norm(x - x_old, np.inf) < tol:
            return x, k + 1, errors

    return x, max_iterations, errors


# --- 主程序 ---
# 从文档【例6.8】中获取方程组信息
A = np.array([
    [10.0, -1.0, -2.0],
    [-1.0, 10.0, -2.0],
    [-1.0, -1.0, 5.0]
])

b = np.array([7.2, 8.3, 4.2])

# 初始猜测向量
x0 = np.zeros_like(b)

# 从文档表6.4中选择不同的松弛因子进行测试
omegas_to_test = [0.8, 1.0, 1.055, 1.2, 1.4, 1.8]

# --- 可视化 ---
# 设置matplotlib以支持中文显示
# 请确保你的环境中已安装支持中文的字体，如'SimHei', 'Microsoft YaHei'等
try:
    plt.rcParams['font.family'] = 'serif'
    plt.rcParams['font.serif'] = ['SimHei','DejaVu Serif', 'Times New Roman', 'Computer Modern Roman']
    plt.rcParams['mathtext.fontset'] = 'cm'
    plt.rcParams['axes.unicode_minus'] = False
except:
    print("未找到中文字体'SimHei', 图表标签可能无法正常显示。")

plt.figure(figsize=(12, 9))

# 对每个omega值运行SOR算法并绘图
for omega in omegas_to_test:
    solution, iterations, errors = sor_solver(A, b, omega, x0)
    label_text = f'ω = {omega} ({iterations} 次迭代)'
    if omega == 1.0:
        label_text = f'高斯-赛德尔 (ω = 1.0) ({iterations} 次迭代)'
    if omega == 1.055:
        label_text = f'最佳因子 ω = {omega} ({iterations} 次迭代)'

    plt.plot(range(1, len(errors) + 1), errors, marker='o', linestyle='-', markersize=4, label=label_text)

# --- 绘图设置 ---
plt.yscale('log')  # Y轴使用对数刻度以更好地观察误差下降
plt.xlabel('迭代次数 (k)', fontsize=14)
plt.ylabel('误差 (无穷范数 ||x^(k) - x*||∞)', fontsize=14)
plt.title('不同松弛因子(ω)对SOR算法收敛速度的影响', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True, which="both", ls="--")
plt.show()

SOR-V2

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties, findfont
from scipy.sparse import diags


# ========== 自动检测中文字体 ==========
def setup_chinese_font():
    """自动检测并设置可用的中文字体"""
    chinese_fonts = ["Microsoft YaHei", "FangSong", "KaiTi", "SimHei", "SimSun"]

    found_font = None
    for font in chinese_fonts:
        try:
            if findfont(FontProperties(family=[font])) not in ['', None]:
                found_font = font
                break
        except:
            continue

    if found_font:
        plt.rcParams['font.sans-serif'] = [found_font]
        plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
        print(f"✅ 使用中文字体: {found_font}")
    else:
        print("⚠️ 未找到可用中文字体，请手动安装 SimHei.ttf 或指定字体路径。")


# 初始化中文字体
setup_chinese_font()


# ========== SOR迭代求解器 ==========
def sor_solver(A, b, omega, x0, tol=1e-5, max_iterations=200):
    """
    使用SOR(逐次超松弛)方法求解线性方程组 Ax = b

    参数:
        A: 系数矩阵 (n x n)
        b: 右侧向量 (n x 1)
        omega: 松弛因子 (0 < omega < 2)
        x0: 初始猜测解 (n x 1)
        tol: 收敛容差
        max_iterations: 最大迭代次数

    返回:
        x: 解向量
        residuals: 残差历史记录
        iterations: 实际迭代次数
    """
    n = len(b)
    x = x0.copy()
    residuals = []

    for k in range(max_iterations):
        x_new = x.copy()
        for i in range(n):
            sigma = 0.0
            # 前向替换 (已经更新的分量)
            for j in range(i):
                sigma += A[i, j] * x_new[j]
            # 后向替换 (尚未更新的分量)
            for j in range(i + 1, n):
                sigma += A[i, j] * x[j]

            # SOR更新公式
            x_new[i] = (1 - omega) * x[i] + (omega / A[i, i]) * (b[i] - sigma)

        # 计算残差
        residual = np.linalg.norm(A @ x_new - b)
        residuals.append(residual)

        # 检查收敛
        if residual < tol:
            break

        x = x_new

    return x_new, np.array(residuals), k + 1


# ========== 测试案例生成 ==========
def generate_test_case(n=100, k=3):
    """
    生成对角占优的三对角矩阵测试案例

    参数:
        n: 矩阵维度
        k: 对角占优系数 (越大对角占优越强)

    返回:
        A: 系数矩阵
        b: 右侧向量
        x_true: 真实解 (随机生成)
    """
    # 主对角线元素
    main_diag = (2 * k + 1) * np.ones(n)
    # 上下对角线元素
    off_diag = -k * np.ones(n - 1)
    # 构造三对角矩阵
    A = diags([off_diag, main_diag, off_diag], [-1, 0, 1]).toarray()

    # 生成随机真实解
    x_true = np.random.rand(n)
    # 计算对应的右侧向量
    b = A @ x_true

    return A, b, x_true


# ========== 可视化分析 ==========
def plot_convergence(residuals, omega):
    """绘制单个松弛因子的收敛曲线"""
    plt.figure(figsize=(10, 6))
    plt.semilogy(residuals, 'o-', label=f'ω={omega:.2f}')
    plt.title('SOR方法收敛曲线', fontsize=14)
    plt.xlabel('迭代次数', fontsize=12)
    plt.ylabel('残差 (log scale)', fontsize=12)
    plt.legend()
    plt.grid(True, which="both", ls="--")
    plt.show()


def plot_omega_comparison(omega_list, all_residuals):
    """绘制不同松弛因子的收敛对比图"""
    plt.figure(figsize=(12, 6))
    for omega, residuals in zip(omega_list, all_residuals):
        plt.semilogy(residuals, 'o-', markersize=3, label=f'ω={omega:.2f}')

    plt.title('不同松弛因子下的SOR收敛对比', fontsize=14)
    plt.xlabel('迭代次数', fontsize=12)
    plt.ylabel('残差 (log scale)', fontsize=12)
    plt.legend()
    plt.grid(True, which="both", ls="--")
    plt.show()


def plot_omega_performance(omega_list, iterations_list):
    """绘制松弛因子与迭代次数的关系图"""
    plt.figure(figsize=(12, 6))
    plt.plot(omega_list, iterations_list, 'bo-')
    plt.title('松弛因子对收敛速度的影响', fontsize=14)
    plt.xlabel('松弛因子 ω', fontsize=12)
    plt.ylabel('达到收敛所需迭代次数', fontsize=12)
    plt.grid(True)
    plt.show()


# ========== 主程序 ==========
def main():
    # 生成测试案例
    n = 100
    A, b, x_true = generate_test_case(n)
    x0 = np.zeros(n)  # 初始猜测

    # 测试不同松弛因子
    omega_list = np.linspace(1.0, 1.9, 10)  # 测试1.0到1.9之间的5个ω值
    all_residuals = []
    iterations_list = []

    print("开始不同松弛因子的对比测试...")
    for omega in omega_list:
        # 使用SOR求解
        x_sor, residuals, iterations = sor_solver(A, b, omega, x0)
        all_residuals.append(residuals)
        iterations_list.append(iterations)

        # 计算误差
        error = np.linalg.norm(x_sor - x_true)
        print(f"ω={omega:.2f}: {iterations}次迭代, 最终误差={error:.2e}")

    # 绘制对比图
    plot_omega_comparison(omega_list, all_residuals)
    plot_omega_performance(omega_list, iterations_list)

    # 单独展示最优ω的结果
    best_idx = np.argmin(iterations_list)
    best_omega = omega_list[best_idx]
    print(f"\n最佳松弛因子: ω={best_omega:.2f}")
    plot_convergence(all_residuals[best_idx], best_omega)


if __name__ == "__main__":
    main()

SOR三图合一

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from matplotlib.font_manager import FontProperties, findfont
from scipy.sparse import diags


# ========== 自动检测中文字体 ==========
def setup_chinese_font():
    """自动检测并设置可用的中文字体"""
    chinese_fonts = ["Microsoft YaHei", "FangSong", "KaiTi", "SimHei", "SimSun"]

    found_font = None
    for font in chinese_fonts:
        try:
            if findfont(FontProperties(family=[font])) not in ['', None]:
                found_font = font
                break
        except:
            continue

    if found_font:
        plt.rcParams['font.sans-serif'] = [found_font]
        plt.rcParams['axes.unicode_minus'] = False  # 解决负号显示问题
        print(f"✅ 使用中文字体: {found_font}")
    else:
        print("⚠️ 未找到可用中文字体，请手动安装 SimHei.ttf 或指定字体路径。")


# 初始化中文字体
setup_chinese_font()


# ========== SOR迭代求解器 ==========
def sor_solver(A, b, omega, x0, tol=1e-5, max_iterations=200):
    """
    使用SOR(逐次超松弛)方法求解线性方程组 Ax = b

    参数:
        A: 系数矩阵 (n x n)
        b: 右侧向量 (n x 1)
        omega: 松弛因子 (0 < omega < 2)
        x0: 初始猜测解 (n x 1)
        tol: 收敛容差
        max_iterations: 最大迭代次数

    返回:
        x: 解向量
        residuals: 残差历史记录
        iterations: 实际迭代次数
    """
    n = len(b)
    x = x0.copy()
    residuals = []

    for k in range(max_iterations):
        x_new = x.copy()
        for i in range(n):
            sigma = 0.0
            # 前向替换 (已经更新的分量)
            for j in range(i):
                sigma += A[i, j] * x_new[j]
            # 后向替换 (尚未更新的分量)
            for j in range(i + 1, n):
                sigma += A[i, j] * x[j]

            # SOR更新公式
            x_new[i] = (1 - omega) * x[i] + (omega / A[i, i]) * (b[i] - sigma)

        # 计算残差
        residual = np.linalg.norm(A @ x_new - b)
        residuals.append(residual)

        # 检查收敛
        if residual < tol:
            break

        x = x_new

    return x_new, np.array(residuals), k + 1


# ========== 测试案例生成 ==========
def generate_test_case(n=100, k=3):
    """
    生成对角占优的三对角矩阵测试案例

    参数:
        n: 矩阵维度
        k: 对角占优系数 (越大对角占优越强)

    返回:
        A: 系数矩阵
        b: 右侧向量
        x_true: 真实解 (随机生成)
    """
    # 主对角线元素
    main_diag = (2 * k + 1) * np.ones(n)
    # 上下对角线元素
    off_diag = -k * np.ones(n - 1)
    # 构造三对角矩阵
    A = diags([off_diag, main_diag, off_diag], [-1, 0, 1]).toarray()

    # 生成随机真实解
    x_true = np.random.rand(n)
    # 计算对应的右侧向量
    b = A @ x_true

    return A, b, x_true


# ========== 可视化分析 ==========
def plot_combined_analysis(omega_list, all_residuals, iterations_list):
    """
    将三种可视化结果整合到一张图中：
        1. 不同 ω 的收敛曲线
        2. 所有 ω 的残差对比
        3. 松弛因子与迭代次数关系图
    """
    # 准备数据：收敛曲线数据
    df_curve = pd.DataFrame()
    for residuals, omega in zip(all_residuals, omega_list):
        df_single = pd.DataFrame({
            'Iteration': range(1, len(residuals) + 1),
            'Residual': residuals,
            'Omega': f'ω={omega:.2f}'
        })
        df_curve = pd.concat([df_curve, df_single], ignore_index=True)

    # 创建画布：3个子图横向排列
    fig, axes = plt.subplots(1, 3, figsize=(20, 6))

    # 子图1: 多个 ω 的收敛曲线对比
    sns.lineplot(data=df_curve, x='Iteration', y='Residual', hue='Omega', style='Omega',
                 markers=True, dashes=False, ax=axes[0])
    axes[0].set_title('不同松弛因子下的收敛曲线')
    axes[0].set_yscale('log')
    axes[0].set_xlabel('迭代次数')
    axes[0].set_ylabel('残差 (log scale)')
    axes[0].grid(True, which="both", linestyle="--")
    axes[0].legend(title='松弛因子 ω')

    # 子图2: 松弛因子 vs 迭代次数关系图
    axes[1].plot(omega_list, iterations_list, 'bo-')
    axes[1].set_title('松弛因子对收敛速度的影响')
    axes[1].set_xlabel('松弛因子 ω')
    axes[1].set_ylabel('达到收敛所需迭代次数')
    axes[1].grid(True)
    axes[1].axvline(x=omega_list[np.argmin(iterations_list)], color='r', linestyle='--', linewidth=1.5)

    # 子图3: 最佳 ω 的收敛曲线放大展示
    best_idx = np.argmin(iterations_list)
    best_omega = omega_list[best_idx]
    best_residuals = all_residuals[best_idx]

    axes[2].semilogy(best_residuals, 'o-', label=f'ω={best_omega:.2f}')
    axes[2].set_title(f'最佳松弛因子收敛曲线 (ω={best_omega:.2f})')
    axes[2].set_xlabel('迭代次数')
    axes[2].set_ylabel('残差 (log scale)')
    axes[2].legend()
    axes[2].grid(True, which="both", ls="--")

    plt.tight_layout()
    plt.show()


# ========== 主程序 ==========
def main():
    # 生成测试案例
    n = 100
    A, b, x_true = generate_test_case(n)
    x0 = np.zeros(n)  # 初始猜测

    # 测试不同松弛因子
    omega_list = np.linspace(1.0, 1.9, 10)  # 测试1.0到1.9之间的5个ω值
    all_residuals = []
    iterations_list = []

    print("开始不同松弛因子的对比测试...")
    for omega in omega_list:
        x_sor, residuals, iterations = sor_solver(A, b, omega, x0)
        all_residuals.append(residuals)
        iterations_list.append(iterations)

        error = np.linalg.norm(x_sor - x_true)
        print(f"ω={omega:.2f}: {iterations}次迭代, 最终误差={error:.2e}")

    # 替换原有多个绘图函数，统一显示三合一图像
    plot_combined_analysis(omega_list, all_residuals, iterations_list)

    # 单独展示最优ω的结果
    best_idx = np.argmin(iterations_list)
    best_omega = omega_list[best_idx]
    print(f"\n最佳松弛因子: ω={best_omega:.2f}")


if __name__ == "__main__":
    main()