数值分析第二次作业-求解系数矩阵为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)

  • 12
    点赞
  • 51
    收藏
    觉得还不错? 一键收藏
  • 1
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值