雅可比迭代和高斯—赛德尔迭代法

本文介绍了线性方程组的两种迭代解法——雅可比迭代法和高斯-赛德尔迭代法,分别提供了MATLAB和Python的代码示例。这两种方法通过对系数矩阵进行分解,然后进行迭代更新,以求解线性方程组。文章还给出了具体的测试用例和求解过程。
摘要由CSDN通过智能技术生成

一、雅可比迭代法

       对于线性方程组AX=b,我们首先将系数矩阵A分解为对角矩阵D、下三角矩阵L和上三角矩阵U:

\left ( D-L-U \right )X= b

DX=\left ( L+U \right )X+b

X^{k+1}=D^{-1}\left ( L+U \right )X^{k}+D^{-1}b

1.1雅可比迭代法的matlab代码

  在这里,我们求解下面的带状方程(以下程序均是以求解该带状方程为例):

12x_{1}-2x_{2}+x_{3} =5x

-2x_{1}+12x_{2}-2x_{3}+x_{4}=5

x_{1}-2x_{2}+12x_{3}-2x_{4}+x_{5}=5

x_{2}-2x_{3}+12x_{4}-2x_{5}+x_{6}=5

.............

x_{46}-2x_{47}+12x_{48}-2x_{49}+x_{50}=5 

x_{47}-2x_{48}+12x_{49}-2x_{50}=5 

x_{48}-2x_{49}+12x_{50}=5

function X0=jacobi(A,b,X0,delta,max1)
%输入  -A代表线性方程组AX=b的系数矩阵
%      -b代表线性方程组AX=b右侧的数值
%      -X0代表线性方程组AX=b进行高斯-赛德尔迭代法求解的迭代初值
%      -delta代表余项AX(k)-B的范数允许误差
%      -max1代表迭代的次数
%输出  -X0代表通过雅可比迭代法求解线性方程组AX=b的解
[N,N]=size(A);
L=-tril(A,-1);
U=-triu(A,1);
D=A+L+U;
B=inv(D)*(L+U)
f=inv(D)*b;
for k=1:max1
    X0=B*X0+f;
    err=abs(norm(A(:,:)*X0(:)-b(:),2))
    if err<delta
        break
    end
end   

 1.2雅可比迭代法的python代码

import numpy as np

def jacobi(A,b,X0,max1):
    '''A代表线性方程组AX=b的系数矩阵
           b代表线性方程组AX=b右边的部分
           X0代表高斯—赛德尔迭代的初始值
           max1代表迭代的次数'''
    n = np.shape(A)[0]
    L = -np.tril(A, -1)
    U = -np.triu(A, 1)
    D = A + L + U
    B = np.dot(np.linalg.inv(D),L+U)
    f = np.dot(np.linalg.inv(D),b)
    for i in range(max1):
        X0 = np.dot(B, X0) + f
    return X0

n=50
#线性方程组AX=b右边的部分
b=np.zeros((n,1))
for i in range(n):
    b[i]=5
#线性方程组AX=b的系数矩阵
A=np.zeros((n,n))
for i in range(n):
    A[i,i]=12
    if i>=0 and i<=48:
        A[i,i+1]=-2
    if i>=0 and i<=47:
        A[i,i+2]=1
    if i>=1:
        A[i,i-1]=-2
    if i>=2:
        A[i,i-2]=1
#迭代的初始值
X0=np.zeros((n,1))
for i in range(n):
    X0[i]=0
#进行迭代的次数
max1=50
#进行高斯—赛德尔迭代求解线性方程组AX=b的解
X=jacobi(A,b,X0,max1)
#输出由高斯—赛德尔迭代求得的线性方程组的解
print(X)

二、高斯—赛德尔迭代法

      高斯—赛德尔迭代法是再雅可比迭代法的基础上,在计算y_{k+1}时尽可能地用最新的x_{k+1}来替换x_{k} 。同样对于线性方程组 AX=b,我们首先将系数矩阵A分解为对角矩阵D、下三角矩阵L和上三角矩阵U:

\left ( D-L-U \right )X= b

DX^{\left ( k+1 \right )}=LX^{\left ( k+1 \right )}+UX^{\left ( k \right )}+b

\left ( D-L \right )X^{\left ( k+1 \right )}=UX^{\left ( k \right )}+b

X^{\left ( k+1 \right )}=\left ( D-L^{} \right )^{-1}UX^{\left ( k \right )}+\left ( D-L^{} \right )^{-1}b 

 2.1高斯—赛德尔迭代法的matlab代码

function X0=Gauss_Saidel(A,b,X0,max1)
%输入  -A代表线性方程组AX=b的系数矩阵
%      -b代表线性方程组AX=b右侧的数值
%      -X0代表线性方程组AX=b进行高斯-赛德尔迭代法求解的迭代初值
%      -max1代表迭代的次数
%输出  -X0代表通过高斯-赛德尔迭代法求解线性方程组AX=b的解
[N,N]=size(A);
L=-tril(A,-1);
U=-triu(A,1);
D=A+L+U;
B=inv(D-L)*U;
f=inv(D-L)*b;
for k=1:max1
    X0=B*X0+f
end

2.2高斯—赛德尔迭代法的python代码

import numpy as np

def Gauss_Saidel(A,b,X0,max1):
    '''A代表线性方程组AX=b的系数矩阵
       b代表线性方程组AX=b右边的部分
       X0代表高斯—赛德尔迭代的初始值
       max1代表迭代的次数'''
    n=np.shape(A)[0]
    L=-np.tril(A,-1)
    U=-np.triu(A,1)
    D=A+L+U
    B=np.dot(np.linalg.inv(D-L),U)
    f=np.dot(np.linalg.inv(D-L),b)
    for i in range(max1):
        X0=np.dot(B,X0)+f
    return X0

n=50
#线性方程组AX=b右边的部分
b=np.zeros((n,1))
for i in range(n):
    b[i]=5
#线性方程组AX=b的系数矩阵
A=np.zeros((n,n))
for i in range(n):
    A[i,i]=12
    if i>=0 and i<=48:
        A[i,i+1]=-2
    if i>=0 and i<=47:
        A[i,i+2]=1
    if i>=1:
        A[i,i-1]=-2
    if i>=2:
        A[i,i-2]=1
#迭代的初始值
X0=np.zeros((n,1))
for i in range(n):
    X0[i]=0
#进行迭代的次数
max1=50
#进行高斯—赛德尔迭代求解线性方程组AX=b的解
X=Gauss_Saidel(A,b,X0,max1)
#输出由高斯—赛德尔迭代求得的线性方程组的解
print(X)

雅可比迭代法高斯-德尔迭代法都是用于求线性方程组的数值方法。它们通常用于决大规模稀疏矩阵的问题。 **雅可比迭代法**(也称直接法)是一种基于方程组导数的迭代过程。假设有一个线性系统 \(Ax = b\),其中 \(A\) 是系数矩阵,\(x\) 是未知向量,\(b\) 是常数项。雅可比迭代法通过构建并求 \(J(A)x^{(k+1)} = b - Ax^{(k)}\) 来迭代逼近,其中 \(J(A)\) 是 \(A\) 的雅可比矩阵(即对 \(A\) 按元素求导后的矩阵)。每一步迭代都要求计算雅可比矩阵乘以当前猜测值,这在矩阵非常大或非对角占优时效率较低。 **高斯-德尔迭代法**(Gauss-Seidel method)也是一种迭代线性方程组的方法,但它采用分块的方式更新每个变量。算法每次只考虑部分已知的变量值来更新下一个变量,顺序通常是自左到右、自上而下。高斯-德尔通常适用于对角占优矩阵,因为它在每个步骤中使用了更精确的信息。 以下是简单的Python代码示例(仅适用于一维情况,实际应用需处理多维数组): ```python import numpy as np def jacobian_iterate(A, b, x0, tolerance=1e-6): def jacobi(A, x): return A @ x x = x0.copy() delta = float('inf') while delta > tolerance: old_x = x x = jacobi(A, (b - A @ x) / A.diagonal()) delta = np.linalg.norm(x - old_x) return x def gauss_seidel(A, b, x0, tolerance=1e-6): n = len(b) for i in range(n): if i == 0: for k in range(i + 1, n): b[i] -= A[i, k] * x0[k] else: for k in range(i): b[i] -= A[i, k] * x0[k] x0[i] = (b[i] - sum(A[i, j] * x0[j] for j in range(i))) / A[i, i] return x0 # 使用示例: A = np.array([[4, 1], [1, 3]]) b = np.array([5, 7]) x0 = np.zeros_like(b) jacobi_solution = jacobian_iterate(A, b, x0) gs_solution = gauss_seidel(A, b, x0) print("雅可比迭代法结果:", jacobi_solution) print("高斯-德尔迭代法结果:", gs_solution) ```
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

打赏作者

天下弈星~

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

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

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

打赏作者

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

抵扣说明:

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

余额充值