逐次超松弛迭代法

对方程组

 进行SOR迭代的步骤如下

1. 分解A

 

=

2. 取松弛因子

 

 

3. 取初始向量

4. 迭代

5. 例题

用SOR迭代法求解方程组

 

Python代码：

import numpy as np
import math
import sys
#分解矩阵
def DLU(A):
    D=np.zeros(np.shape(A))
    L=np.zeros(np.shape(A))
    U=np.zeros(np.shape(A))
    for i in range(A.shape[0]):
        D[i,i]=A[i,i]
        for j in range(i):
            L[i,j]=-A[i,j]
        for k in list(range(i+1,A.shape[1])):
            U[i,k]=-A[i,k]
    L=np.mat(L)
    D=np.mat(D)
    U=np.mat(U)
    return D,L,U

#迭代
def SOR(A,b,x0,w,maxN,p):  #x0为初始值，maxN为最大迭代次数，p为允许误差
    D,L,U=DLU(A)
    if len(A)==len(b):
        M=np.linalg.inv(D-w*L)
        M=np.mat(M)
        B=M * ((1-w)*D + w*U)
        B=np.mat(B)
        f=w*M*b
        f=np.mat(f)
    else:
        print('维数不一致')
        sys.exit(0)  # 强制退出
    
    a,b=np.linalg.eig(B) #a为特征值集合，b为特征值向量
    c=np.max(np.abs(a)) #返回谱半径
    if c<1:
        print('迭代收敛')
    else:
        print('迭代不收敛')
        sys.exit(0)  # 强制退出
#以上都是判断迭代能否进行的过程，下面正式迭代
    k=0
    while k<maxN:
        x=B*x0+f
        k=k+1
        eps=np.linalg.norm(x-x0,ord=2)
        if eps<p:
            break
        else:
            x0=x
    return k,x

A = np.mat([[-4,1,1,1],[1,-4,1,1],[1,1,-4,1],[1,1,1,-4]])
b = np.mat([[1],[1],[1],[1]])
x0=np.mat([[0],[0],[0],[0]])
maxN=100
p=0.00001
w=1.3
print("原系数矩阵a:")
print(A, "\n")
print("原值矩阵b:")
print(b, "\n")
k,x=SOR(A,b,x0,w,maxN,p)
print("迭代次数")
print(k, "\n")
print("数值解")
print(x, "\n")

输出

系数矩阵a:
[[-4  1  1  1]
 [ 1 -4  1  1]
 [ 1  1 -4  1]
 [ 1  1  1 -4]]

原值矩阵b:
[[1]
 [1]
 [1]
 [1]]

迭代收敛
迭代次数
12

数值解
[[-1.00000152]
 [-0.99999922]
 [-1.00000012]
 [-1.00000052]]