数值分析第二次作业-求解系数矩阵为Hilbert 矩阵的线性方程组

1、问题

现要求解系数矩阵由16 阶Hilbert 方程组构成的线性方程组，右端项为

 即要求解方程组Ax = b，其中 A=A0,b=b0

 分别用高斯-赛德尔方法、最速下降法、共轭梯度法求解如下。

2.1. 高斯-赛德尔方法 

 

 2.2. 最速下降法

 

 2.3. 共轭梯度法

 在最速下降法中，搜索方向p取的是函数减少最快的方向，负梯度方
向从局部上看是二次函数的最快下降方向，但是整体上看这并非最好。对于对称
正定矩阵A，共轭梯度法选择关于A 共轭的向量p1, p2, …代替最速下降法中的
负梯度方向，使迭代法对任意给定的初始点x0具有有限步收敛性。共轭梯度法
本质上是利用 A正交定理，反复对残差实施施密特正交化。

 3、对比

 绘制出最速下降法的搜索示意图如图3.1，梯度下降法中，相邻残差向量是
相互正交的，从图中可以看出，每一步的向量都与前后步的向量正交。

 在此次实验中，由于Hilbert 矩阵的病态性，高斯-赛德尔迭代法、最速下降
法、共轭梯度法表现各不相同。高斯-赛德尔迭代法和最速下降法迭代次数较高
且结果不稳定，共轭梯度法以超线性的效率计算出结果，仅迭代3 次。

 4、代码实现

# *coding:utf-8 *
import numpy as np
# A_0是16阶hilbert矩阵
A_0 = 1. / (np.arange(1, 17) + np.arange(0, 16)[:, np.newaxis])
A = np.array(A_0)
b = np.array([2877.0 / 851, 3491.0 / 1431, 816.0 / 409, 2035.0 / 1187,
              2155.0 / 1423, 538.0 / 395, 1587.0 / 1279, 573.0 / 502,
              947.0 / 895, 1669.0 / 1691, 1589.0 / 1717, 414.0 / 475,
              337.0 / 409, 905.0 / 1158, 1272.0 / 1711, 173.0 / 244])
x0 = np.array(np.zeros(16))
x = np.array(np.zeros(16))
r = b - np.dot(A, x)
p = r

def GS(x0):
    k = 0
    while k < 10000:
        for i in range(16):
            temp = 0
            temp_x = x0.copy()
            for j in range(16):
                if i != j:
                    temp += x[j] * A[i][j]
            x[i] = (b[i] - temp) / A[i][i]
            x0[i] = x[i].copy()
        norm = np.linalg.norm(x - temp_x)
        k += 1
        if norm < 1e-4:
            break
        else:
            x0 = x.copy()
            print("高斯-赛德尔迭代法：迭代次数k =", k)
            print(x)

def GD(x):
    k = 0
    while k < 10000:
        x0 = x.copy()
        k += 1
        r = b - np.dot(A, x0)
        a = np.dot(r.T, r) / np.dot(r.T, np.dot(A, r))
        x = x0 + a * r
        norm = np.linalg.norm(x0 - x)
        if norm < 1e-4:
            break
        print("最速下降法：迭代次数k =", k)
        print(x)

def CG(x, r, p):
    k = 0
    while True:
        k += 1
        r1 = r
        a = np.dot(r.T, r) / np.dot(p.T, np.dot(A, p))
        x = x + a * p
        r = r - a * np.dot(A, p)
        er = np.linalg.norm(np.dot(A, x) - b) / np.linalg.norm(b)
        if er < 1e-4:
            break
        else:
            beta = np.linalg.norm(r) ** 2 / np.linalg.norm(r1) ** 2
            p = r + beta * p
        print("共轭梯度法：迭代次数k=", k)
        print(x)

if __name__ == '__main__':
    # GS(x0)
    # GD(x)
    CG(x, r, p)