线性联立方程的雅可比迭代解(jacobi Iteration).(python,数值积分)

第八课 线性联立方程的雅可比迭代

迭代方法

在前面几节中，已经描述了线性代数方程组的“直接”解法。我们所说的“直接”是指求解方法是通过固定次数的运算得到答案。由于存在舍入误差，解决方案在计算机硬件允许的范围内是“精确的”。用户没有机会在解决过程中进行干预。相比之下，本节讨论的是迭代或者是“间接”求解方法。这些都是需要用户首先猜测一个答案，然后一次又一次地进行修正来求解的。下面将描述几种方法，这些方法在校正的技术上有所不同，并且会发现具有对真实解的不同收敛速度。是否会实现收敛的问题将不在我的分享的讨论范围之内。一般来说，具有“对角优势矩阵”的方程很可能有收敛的迭代解。这个明聪的粗略定义为，任何一行中的对角线项都大于非对角线项的总和。正如我们之前看到的，在许多工程应用中，这种对角优势是广泛存在的。
由于用户的初始猜测不太可能是准确的，所以解决方案通常会逐渐变得越来越接近真正的值，这给用户提供了再次选择的机会，因为选择“足够接近”期望的值的尝试是可能的。例如，如果解决方案在某种意义上与之前的迭代相比几乎没有变化，则可以终止该过程。因此，一般而言，在这些方法种，事先是不知道一个方程的解的。

迭代过程

假设我们刚开始去尝试一个解决方程的迭代办法，对下面这个方程：

迭代方法的主要特点是要把未知数放在方程的两边，所以上面方程可以变化为：

这种最简单的方法叫做‘雅可比迭代’，我们假设初始值都为1，然后代入右手边，
等式变成：

所以猜测是错误的。然后，如果这个迭代过程会收敛，我们可以看出上面的值比初始值更接近正确值，再把上面的值代入右手边，结果为

进程继续的话，接下来的值总结在表中

对于上面这个例子，可以从表中看出当x1=1.0，x2=2.0的时候，是收敛的。
这类方程，在开始的时候，修改了[A]{x} = {b}的初始方程，将[A]的每一行和{b}向量除以[A]的相应对角线，使得修改后的系数矩阵的主对角线项为1，因此

假设矩阵[A]和右侧向量{b}保持如上所述通过对角线分割后获得的修正系数，则系数矩阵可以被划分为下面的形式

其中[I]为单位矩阵，

需要注意的是，这里的[L]和[U]和因式分解中的[L]和[U]没有关系。通过这样的定义，矩阵形式可以写成
变化之后未知数就出现在方程的两边，迭代求解的方案为：

具体例子

通过迭代解这个方程

将每个方程除以系数矩阵的主对角项，得到

因此
代入猜测值：

得到

然后

对于这个例子，收敛非常慢。为了在计算机上运行这个方案，我们需要的只是一个计算[[L]+[U]]{ x]的矩阵向量乘法、之后再加上[b]的向量加法，和一些检查收敛性的方法。
代码如下：
其中一个主程序和一个子程序，子程序为检查是否收敛的程序checkit

#使用天际线储存的乔列斯基LLT分解，使用'罚'方法的规定解
import numpy as np
import math
import B
n=3
converged=np.array([False])
xnew=np.zeros((n,1))
a=np.array([[16,4,8],[4,5,-4],[8,-4,22]],dtype=np.float)
b=np.array([[4],[2],[5]],dtype=np.float)
x=np.array([[1],[1],[1]],dtype=np.float)
tol=1.0e-5
limit=100
print('系数矩阵')
print(a[:,0])
print('右手边向量',b[:,0])
print('初始猜测值',x[:,0])
for i in range(1,n+1):
    diag=a[i-1,i-1]
    a[i-1,:]=a[i-1,:]/diag
    b[i-1,0]=b[i-1,0]/diag
a[:]=-a[:]
for i in range(1,n+1):
    a[i-1,i-1]=0
print('第一次迭代')
iters=0
while(True):
    iters=iters+1
    xnew[:]=b[:]+np.dot(a,x)
    if iters<5:
        print(x[:,0])
    B.checkit(xnew,x,tol,converged)
    if converged==True or iters==limit:
        break
    x[:,0]=xnew[:,0]
print('到收敛需要迭代次数',iters)
print('解向量',x[:,0])   

checkit

def checkit(loads,oldlds,tol,converged):
#检查前后两个量的收敛性
  neq=loads.shape[0]
  big=0.0 
  converged[:]=True
  for i in range(1,neq+1):
    if abs(loads[i-1,0])>big:
      big=abs(loads[i-1,0])
  for i in range(1,neq+1):
    if abs(loads[i-1,0]-oldlds[i-1,0])/big>tol:
      converged[:]=False

终端输出结果：

程序结果与计算值相同。