数值分析实验报告——Hilbert矩阵的求解

格式：doc
大小：261.00 KB
文档页数：7

数值分析课程实验报告

题目：病态线性方程组的求解

理论分析表明，数值求解病态线性方程组很困难。考虑求解如下的线性方程组的求解Hx = b，期中H是Hilbert矩阵，()ijnnHh，11ijhij，i，j = 1,2，…，n

1. 估计矩阵的2条件数和阶数的关系

2. 对不同的n，取(1,1,,1)nx，分别用Gauss消去，Jacobi迭代，Gauss-seidel迭代，SOR迭代和共轭梯度法求解，比较结果。

3. 结合计算结果，试讨论病态线性方程组的求解。

解答过程

1.估计矩阵的2-条件数和阶数的关系

矩阵的2-条件数定义为：1222()CondAAA，将Hilbert矩阵带入有：

1222()CondHHH

调用自编的Hilbert_Cond函数对其进行计算，取阶数n = 50，可得从1阶到50阶的2-条件数，以五位有效数字输出，其中前10项见表1。

表1.前十阶Hilbert矩阵的2-条件数

阶数 1 2 3 4 5

2-条件数 1 19.281 524.06 1.5514e+004 4.7661e+005

阶数 6 7 8 9 10

2-条件数 1.4951e+007 4.7537e+008 1.5258e+010 4.9315e+011 1.6025e+013

从表1可以看出，随着阶数每递增1，Hilbert矩阵的2-条件数都至少增加一个数量级，但难以观察出明显的相依规律。故考虑将这些数据点绘制在以n为横轴、Cond（H）2为纵轴的对数坐标系中（编程用Hilbert_Cond函数同时完成了这个功能），生成结果如图1。 2 / 7

图1.不同阶数下Hilbert矩阵的2-条件数分布

由图可见，当维数较小时，在y-对数坐标系中Cond（H）2与n有良好的线性关系；但n超过10后，线性趋势开始波动，n超过14后更是几乎一直趋于平稳。事实上，从n = 12开始，系统便已经开始提出警告：“Warning: Matrix is close to singular or badly scaled.Results

may be inaccurate.”。也就是说，当n较大时，H矩阵已经接近奇异，计算结果可能是不准确的。通过查阅相关资料，我找到了造成这种现象的原因：在matlab中，用inv函数求条件数过大的矩阵的逆矩阵将是不可靠的。而调用系统自带的专门对Hilbert矩阵求逆的invhilb（n）函数则不存在这个问题，生成结果如图2。

图2. 修正后的不同阶数下Hilbert矩阵的2-条件数分布

简便起见，取n不大于10的前十项进行线性拟合，结果如图3。

0123456789101111010010001000010000010000001E71E81E91E101E111E121E13 Cond(H) Linear Fit of Sheet1 Cond(H)Cond(H)nEquationy = a + b*xWeightNo WeightingResidual Sum

of Squares0.05264Pearson's r0.99985Adj. R-Square0.99967ValueStandard ErroCond(H)Intercept-1.651830.05541Cond(H)Slope1.478630.00893

图3.1~10阶Hilbert矩阵2-条件数的线性拟合

由拟合结果知，Cond（H）2与n的关系为：

2lg()1.651831.47863CondHn

其线性相关系数r=0.99985，可见二者具有较好的线性关系。对上式稍作变形得：

1.651831.478632()10nCondH

这就是对Hilbert矩阵的-2条件数与其阶数n的关系估计。可见Hilbert矩阵的2-条件数会随其阶数n的增加呈指数增大趋势，因此当n较大时Hilbert矩阵将是严重病态的，甚至导致matlab中inv求逆运算失真。

2.对不同的n，采用各种方法求解方程

调用自编的Calc主函数（其中包括的Hilbert函数以及b_函数可创建出对应阶数的H矩阵以及向量b，Gauss_Cal函数、Jacobi_Cal函数、Gauss_Seidel_Cal函数、SOR_Cal函数（该函数自动寻找最优松弛因子，然后以最优因子进行求解）以及CG_Cal函数则可完成各自方法的求解），分别取n = 2,5,10,20,50，对于迭代法设定终止计算精度为1010，所得计算结果以16位有效数字输出，分别见表2~表6。

表2.n=2的计算结果

求解方法迭代矩阵谱半径p；

是否收敛迭代次数N 解向量x 误差e

GAUSS消元 — — 1

0.999999999999999 8.00593208497344e-016

Jacobi迭代 0.866025403784439

收敛 169 1.000000000016

1.000000000048 5.05958147990869e-011

GS迭代 0.75

收敛 77 1.00000000015982

0.999999999760273 2.88116190608079e-010

SOR迭代 0.35

（1.35b）

收敛 26 1.00000000001783

0.999999999982226 2.51745780938445e-011

CG迭代 —

收敛 2 1

3.27844930319752e-015

从表2可看出，n=2时，四种迭代法都能够收敛，迭代次数最大为e+2量级（J法），最小仅要2次（CG法），并且五种解法都能给出非常精确的结果，最大误差为e-10量级（GS法）。

表3.n=5的计算结果

求解方法迭代矩阵谱半径p；

是否收敛迭代次数N 误差e

GAUSS消元 — — 3.47269899767798e-012

Jacobi迭代 3.44414219116596

不收敛 — —

GS迭代 0.999957671222957

收敛 231815 2.3623228567773e-006

SOR迭代 0.999190149180609

（1.83b）

收敛 21500 8.86896685436631e-008

CG迭代 —

收敛 7 5.84409588812087e-011

从表3可以看出，n=5时，J法已经不收敛，其余解法依然收敛（但值得注意的是GS法以及SOR法的谱半径以及非常接近1，达到了收敛的边缘），最大迭代次数达到e+6量级（GS法），最小需要7次（CG法），计算精度依然较高，最大误差为e-6量级。SOR法比GS法节省了较多的迭代次数。

表4.n=10的计算结果

求解方法迭代矩阵谱半径p；

是否收敛迭代次数N 误差e

GAUSS消元 — — 0.00070763269657929

Jacobi迭代 7.77981513192998

不收敛 — —

GS迭代 0.999999999997045

收敛 3086843 0.00213101923167426

SOR迭代 0.999999999871931

（1.93b）

收敛 10303801 0.00177481054386451

CG迭代 —

收敛 8 0.00116732746950012

从表4可以看出，n=10时，各种解法的收敛性与n=5时相同（GS法与SOR法的谱半径进一步趋近1），最大迭代次数达e+8量级（SOR法），最小需要8次（CG法），计算精度已经较低，最大误差达到e-3量级。此时有一个异常现象：SOR法的迭代次数不但不比GS法少，反而超出好几倍。但根据计算出的两种算法下的迭代矩阵谱半径，可以看出SOR法为0.999999999871931，而GS法为0.999999999997045，在最优松弛因子之下的SOR法确实具有更小的迭代矩阵谱半径，因此应当更快收敛。经检查，计算程序的编制应该没有错误，但计算结果确实与理论不符，希望老师指点迷津！

表5.n=20的计算结果

求解方法迭代矩阵谱半径p；

是否收敛迭代次数N 误差e

GAUSS消元 — — 105.281872512979

Jacobi迭代 16.4920989837926

不收敛 — —

GS迭代 1

收敛 12545046 0.00166739398859172

SOR迭代 1

（0.51b）

收敛 13118171 0.00111174295388826

CG迭代 —

收敛 10 0.00220319859317489

从表5可以看出，n=20时，各种解法的收敛性依然没变（但GS法以及SOR法的谱半径在matlab的十六位显示精度下已经等于1），最大迭代次数e+8量级（SOR法），最小需要10次（CG法），Gauss消元法的计算结果已失去意义（误差e+2量级），迭代法的计算误差均为e-3量级。且此时SOR的迭代次数依然高于GS。

表6.n=50的计算结果

求解方法迭代矩阵谱半径p；

是否收敛迭代次数N 误差e

GAUSS消元 — — 1697.4503939611

Jacobi迭代 42.6768950976645

不收敛 — —

GS迭代 1

收敛 18145224 0.00294769067781031

SOR迭代 1

（0.06b）

收敛 3957421 0.00150790137625883

CG迭代 —

收敛 15 0.00187865108201776

从表6可以看出，n=50时，计算结果的特点与n=20相似，不同之处在于Gauss消元法的误差进一步增大（e+3量级）。

3. 综合讨论

根据上一部分的计算结果，可以总结出求解病态矩阵的两个特点：

1.收敛性差，收敛速度慢。

对于本例，当阶数n大于2以后J法就不收敛了，这是收敛性差的一个体现。GS法和SOR法虽然一直保持收敛（事实上这是由于Hilbert矩阵是正定对称的，所以决定了GS法以及SOR法一定收敛），但随着n的增大它们的谱半径迅速趋近于1（n=20时matlab的16位显示精度已经无法显示出它们的谱半径与1的差别），根据理论我们知道，趋近于1的谱半径决定了计算的收敛速度将会非常慢，之前计算结果中的迭代次数N也很好地验证了这一点。

2.计算精度低，阶数较高时误差惊人。

根据理论，线性方程组系数矩阵的条件数的直观意义就是初始数据扰动的相对误差的放大倍数。由于数值计算时的舍入误差等因素，初始数据的微小扰动难以避免，因此解病态方程组时，这种微小的初始相对误差将被病态矩阵巨大的条件数所放大，最后造成计算结果的巨大偏差。对于本例，n=20时Gauss消元法的计算误差竟然达到e+2量级，使得计算结果完全失去意义，而其他几种迭代法的计算误差看似能够接受（e-3量级），但仔细分析上一节

数值分析综合实验报告

一、实验目的

通过本次综合实验，掌握数值分析中常用的插值方法、方程求根方法以及数值积分方法，了解这些方法在实际问题中的应用，提高数值计算能力。

二、实验内容

1. 插值方法

（1）拉格朗日插值法：利用已知数据点构造多项式，以逼近未知函数。

（2）牛顿插值法：在拉格朗日插值法的基础上，通过增加基函数，提高逼近精度。

2. 方程求根方法

（1）二分法：适用于函数在区间内有正负值的情况，通过不断缩小区间来逼近根。

（2）Newton法：利用函数的导数信息，通过迭代逼近根。

（3）不动点迭代法：将方程转化为不动点问题，通过迭代逼近根。

3. 数值积分方法

（1）矩形法：将积分区间等分，近似计算函数值的和。

（2）梯形法：将积分区间分成若干等分，用梯形面积近似计算积分。

（3）辛普森法：在梯形法的基础上，将每个小区间再等分，提高逼近精度。

三、实验步骤

1. 拉格朗日插值法

（1）输入已知数据点，构造拉格朗日插值多项式。

（2）计算插值多项式在未知点的函数值。

2. 牛顿插值法

（1）输入已知数据点，构造牛顿插值多项式。

（2）计算插值多项式在未知点的函数值。

3. 方程求根方法

（1）输入方程和初始值。（2）选择求解方法（二分法、Newton法、不动点迭代法）。

（3）迭代计算，直到满足精度要求。

4. 数值积分方法

（1）输入被积函数和积分区间。

（2）选择积分方法（矩形法、梯形法、辛普森法）。

（3）计算积分值。

四、实验结果与分析

1. 插值方法

（1）拉格朗日插值法：通过构造多项式，可以较好地逼近已知数据点。

（2）牛顿插值法：在拉格朗日插值法的基础上，增加了基函数，提高了逼近精度。

2. 方程求根方法

（1）二分法：适用于函数在区间内有正负值的情况，计算简单，但收敛速度较慢。

（2）Newton法：利用函数的导数信息，收敛速度较快，但可能存在数值不稳定问题。

（3）不动点迭代法：将方程转化为不动点问题，收敛速度较快，但可能存在初始值选择不当的问题。

实验二MATLAB矩阵分析与处理答案

1 / 3 实验二 MATLAB矩阵分析与处理

班级姓名学号

1. 设有分块矩阵33322322ERAOS，其中E、R、O、S分别为单位矩阵、随机矩阵、零矩阵和对角矩阵，试通过数值计算验证22+ERRSAOS。

>> E=eye(3);

>> R=rand(3,2);

>> O=zeros(2,3);

>> S=diag([1 2]);

>> A=[E R;O S];

>> A^2

>> B=[E R+R*S;O S^2]

2. 产生5阶希尔伯特矩阵H和5阶帕斯卡矩阵P，且求其行列式的值Hh和Hp以及它们的条件数Th和Tp，判断哪个矩阵性能更好。为什么？

>> H=hilb(5);

>> P=pascal(5);

>> Hp=det(P);

>> Hh=det(H);

>> Th=cond(H,2)

>> Tp=cond(P,2)

3. 建立一个55矩阵，求它的行列式值、迹、秩和范数。

>> A=rand(5);

>> det(A);

>> trace(A);

>> rank(A);

>> norm(A,1);

>> norm(A,2);

>> norm(A,inf);

4. 已知

2961820512885A

求A的特征值及特征向量。

>> A=[-29 6 18;20 5 12;-8 8 5];

>> [V,D]=eig(A)

V =

0.7130 0.2803 0.2733

-0.6084 -0.7867 0.8725

0.3487 0.5501 0.4050 2 / 3 D =

-25.3169 0 0

0 -10.5182 0

基于Matlab的数值分析实验的设计r——以矩阵的奇异值分解为例

杨海涛

【摘要】基于Matlab的数学实验,目的是为了提高普通高校学生的数学兴趣和动手能力,为培养解决实际数学问题的能力而设计的实验方法和步骤.利用矩阵的奇异值分解svd函数的数值算法设计了一种在实验课上用svd算法仿真的实验方案.

【期刊名称】《内蒙古民族大学学报（自然科学版）》

【年(卷),期】2016(031)006

【总页数】4页(P465-468)

【关键词】数学实验;数值计算;Matlab;奇异值分解

【作者】杨海涛

【作者单位】内蒙古民族大学数学学院,内蒙古通辽 028043

【正文语种】中文

【中图分类】O241.4

利用计算机计算解决实际问题时需要建立数学模型，同时也需要相应数学软件的支持，为了更好的掌握数学软件这一有利的工具，本文通过Matlab中的矩阵的奇异值分解〔1〕svd函数设计了如何学习和掌握该函数的实验方法，从而可以进一步学习其它函数.目前，Matlab在计算仿真和图形图像的压缩处理〔2～4〕有着广泛的应用，希望这个实验方法能够给正在学习数学软件的学生起到抛砖引玉的作用. 1.1 理论基础

任意复长方矩阵的奇异值分解定理如下（Autonne-Eckart-Young定理）：〔5〕

设A为一个m×n的实矩阵（m≥n），则存在正交矩阵V和U，使得：

其中∑=diag(σ1,σ2,…,σr)，而且σ1≥σ2≥…≥σr≥0.

奇异值分解在统计分析、信号与图像处理、系统理论和控制中有着广泛的应用.

1.2 SVD的数值算法〔6〕

求解一般矩阵奇异值问题最常用的算法可分为两大类：一类是QR分解的方法；另一类是Jacobi旋转的方法.

奇异值分解的QR分解法分为两个阶段：

一是矩阵的二重对角化，即利用Householder变换把矩阵A变换成二重对角矩阵.这个过程就是追赶（chasing）过程〔7〕.

二是QR分解，即保持二重对角矩阵形式不变，利用正交变换使对角线元素逐渐减小，使得矩阵接近对角矩阵.

矩阵式实验报告

一、实验目的

1. 理解矩阵的基本概念和性质。

2. 掌握矩阵的运算方法，包括加法、减法、乘法等。

3. 学习矩阵的应用，如线性方程组的求解。

4. 提高数学建模和解决问题的能力。

二、实验内容

本次实验主要围绕矩阵的运算和应用展开，具体内容包括：

1. 矩阵的加法与减法

2. 矩阵的乘法

3. 矩阵的逆

4. 线性方程组的求解

三、实验步骤

1. 矩阵的加法与减法

（1）选择两个矩阵A和B，确保它们具有相同的行数和列数。

（2）将矩阵A和B对应位置的元素相加或相减，得到新的矩阵C。

（3）验证矩阵C的行数和列数与矩阵A和B相同。

2. 矩阵的乘法

（1）选择两个矩阵A和B，确保矩阵A的列数等于矩阵B的行数。

（2）计算矩阵A的每一行与矩阵B的每一列的点积，得到新的矩阵C。

（3）验证矩阵C的行数等于矩阵A的行数，列数等于矩阵B的列数。

3. 矩阵的逆

（1）选择一个可逆矩阵A。

（2）使用高斯-约当消元法求解矩阵A的逆。（3）验证矩阵A与其逆矩阵的乘积为单位矩阵。

4. 线性方程组的求解

（1）选择一个线性方程组，例如：AX = B，其中A是系数矩阵，X是未知数矩阵，B是常数矩阵。

（2）使用高斯-约当消元法求解线性方程组。

（3）验证求解得到的X矩阵是否满足原方程组。

四、实验结果与分析

1. 矩阵的加法与减法

通过实验，我们发现矩阵的加法与减法运算满足交换律和结合律，且结果矩阵的行数和列数与原矩阵相同。

2. 矩阵的乘法

实验结果表明，矩阵的乘法运算满足交换律和结合律，且结果矩阵的行数等于第一个矩阵的行数，列数等于第二个矩阵的列数。

3. 矩阵的逆

实验发现，对于可逆矩阵，其逆矩阵存在，且满足A A^(-1) = A^(-1) A = E（单位矩阵）。

4. 线性方程组的求解

通过高斯-约当消元法，我们成功求解了线性方程组，并验证了求解结果的正确性。

五、实验结论

1. 理解了矩阵的基本概念和性质，掌握了矩阵的运算方法。

下载提示

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

数值分析实验报告——Hilbert矩阵的求解

合集下载

数值分析综合实验报告

实验二MATLAB矩阵分析与处理答案

基于Matlab的数值分析实验的设计r——以矩阵的奇异值分解为例

矩阵式实验报告

文档推荐

最新文档