python:共轭梯度法（以希尔伯特矩阵为例）

共轭梯度法（以希尔伯特矩阵为例）

实现共轭梯度算法并且使用它解决对称矩阵——希尔伯特矩阵为系数矩阵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)