在python中实现Gauss-Jacobi算法时,我发现两种不同的实现需要显著不同的迭代次数才能收敛
第一个实现是我最初提出的
import numpy as np
def GaussJacobi(A, b, x, x_solution, tol):
k = 0
N = A.shape[0]
D = np.diag(A)
R = A-np.diagflat(D);
while(checkTol(tol, x, x_solution)):
x_new = np.zeros(N, dtype=np.double) #x(k+1)
for i in range(N):
aii = D[i]
bi = b[i]
s = np.dot(R[i], x)
x_n[i] = (1/aii)*(bi - s)
x = x_new
k+=1
print('x(%d) =' % k, x)
return k
第二个实现基于this article.
def GaussJacobi(A, b, x, x_solution, tol):
k = 0
N = A.shape[0]
D = np.diag(A)
R = A-np.diagflat(D);
while(checkTol(tol, x, x_solution)):
for i in range(N):
x = (b - np.dot(R, x)) / D
k+=1
print('x(%d) =' % k, x)
return k
当解决以下问题时
A = [ 4, -1, 0, -1, 0, 0]
[-1, 4, -1, 0, -1, 0]
[ 0, -1, 4, 0, 0, -1]
[-1, 0, 0, 4, -1, 0]
[0, -1, 0, -1, 4, -1]
[0, 0, -1, 0, -1, 4]
b = [2, 1, 2, 2, 1, 2]
x_solution =[1, 1, 1, 1, 1, 1]
x0 = [0, 0, 0, 0, 0, 0]
第一个实现需要37次迭代才能收敛,误差为1e-8,而第二个实现只需要7次迭代才能收敛
是什么使得第二个实现比第一个实现快得多
编辑:
我实现了另外两种方法,Gauss-Seidel方法和SOR方法。这两种方法的实现方式与我最初的慢高斯-雅可比方法类似
我对100个NxN对角占优矩阵进行了随机测试,每个N=4…20,以获得收敛前的平均迭代次数
N Gauss-Jacobi Gauss-Jacobi Fast Gauss Seidel SOR -- w=1.5
--- -------------- ------------------- -------------- --------------
4 40.96 17.04 40.6804 40.9204
5 49.11 17.25 48.7489 48.9389
6 56.11 16.04 55.6789 55.9089
7 70.26 18 69.6774 70.0074
8 76.4 16.54 75.756 76.236
9 83.56 17.03 82.8344 83.1044
10 92.33 16.24 91.5267 91.7267
11 98.02 16.59 97.1598 97.4598
12 107.39 15.98 106.436 106.756
13 123.48 17.75 122.375 122.655
14 125.07 16.04 123.949 124.239
15 132.41 16.68 131.206 131.496
16 145 16.31 143.67 143.91
17 149.66 16.75 148.283 148.493
18 154.21 15.58 152.788 153.078
19 163.18 16.51 161.668 161.918
20 167.58 15.38 166.014 166.254
更快的Gauss-Jacobi实现不仅比其他任何实现都要快,而且似乎不像其他方法那样随着数组大小的增加而增加
当检查正在运行的方法时,快速方法似乎在其第一次迭代中创建了一个非常好的猜测
我的猜测是,它必须与np.dot
函数有关,但我不明白为什么这与独立处理每个点积不同
您的第二个实现对
N
每增加k
个实际的迭代,因为对x
的赋值已经覆盖了整个向量。因此,它的“优势”随着问题的规模而增加相关问题 更多 >
编程相关推荐