解线性方程组的迭代法-Gauss-Seidel迭代和SOR迭代的python实现

直接上代码
1） 通用程序
本程序主要包含三个程序
GaussS(a,b, x): Gauss-Seidel的子程序,输入量包括系数矩阵a，等式右侧矩阵b，未知数x的初始值矩阵0
Relax(a,b,x, omiga)SOR法迭代的通用程序，omiga是松弛因子
Reform(a,b)列主元法将方程进行变换，以避免对角线元素为0。
2） 程序里面包含一个例子。 Gauss-Seidel方法求根，误差取方程ax-b的二范数err<10-6迭代步数为22，SOR方法求根，误差同样取方程ax-b的二范数err<10-6，松弛因子取1.2，最后迭代步数为13次，可见SOR方法更快

import math
import numpy as np  
from numpy import * 

def Reform(a,b):  #列主元
    p=np.column_stack((a,b))
    row=p.shape[0]
    for m in range(0,row):
        if m<row:
            big=np.argmax(abs(p[m:,m]))     #找到最大元素对应的行。最后一个元素不用找最大值；只找对角线下方的主元
        else:
            big=0
        b1=big+m                      #   下面几行是主元和对角线所在行交换
        c= np.copy(p[b1,:])  #  更方便的方法：   a[[b1,j], :] = a[[j,b1], :]
        p[b1,:]=p[m,:]
        p[m,:]=c
    return p
    
def GaussS(a,b, x):
    p=Reform(a,b)     #预处理，将每一列最大值移到对角线，这里的p最后一行是b
    row=p.shape[0]    #获取行数
    a0=p[:,0:row]     #系数矩阵
    b0=p[:,row]       #b矩阵
 #   print('p',p)
    j=0
    err=100.        #初始化
    print("Gauss-S method")
    while(err>1.e-6 and j <2500): #控制迭代次数
        i=0
        while( i < x.size):  #控制循环次数
            if a0[i,i]==0:
                print('a[i,i]=0, i=',i)
            x[i]=-(np.dot(a0[i,:],x)-b0[i]-a0[i,i]*x[i])/a0[i,i]
            i=i+1
        j=j+1
        err=Norm(a0,b0, x)  #计算ax-b二范数
 #       print(j,x)
    return x,j

def Relax(a,b, x,omiga):
    p=Reform(a,b)     #预处理，将每一列最大值移到对角线，这里的p最后一行是b
    row=p.shape[0]    #获取行数
    a0=p[:,0:row]     #系数矩阵
    b0=p[:,row]       #
    err=100.
    j=0
    while(err>1.e-6 and j < 2500): 
        i=0
        while( i < x.size):  #控制循环次数
            if a0[i,i]==0:
                print('a[i,i]=0, i=', i)
            x[i]=(1-omiga)*x[i]-omiga*(np.dot(a0[i,:],x)-b0[i]-a0[i,i]*x[i])/a0[i,i]
            i=i+1
        j=j+1
        err=Norm(a0,b0, x) #计算ax-b二范数
 #       print(j, x)
    return x, j

def Norm(a,b, x):  #   计算范数
    axb=np.dot(a,x)-b
    normaxb=np.linalg.norm(axb, ord=2)   #二范数
    return normaxb

#主程序

n=9
a=np.array([
[31.,-13.,0. ,0. , 0.,-10., 0. , 0. , 0. ],
[-13.,35.,9. ,0. ,-11,  0., 0. , 0. , 0. ],
[ 0., -9.,31 ,-10.,0. ,  0., 0. , 0. , 0. ],
[ 0., 0. ,-10.,79.,-30.,  0., 0. , 0. ,-9.],
[ 0., 0. ,0. ,-30.,57., -7., 0. , -5., 0. ],
[ 0., 0. ,0. ,0. ,-7.,  47., -30., 0. , 0.],
[ 0., 0. ,0. ,0. ,0. , -30., 41., 0. , 0. ],
[ 0., 0. ,0. ,0. ,-5. ,  0., 0. , 27, -2. ],
[ 0., 0. ,0. ,-9.,0. ,  0., 0. , -2. , 29.],
]) 
b=np.array([-15.,27.,-23.,0.,-20.,12.,-7.,7.,10])
x=np.zeros(n,dtype=float)
print('a',a)
print('b',b)

loosefactor=1.2
result,j=Relax(np.copy(a),np.copy(b),np.copy(x), 1.2)

print('SOR method j=',j)
print('x=',result)
result2,j=GaussS(np.copy(a),np.copy(b),np.copy(x))

print('GS method j=',j)
print('x=',result2)
print(x)