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)))
``````

``````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)
``````

3.奇怪的稀疏矩阵

``````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
``````

``````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)
``````