解线性方程组python实现迭代法(Jacobi迭代、Gauss-Seidel迭代、松弛迭代)

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)

  • 3
    点赞
  • 42
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

安心不心安

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值