我正在写一个代码,用高斯-蔡德尔方法解决方程组。结果不正确,但我不知道为什么。可能是什么问题?
import numpy as np
A = np.matrix([[10, 3, 0],
[3, 15, 1],
[0, 1, 7]]).astype(float)
f = np.matrix([2, 12, 5]).astype(float)
f = np.reshape(f, (3, 1))
x_new = np.matrix([0, 0, 0]).astype(float)
x_new = np.reshape(x_new, (3, 1))
N = 100
D = np.matrix(np.zeros(A.shape))
for n in range(A.shape[0]):
D[n, n] = A[n, n]
A1 = np.tril(A, -1)
A2 = np.triu(A, 1)
D_inv = np.linalg.inv(D)
x = np.matrix([0, 0, 0]).astype(float)
x = np.reshape(x, (3, 1))
eps = 1e-5
dim = len(A)
def gauss_zeidel(x):
i = 0
t = 0
s1 = 0
s2 = 0
while t <= N:
x_new = np.copy(x)
for i in range(dim):
print("--new cycle--")
print("===============================")
for j in range (i-1):
s1 += A[i,j] * x_new.item(j)
for j in range(i+1, dim):
s2 += A[i,j] * x.item(j)
print("iteration result")
x_new.itemset(i, (f.item(i) - s1 - s2)/A[i,i])
print((x.item(i)), i)
print(t)
print("===============================")
x = x_new
t+=1
print(x)
1条答案
按热度按时间cgh8pdjw1#
两个错误(除了非常低效之外,因为numpy的要点是永远不要迭代数组。但我认为这是一个教学目的。我很确定numpy或scipy系统求解已经在合适的时候使用Gauss-Seidel,除非他们使用更好的方法;但更快)
1.仅将s1和s2设置为0一次。显然,这两个for循环应该计算Σ ¹ a.x ¹和Σ a x。所以s1和s2应该在这些循环之前设置为0。换句话说,你计算的是这个总和加上前面的结果。所以呢
1.我刚刚使用了和你一样的基于1的计数(从1到n的总和),这是通常在数学和大多数高斯-赛德尔算法描述中使用的。但是numpy是基于0的索引。
range
排除上限。与初学者的想法相反,哪一个更容易阅读:这意味着range(a,b)
正好迭代b-a
项。range(a,c)
然后range(c,b)
将从a到B的所有项迭代一次且仅一次(包括c,它在第二个范围内迭代,而不是第一个)。这意味着当你的s1是从range(i-1)
计算出来的(隐含的range(0,i-1)
和s2是从range(i+1,dim)
计算出来的),s1和s2都是由i-1+n-(i+1) = n-2
迭代组成的。当你知道它应该是n-1
(除了对角元素之外的所有元素)时。因为i-1
和i+1
之间有2个元素,所以你不会只跳过s1和s2之间的1个对角元素。所以,长话短说,第一个
range
应该是range(i)
(意思是“从0包括到i排除”,或者“从0到i-1”)通过这2个修改(s1=0,s2=0;范围(i)而不是范围(i-1)),结果正确。
当然,有许多可能的优化。
只有几个