python实现希尔伯特矩阵的一般迭代及SOR迭代

数值分析 李庆阳

x*=(1,1,1,1,……，1)T，b=Hn·x*

n=6,8,10 w=1,1.25,1.5,求解

使用矩阵方法求解,不适用于实际案例，仅用于方法复现

import numpy as np            #222100200668

P=0.01                        #设置精度要求

def Hn(n):                    #希尔伯特矩阵生成函数
    H=np.empty([n,n])
    for i in range(0,n):
        for j in range(0,n):
            H[i][j]=(1/(i+j+1))
    return H

def SR(M):                   #谱半径求解函数
    a,b=np.linalg.eig(M)     #a为特征值集合，b为特征值向量
    return np.max(np.abs(a)) #返回最大特征值

def ND(n,P):                #一般迭代
    H=Hn(n)
    B=np.empty([n,n])
    X=np.ones([n,1])         #X为x的精确值
    b=np.dot(H,X)
    x0=np.zeros([n,1])       #x0为初始值，赋0
    f=np.empty([n,1])
    for i in range(n):
        for j in range(n):
            if i==j:
                B[i][j]=0
            else:
                B[i][j]=(-1)*(H[i][j]/H[i][i])
    M=SR(B)                  #谱半径不小于1发散
    if M>=1:
            print(f'n={n}时，Lw谱半径={M}一般迭代法发散')
    else:
        for i in range(n):
            f[i][0]=b[i][0]/H[i][i]
        x=np.dot(B,x0)+f
        for k in range(200000):
            if np.linalg.norm(X-x,ord=np.inf)<P:       #无穷范数小于精度值终止迭代
                break
            if np.linalg.norm(X-x,ord=np.inf)>=P:
                x=np.dot(B,x)+f
            print(f'n={n},w={w}时，SOR迭代{k}次得到\nx={x}\n误差为{np.linalg.norm(X-x,ord=np.inf)}\n')
    

def SOR(n,P):
    H=Hn(n)
    X=np.ones([n,1])
    b=np.dot(H,X)
    D=np.zeros([n,n])        #D对角元素矩阵
    L=np.zeros([n,n])        #L下三角矩阵
    U=np.zeros([n,n])        #U上三角矩阵
    x0=np.zeros([n,1])
    for i in range(n):
        for j in range(n):
            if i==j:
                D[i][i]=H[i][i]
            if i>j:
                L[i][j]=-1*H[i][j]
            if i<j:
                U[i][j]=-1*H[i][j]
    
    for w in np.arange(1,1.6,0.25):
        Lw=np.dot(np.linalg.inv(D-w*L),(1-w)*D+w*U)
        M=SR(Lw)
        if M>=1:
            print(f'n={n}时，Lw谱半径={M}SOR迭代法发散')
        else:
            f=np.dot(w*np.linalg.inv(D-w*L),b)
            x=np.dot(Lw,x0)+f
        for k in range(200000):
            if np.linalg.norm(X-x,ord=np.inf)<P:
                break
            if np.linalg.norm(X-x,ord=np.inf)>=P:
                x=np.dot(Lw,x)+f
        print(f'n={n},w={w}时，SOR迭代{k}次得到\nx={x}\n误差为{np.linalg.norm(X-x,ord=np.inf)}\n')


for n in range(6,11,2):
    ND(n,P)

  
for n in range(6,11,2):
    SOR(n,P)

结果为：