共轭梯度法(以希尔伯特矩阵为例)
实现共轭梯度算法并且使用它解决对称矩阵——希尔伯特矩阵为系数矩阵A时,右边常数项列向量为单位列向量,初始点为零向量。展示维度分别为5,8,12,20时的迭代次数与结果。
迭代终止条件为残差的二范数小于1e-6
import numpy as np
def cg(x0, A, b):
r0 = np.dot(A, x0) - b
p0 = -r0
rk = r0
pk = p0
xk = x0
t = 0 #记录迭代次数
while np.linalg.norm(rk) >= 1e-6:
rr = np.dot(rk.T, rk)
ak = rr / np.dot(np.dot(pk.T, A), pk)
xk = xk + ak * pk
rk = rk + ak * np.dot(A, pk)
bk = np.dot(rk.T, rk) / rr
pk = -rk + bk * pk
t += 1
return xk, t
def Hilbert(n):
sq = np.zeros((n,n))
for i in range(n):
for j in range(n):
sq[i][j] = 1 / (i + j + 1)
return sq
if __name__=='__main__':
sque = [5, 8, 12, 20]
for n in sque:
b = np.ones((n,1))
A = Hilbert(n)
x0 = np.zeros((n,1))
x,t = cg(x0, A, b)
print(n,"阶时,近似值为:",x," 迭代次数为:",t)