线性方程组的迭代法 python代码实现

第一次写博客，时间不多，所以数学原理什么的就不写了。

直接贴代码和实例演示，以下代码基于python和numpy。

解释就一张图。

设有线性方程组 Ax = b，其中A非奇异, 且 aii不等于0, (i = 1, 2, …, n )。

这里十分详细的讲解了“迭代法求解线性方程组”。

往期博客：

函数插值法之牛顿插值法 python代码实现

数值积分 python代码实现

数值微分 python代码实现

python计算信息熵、信息增益和信息增益率

Jacobi迭代法

定义函数

首先

import numpy as np

以下便是我定义的函数：

def jacobi(a,b,c=0.0001,d=30):
    x1=np.zeros(a.shape[1])
    x2=np.zeros(a.shape[1])
    k=0
    while k<d:
        k=k+1
        print('k=',k)
        for i in range(a.shape[1]):
            x2[i]=(-a[i].dot(x1)+b[i]+a[i,i]*x1[i])/a[i,i]
        if np.max(np.abs(x2-x1))<=c:
            print("x%d=" % k,x2)
            print(np.max(np.abs(x2-x1)))
            break
        print("x%d=" % k,x2)
        x1=x2.copy()
	return x2

参数说明

“a”就是“A”,“b”就是“b”，和上图对应。

c是一个跳出循环的限定值。

d是最大迭代次数。

实例运行

拿个实例运行一下。

a=np.array([[2,-1,0,0],
          [-1,2.5,-1,0],
           [0,-1,2.5,-1],
           [0,0,-1,2]])
b=np.array([-4,4,4,-4])
jacobi(a,b)

得出以下结果：

k= 1
x1= [-2. 1.6 1.6 -2. ]
k= 2
x2= [-1.2 1.44 1.44 -1.2 ]
k= 3
x3= [-1.28 1.696 1.696 -1.28 ]
k= 4
x4= [-1.152 1.7664 1.7664 -1.152 ]
k= 5
x5= [-1.1168 1.84576 1.84576 -1.1168 ]
k= 6
x6= [-1.07712 1.891584 1.891584 -1.07712 ]
k= 7
x7= [-1.054208 1.9257856 1.9257856 -1.054208 ]
k= 8
x8= [-1.0371072 1.94863104 1.94863104 -1.0371072 ]
k= 9
x9= [-1.02568448 1.96460954 1.96460954 -1.02568448]
k= 10
x10= [-1.01769523 1.97557002 1.97557002 -1.01769523]
k= 11
x11= [-1.01221499 1.98314992 1.98314992 -1.01221499]
k= 12
x12= [-1.00842504 1.98837397 1.98837397 -1.00842504]
k= 13
x13= [-1.00581301 1.99197957 1.99197957 -1.00581301]
k= 14
x14= [-1.00401021 1.99446662 1.99446662 -1.00401021]
k= 15
x15= [-1.00276669 1.99618256 1.99618256 -1.00276669]
k= 16
x16= [-1.00190872 1.99736635 1.99736635 -1.00190872]
k= 17
x17= [-1.00131683 1.99818305 1.99818305 -1.00131683]
k= 18
x18= [-1.00090847 1.99874649 1.99874649 -1.00090847]
k= 19
x19= [-1.00062675 1.99913521 1.99913521 -1.00062675]
k= 20
x20= [-1.0004324 1.99940338 1.99940338 -1.0004324 ]
k= 21
x21= [-1.00029831 1.99958839 1.99958839 -1.00029831]
k= 22
x22= [-1.0002058 1.99971603 1.99971603 -1.0002058 ]
k= 23
x23= [-1.00014198 1.99980409 1.99980409 -1.00014198]
8.805852909499201e-05

Gauss-Seidel迭代法

定义函数

定义的函数如下：

def Gauss_Seidel(a,b,c=0.0001,d=30):
    x1=np.zeros(a.shape[1])
    x2=np.zeros(a.shape[1])
    k=0
    while k<d:
        k=k+1
        print('k=',k)
        for i in range(a.shape[1]):
            x2[i]=(-a[i].dot(x2)+b[i]+a[i,i]*x2[i])/a[i,i]
        if np.max(np.abs(x2-x1))<=c:
            print("x%d=" % k,x2)
            print(np.max(np.abs(x2-x1)))
            break
        print("x%d=" % k,x2)
        x1=x2.copy()
	return x2

实例运行

拿同一个实例运行一下。

a=np.array([[2,-1,0,0],
          [-1,2.5,-1,0],
           [0,-1,2.5,-1],
           [0,0,-1,2]])
b=np.array([-4,4,4,-4])
Gauss_Seidel(a,b)

得到如下结果：

k= 1
x1= [-2. 0.8 1.92 -1.04]
k= 2
x2= [-1.6 1.728 1.8752 -1.0624]
k= 3
x3= [-1.136 1.89568 1.933312 -1.033344]
k= 4
x4= [-1.05216 1.9524608 1.96764672 -1.01617664]
k= 5
x5= [-1.0237696 1.97755085 1.98454968 -1.00772516]
k= 6
x6= [-1.01122458 1.98933004 1.99264195 -1.00367902]
k= 7
x7= [-1.00533498 1.99492279 1.99649751 -1.00175125]
k= 8
x8= [-1.0025386 1.99758356 1.99833293 -1.00083354]
k= 9
x9= [-1.00120822 1.99884988 1.99920654 -1.00039673]
k= 10
x10= [-1.00057506 1.99945259 1.99962234 -1.00018883]
k= 11
x11= [-1.0002737 1.99973946 1.99982025 -1.00008987]
k= 12
x12= [-1.00013027 1.99987599 1.99991445 -1.00004278]
k= 13
x13= [-1.000062 1.99994098 1.99995928 -1.00002036]
6.826783096514077e-05

SOR迭代法

定义函数

定义的函数如下：

def SOR(a,b,w=1,c=0.0001,d=30):
    x1=np.zeros(a.shape[1])
    x2=np.zeros(a.shape[1])
    k=0
    while k<d:
        k=k+1
        print('k=',k)
        for i in range(a.shape[1]):
            x2[i]=(-a[i].dot(x2)+b[i])*w/a[i,i]+x2[i]
        if np.max(np.abs(x2-x1))<=c:
            print("x%d=" % k,x2)
            print(np.max(np.abs(x2-x1)))
            break
        print("x%d=" % k,x2)
        x1=x2.copy()
	return x2

参数说明

w是松弛因子。

实例运行

拿同一个实例运行一下。

a=np.array([[2,-1,0,0],
          [-1,2.5,-1,0],
           [0,-1,2.5,-1],
           [0,0,-1,2]])
b=np.array([-4,4,4,-4])
w=1.16
SOR(a,b,w)

得到如下结果：

k= 1
x1= [-2.32 0.77952 2.21769728 -1.03373558]
k= 2
x2= [-1.4966784 2.06582956 1.98006004 -1.00616748]
k= 3
x3= [-0.88235031 2.03480459 2.01647801 -0.98945596]
k= 4
x4= [-0.99863729 2.00270936 2.0035131 -0.99964945]
k= 5
x5= [-0.9986466 2.00182455 2.00044715 -0.99979674]
k= 6
x6= [-0.9991583 2.0003061 2.0001648 -0.99993694]
k= 7
x7= [-0.99995713 2.00004738 2.00002488 -0.99999566]
k= 8
x8= [-0.99997938 2.00001353 2.00000431 -0.99999819]
3.3849332815805155e-05

总结

综上可知，同一个实例，在速度上，SOR迭代法>Gauss-Seidel迭代法>Jacobi迭代法。当然，前提是三者均需满足收敛条件。

迭代法收敛的充要条件如下图所示。

以上便是我（为了不用一步一步按计算器来写作业）写的三种迭代方法的定义函数，可能会有更简单，可以评论一番，多多加交流。

本来是为了补数学作业，想GKD，没想到出了点bug（根本就没几行代码，还出了bug。。。），吃了个中午饭才大概想到了“x1=x2”应该有问题。

不过塞翁司马，焉知祸福，填补了空白也不错。

可谓是前人栽树，后人乘凉。

扩展

2020.3.27

不少同学注意到了我在函数中，默认的初始值为0向量，问我能不能自己定义初始值，怎么定义初始值。
所以我对Jacobi函数做了如下修改。

def jacobi_pro(a,b,x0=np.zeros(a.shape[1]),c=0.0001,d=30):
    print('x0=',x0)
    x1=np.zeros(a.shape[1])
    k=0
    while k<d:
        k+=1
        print('k=',k)
        for i in range(a.shape[1]):
            x1[i]=(-a[i].dot(x0)+b[i]+a[i,i]*x0[i])/a[i,i]
        if np.max(np.abs(x0-x1))<=c:
            print("x%d=" % k,x1)
            print(np.max(np.abs(x0-x1)))
            break
        print("x%d=" % k,x1)
        x0=x1.copy()
    return x1

拿同一个实例运行一下，将初始值设为单位向量。

a=np.array([[2,-1,0,0],
          [-1,2.5,-1,0],
           [0,-1,2.5,-1],
           [0,0,-1,2]])
b=np.array([-4,4,4,-4])
x=np.ones(a.shape[1])
jacobi_pro(a,b,x)

得到如下结果：

x0= [1. 1. 1. 1.]
k= 1
x1= [-1.5 2.4 2.4 -1.5]
k= 2
x2= [-0.8 1.96 1.96 -0.8 ]
k= 3
x3= [-1.02 2.064 2.064 -1.02 ]
k= 4
x4= [-0.968 2.0176 2.0176 -0.968 ]
k= 5
x5= [-0.9912 2.01984 2.01984 -0.9912 ]
k= 6
x6= [-0.99008 2.011456 2.011456 -0.99008 ]
k= 7
x7= [-0.994272 2.0085504 2.0085504 -0.994272 ]
k= 8
x8= [-0.9957248 2.00571136 2.00571136 -0.9957248 ]
k= 9
x9= [-0.99714432 2.00399462 2.00399462 -0.99714432]
k= 10
x10= [-0.99800269 2.00274012 2.00274012 -0.99800269]
k= 11
x11= [-0.99862994 2.00189497 2.00189497 -0.99862994]
k= 12
x12= [-0.99905251 2.00130601 2.00130601 -0.99905251]
k= 13
x13= [-0.99934699 2.0009014 2.0009014 -0.99934699]
k= 14
x14= [-0.9995493 2.00062176 2.00062176 -0.9995493 ]
k= 15
x15= [-0.99968912 2.00042899 2.00042899 -0.99968912]
k= 16
x16= [-0.99978551 2.00029595 2.00029595 -0.99978551]
k= 17
x17= [-0.99985203 2.00020418 2.00020418 -0.99985203]
9.17709429146818e-05

这个结果与最初的结果相比，迭代次数少了26.1%，所以初始值对迭代次数的影响也是蛮大的。

其它两种迭代法的修改方式和Jacobi是一样，在此我就不再赘述了。