10-16：数学课小结

1.Hilbert矩阵

Hilbert矩阵条件数随阶数变化.png

2.SOR迭代法抢救Hilbert矩阵

def MySOR(A,b,initial,delta,w,isw=False,isPrint=False):
    '''
    @author:zengwei
    在高斯赛德尔迭代法基础上多加了一行
    '''
    D = np.diag(np.diag(A)) 
    L = -np.tril(A,-1)
    U = -np.triu(A,1)
    
    if isw==True:
        BJ = np.dot(np.linalg.inv(D),L+U)
        lamdaBJ,_ = np.linalg.eig(BJ)
        rouBJ = np.max(np.abs(lamdaBJ))
        w = 2./(1+np.sqrt(1-rouBJ**2))
        print('最优松弛因子为%f'%w)
    else:
        print('人为设置松弛因子')
        
    d = np.linalg.inv( D - L )           

    BG = np.dot(d,U)                     # 迭代矩阵BG
    lamda,_ = np.linalg.eig(BG)
    if np.max(np.abs(lamda))<1:          # 谱半径小于1
        f = np.dot(d,b)
        X = np.dot( BG ,initial ) + f

        k = -np.log(delta)/-np.log(np.max(np.abs(lamda)))

        BGnorm = np.linalg.norm(BG)
        times = 1                        # 因为前面有了一次迭代
        while np.linalg.norm(X - initial,ord=np.inf) > delta:
            initial = X
            X = np.dot( BG,initial )+f  
            X = w*X + (1-w)*initial      # 为实现SOR多加的一行
            times = times +1
        if isPrint==True:
            print('谱半径为：',np.max(np.abs(lamda)))
            print('理论上的最大迭代次数为：%d次' %(int(k)+1))
            print("实际上的最大迭代次数为：%d次" %times)
        return X
    else:
        print('Sorry,不可收敛。')
        print('谱半径为：',np.max(np.abs(lamda)))

求解一个8阶的Hilbert矩阵，解为 $[0.2 0.2 ... 0.2]$ ，初始解设为0，看看不同的松弛因子 $w$ 会不会使得迭代加快，和能否迭代到我们想要的精度。

N = 8
A = Hilbert(N)
x = np.ones(N)*0.2
initial = np.zoros(N)
b = np.dot(A,x)
delta = 10**(-6)

w = 1
MySOR(A,b,initial,delta,w,isPrint=True)

w = 1.5
MySOR(A,b,initial,delta,w,isPrint=True)

w = 1.99
MySOR(A,b,initial,delta,w,isPrint=True)

结论：可以得到最优的收敛速度，但是近似值无法得到 $10^{-6}$ 的精度。
这是想说明什么呢？我认为这是在说明Hilbert矩阵是个病态矩阵，想让它收敛到想要的精度太难了。雅克比迭代直接歇菜，高斯赛德尔迭代需要巨大的迭代次数，超松弛SOR法虽然可以通过调节松弛因子来加快收敛速度，但依然得不到想要精度的解。由此可见，这三种迭代方法分别为狗、虎、龙。

3.奇怪的稀疏矩阵

稀疏矩阵.png

这个矩阵是严格对角最优矩阵，也是对称正定三对角矩阵。那我们使用SOR求解时可以找到最优的松弛因子

w

。对于对称正定三对角矩阵，还有一个有意思的结论：

p(B_G)=[p(B_j)]^2

，高斯赛德尔迭代比雅克比迭代法收敛速度快

1/2

。
生成矩阵：

Aii = np.array([[4,-1,0],[-1,4,-1],[0,-1,4]])
n = Aii.shape[0]
I = -np.identity(n)
O = np.zeros((n,n))
col1 = np.array([Aii,I,O]).reshape(n**2,n)
col2 = np.array([I,Aii,I]).reshape(n**2,n)
col3 = np.array([O,I,Aii]).reshape(n**2,n)
A = np.concatenate((np.concatenate((col1,col2),axis=1),col3),axis=1)

b = np.zeros(n**2)
b[int(np.floor(n**2/2))]=1

用SOR求解：

initial = np.zeros(n**2)
MySOR(A,b,initial,delta,w,isw=True,isPrint=True)

w = 0.5
MySOR(A,b,initial,delta,w,isw=False,isPrint=True)

计算得到最优的松弛因子为1.171573，是需迭代15次。分别用松弛因子为0.5,1.0,1.5,1.9来测试，迭代次数分别为38次、18次、19次和113次。

最后编辑于：2020.10.17 10:50:35