运用雅可比(Jacobi)和高斯-赛德尔(GS公式)求解方程组

实验名称: 运用雅克比和高斯赛德尔公式求解方程组

实验目的:

a. 比较两种方法的收敛速度
b. 验证收敛条件的正确性

实验内容

1.比较两种方法的收敛速度

Jacobi公式原理

xi(k+1)=1aii(bi−∑j=1,j!=inaijxj(k)),i=1,2,3,...nx_{i}^{(k+1)}=\frac {1} {a_{ii}} (b_i-\sum _{j=1,j!=i}^na_{ij}x_j^{(k)}), i = 1,2,3,...n xi(k+1)=aii1(bi−j=1,j!=i∑naijxj(k)),i=1,2,3,...n

GS公式原理

xi(k+1)=1aii(bi−∑j=1i−1aijxj(k+1)−∑j=i+1naijxj(k)),i=1,2,3,...nx_{i}^{(k+1)}=\frac {1} {a_{ii}} (b_i-\sum _{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum _{j=i+1}^na_{ij}x_j^{(k)}), i = 1,2,3,...n xi(k+1)=aii1(bi−j=1∑i−1aijxj(k+1)−j=i+1∑naijxj(k)),i=1,2,3,...n

选用的问题1:
A=[10−1−2−110−2−1−15]y=[7.28.34.2]TA=\begin{bmatrix} 10 & -1 & -2 \\ -1 & 10 & -2 \\ -1 & -1 & 5 \end{bmatrix} \\ y = \begin {bmatrix} 7.2 & 8.3 & 4.2 \end{bmatrix}^T A=⎣⎡10−1−1−110−1−2−25⎦⎤y=[7.28.34.2]T
A为一个行对角占优矩阵，故其必然对Jacobi公式、GS公式收敛。

选用的问题2:
A=[1001−518−3032−811−227]y=[−7112317]TA=\begin{bmatrix} 10 & 0 & 1 & -5 \\ 1 & 8 & -3 & 0 \\ 3 & 2 & -8 & 1 \\ 1 & -2 & 2 & 7 \end{bmatrix} \\ y = \begin {bmatrix} -7 & 11 & 23 & 17 \end{bmatrix}^T A=⎣⎢⎢⎡10131082−21−3−82−5017⎦⎥⎥⎤y=[−7112317]T
处理结果:

问题1

设定精度10(−15)10^{(-15)}10(−15),Jacobi迭代了32次

设定精度10(−15)10^{(-15)}10(−15),GS迭代了18次

问题2

设定精度10(−15)10^{(-15)}10(−15),Jacobi迭代了47次

设定精度10(−15)10^{(-15)}10(−15),GS迭代了28次

结论: GS的收敛速度总是优于Jacobi

2.验证收敛条件的正确性

收敛定理:

对于迭代格式 x(k+1)=Gx(k)+dx^{(k+1)}=Gx^{(k)}+dx(k+1)=Gx(k)+d, 如果迭代矩阵GGG满足某一范数∣∣G∣∣<1||G||<1∣∣G∣∣<1, 那么该迭代格式对于任何初值x(0)x^{(0)}x(0)均收敛

因此，找出某一迭代矩阵的某一范数∣∣G∣∣>=1||G||>=1∣∣G∣∣>=1，该迭代格式则不收敛。

选用的矩阵:
C=[−2154−814−11]C=\begin{bmatrix} -2 & 1 & 5 \\ 4 & -8 & 1 \\ 4 & -1 & 1 \\ \end{bmatrix} C=⎣⎡−2441−8−1511⎦⎤
Jacobi迭代为例:
G=−D−1(L+U)G=-D^{-1}(L+U) G=−D−1(L+U)
计算得:
G=[0−0.5−2.5−0.50−0.1254−10]G=\begin{bmatrix} 0 & -0.5 & -2.5 \\ -0.5 & 0 & -0.125 \\ 4 & -1 & 0 \\ \end{bmatrix} G=⎣⎡0−0.54−0.50−1−2.5−0.1250⎦⎤
计算G的无穷范数得:
∣∣G∣∣=5.0||G||=5.0 ∣∣G∣∣=5.0
最终结果为经过10000次迭代都不能满足精度要求，发散。

心得体会

1.从收敛速度的角度出发，GS总是优于Jacobi，因为它用到了使用已经计算的“新值”替代“旧值"的思想。

2.在计算效率方面，Jacobi非常适合并行计算，因为每次迭代中的变量都不会改变，直到迭代完成。相反，只要有新的更新可用，GS就可以替换每个变量。

3.总之，Jacobi逐渐靠近正确解；GS需要消耗大量内存以达到更快的收敛速度。

附录

Jacobi、GS的Python实现

def Jacobi(A,y,n=100,c=1e-15):if len(A) == len(y):   #保证维度正确all_x = []#迭代初值 初始化为单行全0矩阵x = [] for i in range(len(y)):x.append([0])count = 0 #迭代次数计数while count < n:nx = [] #保存单次迭代后的值的集合for i in range(len(x)):nxi = y[i][0]for j in range(len(A[i])):if j!=i:nxi = nxi+(-A[i][j])*x[j][0]nxi = nxi / A[i][i]nx.append([nxi]) #迭代计算得到的下一个xi值lc = [] #存储两次迭代结果之间的误差的集合for i in range(len(x)):lc.append(abs(x[i][0]-nx[i][0]))if max(lc) < c:print("迭代了{}次".format(count))print("结果的精度是{}".format(c))return nx,all_x #当误差满足要求时 返回计算结果x = nxall_x.append(x)count = count + 1print("经过{}次迭代未达到精度要求!".format(count))return False #若达到设定的迭代结果仍不满足精度要求 则方程无解else:return False

def GS(A,y,n=100,c=1e-15):if len(A) == len(y):   #保证维度正确all_x = []#迭代初值 初始化为单行全0矩阵x = [] for i in range(len(y)):x.append([0])count = 0 #迭代次数计数while count < n:nx = [] #保存单次迭代后的值的集合for i in range(len(x)):nxi = y[i][0]for j in range(len(A[i])):if j < i:                           #j:1~i-1用新值nxi = nxi - A[i][j] * nx[j][0]elif j != i and j > i:             #j:i+1~n用旧值nxi = nxi - A[i][j] * x[j][0]nxi = nxi / A[i][i]nx.append([nxi]) #迭代计算得到的下一个xi值lc = [] #存储两次迭代结果之间的误差的集合for i in range(len(x)):lc.append(abs(x[i][0]-nx[i][0]))if max(lc) < c:print("迭代了{}次".format(count))print("结果的精度是{}".format(c))return nx,all_x #当误差满足要求时 返回计算结果x = nxall_x.append(x)count = count + 1print("经过{}次迭代未达到精度要求!".format(count))return False #若达到设定的迭代结果仍不满足精度要求 则方程无解else:return False

求矩阵的无穷范数

#矩阵的无穷范数
norm = 0
for i in range(len(G)):res = 0for j in range(len(G[i])):res += abs(G[i][j])if res > norm:norm = res
norm