1. Jacobi迭代

        Jacobi迭代法是一种用于求解线性方程组的迭代算法。它属于迭代法中的直接迭代法，通过不断迭代更新解向量来逼近线性方程组的解。

        Jacobi迭代法的基本概念如下：

1. 给定线性方程组的系数矩阵A和右侧常数向量b。

2. 将系数矩阵A进行对角分解，得到三个矩阵D、L和U，其中D是A的对角矩阵，L是A的严格下三角矩阵（即主对角线以下元素为0），U是A的严格上三角矩阵（即主对角线以上元素为0）。

3. 初始化解向量x为任意初始值（可以是全0向量）。

4. 迭代更新解向量x，直到满足收敛条件：

   • 计算残差向量r = b - Ax，其中A是系数矩阵，x是当前解向量。
   • 使用更新公式：x = D^(-1) * (b - (L + U)x)，其中D^(-1)表示D的逆矩阵。
   • 重复进行以上两步操作，直到残差向量r的范数小于设定的收敛阈值或达到最大迭代次数。

5. 返回近似解向量x。

        Jacobi迭代法的收敛性和速度与系数矩阵A的特性有关。当系数矩阵A是对角占优或严格对角占优时，Jacobi迭代法具有收敛性。然而，在某些情况下，Jacobi迭代法的收敛速度较慢，因此可以结合其他迭代方法，如Gauss-Seidel迭代法或逐次超松弛法，来改进迭代效果。

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

        Jacobi 雅可比迭代法
"""
import numpy as np


def Jacobi(A, b, max_iter=100, epsilon=1e-10):
    n = len(b)
    x = np.zeros(n)  # 初始解
    x_new = np.zeros(n)  # 存储更新后的解向量

    for iter in range(max_iter):
        for i in range(n):
            sum_ax = np.dot(A[i], x) - A[i][i] * x[i]  # 计算除去对角元素的和
            x_new[i] = (b[i] - sum_ax) / A[i][i]  # 更新解向量的第i个分量

        # 判断终止条件
        if np.linalg.norm(x_new - x) < epsilon:
            break

        x = np.copy(x_new)  # 更新解向量
    return x


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

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

x = Jacobi(A, b)
print("解：", x)

2. Gauss-Seidel迭代

        Gauss-Seidel迭代法是一种求解线性方程组的迭代算法，它是Jacobi迭代法的改进版。该方法通过不断迭代更新解向量的每个分量来逼近线性方程组的解。

        Gauss-Seidel迭代法的基本概念如下：

1. 给定线性方程组的系数矩阵A和右侧常数向量b。

2. 将系数矩阵A进行对角分解，得到三个矩阵D、L和U，其中D是A的对角矩阵，L是A的严格下三角矩阵（即主对角线以下元素为0），U是A的严格上三角矩阵（即主对角线以上元素为0）。

3. 初始化解向量x为任意初始值（可以是全0向量）。

4. 迭代更新解向量x，直到满足收敛条件：

   • 对于第i个未知数，计算其新值：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))，其中a(i, j)表示系数矩阵A的第i行第j列元素。
   • 通过不断更新解向量x中的每个分量，得到新的解向量。
   • 重复进行以上两步操作，直到满足收敛条件：残差向量的范数小于设定的收敛阈值或达到最大迭代次数。

5. 返回近似解向量x。

        Gauss-Seidel迭代法相比于Jacobi迭代法的改进之处在于，在每次迭代中，它使用了已经更新过的解向量的新分量来计算下一个未知数的新值，从而加快了收敛速度。然而，与Jacobi迭代法类似，Gauss-Seidel迭代法也需要满足一定的收敛条件才能得到有效的结果。

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

        Gauss-Seidel 高斯赛德尔迭代
"""
import numpy as np


def Jacobi(A, b, max_iter=100, epsilon=1e-10):
    n = len(b)
    x = np.zeros(n)  # 初始解
    x_new = np.zeros(n)  # 存储更新后的解向量

    for iter in range(max_iter):
        for i in range(n):
            sum_ax = np.dot(A[i], x_new) - A[i][i] * x[i]  # 计算除去对角元素的和,注意：高斯-赛德尔，这里使用x_new而不是x
            x_new[i] = (b[i] - sum_ax) / A[i][i]  # 更新解向量的第i个分量

        # 判断终止条件
        if np.linalg.norm(x_new - x) < epsilon:
            break

        x = np.copy(x_new)  # 更新解向量
    return x


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

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

x = Jacobi(A, b)
print("解：", x)

3. 松弛迭代

        松弛迭代（Relaxation Iteration）是一种用于求解线性方程组的迭代算法，常用于解决稀疏矩阵和大规模线性方程组的问题。它在Gauss-Seidel迭代法的基础上引入了一个松弛因子，用于加速收敛过程。

        松弛迭代的基本概念如下：

1. 给定线性方程组的系数矩阵A和右侧常数向量b。

2. 将系数矩阵A进行对角分解，得到三个矩阵D、L和U，其中D是A的对角矩阵，L是A的严格下三角矩阵（即主对角线以下元素为0），U是A的严格上三角矩阵（即主对角线以上元素为0）。

3. 初始化解向量x为任意初始值（可以是全0向量）。

4. 迭代更新解向量x，直到满足收敛条件：

   • 对于第i个未知数，计算其新值：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))，其中a(i, j)表示系数矩阵A的第i行第j列元素，w是松弛因子（取值范围为0 < w < 2）。
   • 通过不断更新解向量x中的每个分量，得到新的解向量。
   • 重复进行以上两步操作，直到满足收敛条件：残差向量的范数小于设定的收敛阈值或达到最大迭代次数。

5. 返回近似解向量x。

        松弛迭代法引入了松弛因子w，用于调整每次迭代更新的幅度。当松弛因子接近于1时，算法的收敛速度较慢；当松弛因子接近于0或2时，算法的收敛速度较快。恰当选择松弛因子可以显著提高算法的收敛速度，但过大或过小的松弛因子可能会导致算法发散。

        需要注意的是，在松弛迭代法中，松弛因子的选择是一个重要的问题，需要根据具体问题进行调试和优化。通常情况下，可以通过试验和经验来确定合适的松弛因子。

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

        松弛迭代
"""
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


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

b = np.array([7, 7, 8])
x0 = np.zeros(3)  # 初始解向量
omega = 1.2  # 松弛因子
x = relaxation_iteration(A, b, x0, omega=omega)

print("解：", x)