matlab计算方法实验指导误差分析

实验一 误差分析
实验1（病态问题）
实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。

对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。

通过本实验可获得一个初步体会。

数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。

病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。

问题提出：考虑一个高次的代数多项式
)1.1()
()20()2)(1()(20
1∏=-=---=k k x x x x x p
显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。

现考虑该多项式的一个扰动
)2.1(0
)(19=+x x p ε
其中ε是一个非常小的数。

这相当于是对（1.1）中19x 的系数作一个小的扰动。

我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。

实验内容：为了实现方便，我们先介绍两个MATLAB 函数：“roots ”和“poly ”。

roots(a)u =
其中若变量a 存储n+1维的向量，则该函数的输出u 为一个n 维的向量。

设a 的元素依次为121,,,+n a a a ，则输出u 的各分量是多项式方程
01121=+++++-n n n n a x a x a x a
的全部根；而函数 poly(v)b =
的输出b 是一个n+1维向量，它是以n 维向量v 的各分量为根的多项式的系数。

可见“roots ”和“poly ”是两个互逆的运算函数。

)
)20:1((;
)2();
21,1(;000000001.0ve poly roots ess ve zeros ve ess +===
上述简单的MATLAB 程序便得到（1.2）的全部根，程序中的“ess ”即是（1.2）中的ε。

实验要求：
（1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。

如果扰动项的系数ε很小，我们自然感觉（1.1）和（1.2）的解应当相差很小。

计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（1.2）中的扰动项改成18x ε或其它形式，实验中又有怎样的现象出现？
（3）（选作部分）请从理论上分析产生这一问题的根源。

注意我们可以将
方程（1.2）写成展开的形式， )
3.1(0
),(1920=+-= x x x p αα
同时将方程的解x 看成是系数α的函数，考察方程的某个解关于α的扰动是否敏感，与研究它关于α的导数的大小有何关系？为什么？你发现了什么现象，哪些根关于α的变化更敏感？
思考题一：（上述实验的改进）
在上述实验中我们会发现用roots 函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve 来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考MATLAB 的帮助。

思考题二：（二进制产生的误差）
用MATLAB 计算1001.01000
1-∑=i 。

结果居然有误差！因为从十进制数角度分析，
这一计算应该是准确的。

实验反映了计算机内部的二进制本质。

思考题三：（一个简单公式中产生巨大舍入误差的例子） 可以用下列式子计算自然对数的底数
n n n
e e )1
1(lim 1+==∞→
这个极限表明随着n 的增加，计算e 值的精度是不确定的。

现编程计算
n n
n f )1
1()(+=与exp(1)值的差。

n 大到什么程度的时候误差最大？你能解释其中
的原因吗？
poly(a) 求给定的根向量a 生成其对应的多项式系数（降序）向量 roots(p) 求解以向量p 为系数的多项式（降序）的所有根 poly2sym(p) 将多项式向量p 表示成为符号多项式（降序） sym(arg) 将数字、字符串或表达式arg 转换为符号对象 syms arg1 arg2 argk 将字符arg1,arg2,argk 定义为基本符号对象 solve('eq1') 求符号多项式方程eq1的符号解
实验二 非线性方程求根
实验2（迭代法、初始值与收敛性）
实验目的：初步认识非线性问题的迭代法与线性问题迭代法的差别，探讨迭代法及初始值与迭代收敛性的关系。

问题提出：迭代法是求解非线性方程的基本思想方法，与线性方程的情况一样，其构造方法可以有多种多样，但关键是怎样才能使迭代收敛且有较快的收敛速度。

实验内容：考虑一个简单的代数方程
012=--x x
针对上述方程，可以构造多种迭代法，如
)1.7(1
2
1-=+n n x x
)2.7(1
11n
n x x +
=+
)3.7(1
1+=+n n x x
在实轴上取初始值x 0，请分别用迭代（7.1）-（7.3）作实验，记录各算法的迭代过程。

实验要求：
（1）取定某个初始值，分别计算（7.1）-（7.3）迭代结果，它们的收敛性如何？重复选取不同的初始值，反复实验。

请自选设计一种比较形象的记录方式（如利用MATLAB 的图形功能），分析三种迭代法的收敛性与初值选取的关系。

（2）对三个迭代法中的某个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？
（3）线性方程组迭代法的收敛性是不依赖初始值选取的。

比较线性与非线性问题迭代的差异，有何结论和问题。

思考题一：用Newton 法求方程
013=--x x
在区间[-3,3]上误差不大于510-的根，分别取初值1,0,5.10-=x 进行计算，比较它们的迭代次数。

x=fzero(fun,x0) 返回一元函数fun 的一个零点，其中fun 为函数句柄，x0为标
量时，返回在x0附近的零点；x0为向量[a,b]时，返回函数在[a,b]中的零点
[x,f,h]=fsolve(fun,x0) 返回一元或多元函数x0附近fun 的一个零点，其中fun 为
函数句柄，x0为迭代初值；f 返回fun 在x 的函数值，应该接近0； h 返回值如果大于0，说明计算结果可靠，否则不可靠
实验三 解线性方程组的迭代法
实验3.1（病态的线性方程组的求解）
问题提出：理论的分析表明，求解病态的线性方程组是困难的。

实际情况是否如此，会出现怎样的现象呢？
实验内容：考虑方程组Hx=b 的求解，其中系数矩阵H 为Hilbert 矩阵，
n j i j i h h H j i n n j i ,,2,1,,
1
1
,
)(,, =-+=
=⨯
这是一个著名的病态问题。

通过首先给定解（例如取为各个分量均为1）再计算出右端b 的办法给出确定的问题。

实验要求：
（1）选择问题的维数为6，分别用Gauss 消去法、J 迭代法、GS 迭代法和SOR 迭代法求解方程组，其各自的结果如何？将计算结果与问题的解比较，结论如何？
（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？
思考题一：讨论病态问题求解的算法 Jacobi 迭代法与Gauss-seidel 迭代法的比较：
用Jacobi 迭代法与Gauss-seidel 迭代法解下列方程组：
（1）⎪⎪⎪⎭⎫
⎝⎛--=⎪⎪⎪⎪⎭
⎫ ⎝⎛
⎪⎪⎪⎭⎫ ⎝⎛--17752
32101110131
x x x （2）⎪⎪⎪⎭
⎫
⎝⎛-=⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛5.25.002
15.05.05.015.05.05.013
1
x x x 取6)0(10,)0,0,0(-==εT x ，你能得出什么结论？
思考题二：SOR （超松驰）迭代法松驰因子对收敛性及速度的影响 试用SOR （超松驰）迭代法求解下列方程组：
⎪⎪⎪⎭⎫
⎝⎛-=⎪⎪⎪⎪⎭
⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--243024241014303
431x x x 取6)0(10,)1,1,1(-==εT x ，选择松驰因子=ω0.8,0.9,1,1.1,1.2等，试看对算法收敛性的影响，并找出所选用的松驰因子的最佳者。

要求编制矩阵迭代求解的函数
文件。

实验3.2方程组性态讨论
（1） 求b Ax =的解向量，其中：
⎪⎪⎪⎭⎫
⎝
⎛=66
610/210/310/4100020003000
9
21A ⎪⎪⎪⎭
⎫ ⎝⎛=610/320001b
（2） 求系数矩阵A 的条件数；
（3） 将33a 改为610/3，3b 改为610/4，求解向量x
~; （4） 令)10,10,1(63-=diag P ，求解Pb PAx =，并求系数矩阵PA 的条件数；
（5） 对PA 中的33a 和Pb 中的3b 给以610-的扰动，求解向量x
ˆ 结合上述结果，讨论两个方程组的性态
实验3.2大型稀疏方程组的数值解法 设n 阶方阵：
⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭
⎫ ⎝
⎛--------
------=32
14
1213214
1412132
14
1412132141213
A B 为A 的各元素之和，显然Ax =b 的解为T x )1,...,1,1(=，用下面三种方法对于阶数n =100,200,…,500，误差限为63210,...,10,10---=ε的各种组合求解，分析收敛速度
（1） 选取列主元Gauss 消元法； （2） Jacobi 迭代法
（3） Gauss-seidel 迭代法 算法的改进：
由于是稀疏矩阵，请考虑改进矩阵的存储方式，减少存储空间的使用来提高计算效率。

也就是用稀疏存储方式。

实验四 解线性方程组的直接方法
实验4 （主元的选取与算法的稳定性）
问题提出：Gauss 消去法是我们在线性代数中已经熟悉的。

但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss 消去法作为数值算法的稳定性呢？Gauss 消去法从理论算法到数值算法，其关键是主元的选择。

主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验内容：考虑线性方程组
n n n R b R A b Ax ∈∈=⨯,,
编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss 消去过程。

实验要求：
（1）取矩阵⎥⎥
⎥
⎥⎥
⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ，则方程有解T x )1,,1,1(* =。

取n=10计算矩阵的条件数。

让程序自动选取主元，结果如何？
（2）现选择程序中手动选取主元的功能。

每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。

若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

（3）取矩阵阶数n=20或者更大，重复上述实验过程，观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异，说明主元素的选取在消去过程中的作用。

（4）将上述矩阵A 中的主元改为0.00006再重新作一次数值实验看看。

（5）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。

重复上述实验，观察记录并分析实验结果。

zeros(m,n) 生成m 行，n 列的零矩阵
ones(m,n) 生成m 行，n 列的元素全为1的矩阵 eye(n) 生成n 阶单位矩阵
rand(m,n) 生成m 行,n 列(0,1)上均匀分布的随机矩阵 diag(x) 返回由向量x 的元素构成的对角矩阵
tril(A) 提取矩阵A 的下三角部分生成下三角矩阵
triu(A) 提取矩阵A的上三角部分生成上三角矩阵rank(A) 返回矩阵A的秩
det(A) 返回方阵A的行列式
inv(A) 返回可逆方阵A的逆矩阵
[V，D]=eig(A) 返回方阵A的特征值和特征向量norm(A,p) 矩阵或向量的p范数
cond(A,p) 矩阵的条件数
[L，U，P]=lu(A) 选列主元LU分解
R=chol(X) 平方根分解
Hi=hilb(n) 生成n阶Hilbert矩阵。