高斯赛德尔迭代
与雅克比迭代对比:Python——雅克比迭代求线性方程组的根
计算代码
矩阵A——方程组的系数。代码中只需改变矩阵系数即可
矩阵B——方程组等号右侧的常数。
#高斯赛德尔迭代
#求解方程组:
#10*x1 - 2*x2 - 1*x3 = 3
#-2*x1 + 10*x2 - 1*x3 = 15
#-1*x1 - 2*x2 + 5*x3 = 10
#AX = B
#只用改变矩阵A、B的系数即可
A = [[10, -2, -1], [-2, 10, -1], [-1, -2, 5]]
B = [[3], [15], [10]]
#初始值
x = [0, 0, 0]
aij = None
n = 0
while True:
v_max = 0
n += 1
print(f'第{n}次迭代')
for i in range(len(A)):
a = A[i]
ax = 0
for j in range(len(a)):
if j == i :
aij = a[j]
else:
ax += a[j]*x[j]
x_new = (B[i][0] - ax) / aij
v = abs(x[i] - x_new)
v_max = max(v_max, v)
x[i] = x_new
if v_max <= 10**-3:
print(x)
print('e=',v_max)
break
else:
print(x)
print('e=',v_max)
结果:
第1次迭代
[0.3, 1.56, 2.684]
e= 2.684
第2次迭代
[0.8804000000000001, 1.94448, 2.9538719999999996]
e= 0.5804
第3次迭代
[0.9842831999999999, 1.9922438399999998, 2.993754176]
e= 0.10388319999999984
第4次迭代
[0.9978241856000001, 1.9989402547200001, 2.999140939008]
e= 0.013540985600000166
第5次迭代
[0.9997021448448, 1.9998545228697602, 2.999882238116864]
e= 0.0018779592447999427
第6次迭代
[0.9999591283856384, 1.9999800494888142, 2.9999838454726535]
e= 0.0002569835408383625
高斯赛德尔代收敛的充分必要条件
L——矩阵A的下三角矩阵
U——矩阵A的上三角矩阵
D——矩阵A的对角矩阵
例如:
from numpy import *
# numpy.linalg模块包含线性代数的函数。使用这个模块,可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等
from numpy.linalg import *
A = mat([[25, 5, 1], [64, 8, 1], [144, 12, 1]])
B = mat([[106.8], [177.2], [279.2]])
L = tril(A, -1) #下三角矩阵
U = triu(A,1) #上三角矩阵
D = diag(diag(A)) #对家矩阵
G = dot(inv(D+L),-U)
print(eigvals(G)) #特征值
结果:
[0. 0.8668019 4.1531981]
特征值中只要有一个绝对值>1,用高斯赛德尔迭代求解就不收敛
所以上面的方程组用高斯赛德尔迭代不收敛
再例如:
from numpy import *
# numpy.linalg模块包含线性代数的函数。使用这个模块,可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等
from numpy.linalg import *
A = mat([[2, 1, 1], [1, 2, 1], [1, 1, 2]])
B = mat([[4], [2], [0]])
L = tril(A, -1) #下三角矩阵
U = triu(A,1) #上三角矩阵
D = diag(diag(A)) #对角矩阵
G = dot(inv(D+L),-U)
print(eigvals(G)) #特征值
结果:
[0. +0.j 0.3125+0.16535946j 0.3125-0.16535946j]
三个特征值的绝对值都<1,所以用高斯赛德尔迭代可以求解这个方程组。
先判断是否收敛再迭代的综合python代码
#高斯赛德尔迭代
#求解方程组:
#10*x1 - 2*x2 - 1*x3 = 3
#-2*x1 + 10*x2 - 1*x3 = 15
#-1*x1 - 2*x2 + 5*x3 = 10
#AX = B
from numpy import *
# numpy.linalg模块包含线性代数的函数。使用这个模块,可以计算逆矩阵、求特征值、解线性方程组以及求解行列式等
from numpy.linalg import *
#只用改变矩阵A、B的系数即可
A = [[10, -2, -1], [-2, 10, -1], [-1, -2, 5]]
B = [[3], [15], [10]]
#判断该矩阵用高斯-赛德尔迭代是否收敛
L = tril(A, -1) #下三角矩阵
U = triu(A,1) #上三角矩阵
D = diag(diag(A)) #对角矩阵
G = dot(inv(D+L),-U)
root = eigvals(G)
flag = True
for r in root:
if r > 1:
print('高斯-赛德尔迭代不收敛')
flag = False
break
else:
pass
#如果迭代收敛,求迭代
if flag:
# 初始值
x = [0, 0, 0]
aij = None
n = 0
while True:
v_max = 0
n += 1
print(f'第{n}次迭代')
for i in range(len(A)):
a = A[i]
ax = 0
for j in range(len(a)):
if j == i:
aij = a[j]
else:
ax += a[j] * x[j]
x_new = (B[i][0] - ax) / aij
v = abs(x[i] - x_new)
v_max = max(v_max, v)
x[i] = x_new
if v_max <= 10 ** -3:
print(x)
print('e=', v_max)
break
else:
print(x)
print('e=', v_max)
结果:
第1次迭代
[0.3, 1.56, 2.684]
e= 2.684
第2次迭代
[0.8804000000000001, 1.94448, 2.9538719999999996]
e= 0.5804
第3次迭代
[0.9842831999999999, 1.9922438399999998, 2.993754176]
e= 0.10388319999999984
第4次迭代
[0.9978241856000001, 1.9989402547200001, 2.999140939008]
e= 0.013540985600000166
第5次迭代
[0.9997021448448, 1.9998545228697602, 2.999882238116864]
e= 0.0018779592447999427
第6次迭代
[0.9999591283856384, 1.9999800494888142, 2.9999838454726535]
e= 0.0002569835408383625