数值分析实验指导

数值分析实验指导

2011年8月

实验一 误差分析

实验1.1（病态问题）

实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。

数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。

问题提出：考虑一个高次的代数多项式

20p(x)(x1)(x2)(x20)(xk)k1(1.1)

显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动

p(x)x190(1.2)

其中是一个非常小的数。这相当于是对（1.1）中x19的系数作一个小的扰动。我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。

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

uroots(a)

其中若变量a存储n+1维的向量，则该函数的输出u为一个n维的向量。设a的元素依次为a1,a2,,an1，则输出u的各分量是多项式方程

a1xa2xnn1anxan10

的全部根；而函数 bpoly(v )的输出b是一个n+1维向量，它是以n维向量v的各分量为根的多项式的系数。可见“roots”和“poly”是两个互逆的运算函数。

ess0.000000001;vezeros(1,21);ve(2)ess;roots(poly(1:20)ve)

上述简单的MATLAB程序便得到（1.2）的全部根，程序中的“ess”即是（1.2）中的。

实验要求：

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

如果扰动项的系数很小，我们自然感觉（1.1）和（1.2）的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（1.2）中的扰动项改成x18或其它形式，实验中又有怎样的现象出现？

（3）（选作部分）请从理论上分析产生这一问题的根源。注意我们可以将

方程（1.2）写成展开的形式， p(x,)x20x190(1.3)

同时将方程的解x看成是系数的函数，考察方程的某个解关于的扰动是否敏感，与研究它关于的导数的大小有何关系？为什么？你发现了什么现象，哪些根关于的变化更敏感？

思考题一：（上述实验的改进）

在上述实验中我们会发现用roots函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考MATLAB的帮助。

思考题二：（二进制产生的误差）

1000用MATLAB计算0.1100。结果居然有误差！因为从十进制数角度分析，

i1这一计算应该是准确的。实验反映了计算机内部的二进制本质。 思考题三：（一个简单公式中产生巨大舍入误差的例子） 可以用下列式子计算自然对数的底数

ee1lim(1)n

n1n这个极限表明随着n的增加，计算e值的精度是不确定的。现编程计算

f(n)(11n)n与exp(1)值的差。n大到什么程度的时候误差最大？你能解释其中

的原因吗？

相关MATLAB函数提示： poly(a) 求给定的根向量a生成其对应的多项式系数（降序）向量 roots(p) 求解以向量p为系数的多项式（降序）的所有根 poly2sym(p) 将多项式向量p表示成为符号多项式（降序） sym(arg) 将数字、字符串或表达式arg转换为符号对象 syms arg1 arg2 argk 将字符arg1,arg2,argk定义为基本符号对象 solve('eq1') 求符号多项式方程eq1的符号解 实验1.2 误差传播与算法稳定性

实验目的：体会稳定性在选择算法中的地位。误差扩张的算法是不稳定的，是我们所不期望的；误差衰减的算法是稳定的，是我们努力寻求的，这是贯穿本课程的目标。

问题提出：考虑一个简单的由积分定义的序列

In10xenx1dx,n1,2...xex1

。而对于n2时，利用分

显然In0,n1,2,...。当n1时，部积分易得

In1nx1nI1110dx1/e0xedxxex1x1|0nx0n1n1ex1dx1nIn1,n2,3,...

另一方面，我们有

In10xendx10xdx1/(n1)

实验内容：由以上递推关系，我们可得到计算序列{In}的两种方法。

., （I）：I11/e,In1nIn1,n2.,.3（II）：

EN0,En11Enn,nN,N1,N.2..,3,2

实验要求：

（1）分别用算法（I）、（II）并在计算中分别采用5位、6位和7位有效数字，请判断哪种算法能给出更精确的结果。

（2）两种算法的优劣，与你的第一感觉是否吻合。请从理论上证明你实验得出的结果，解释实验的结果。算法（I）中的I1的计算误差为e1，由I1递推计算In的误差为en；算法（II）中IN的计算误差为N，由IN向前递推计算In(nN)的误差为n。如果在上述两算法中都假定后面的计算不再引入其他误差， 试给出en与e1的关系和n与N的关系。

（3）算法（I）中通常e1会很小，当n增大时，en的变化趋势如何？算法（II）中N通常相对比较大，当n减小时，误差n又是如何传播的？也就是说比较一下上述两个算法，当某一步产生误差后，该误差对后面的影响是衰减还是扩张的。 （4）通过理论分析与计算实验，针对算法（I）和（II）的稳定性，给出你的结论。

实验二 插值法

实验2.1（多项式插值的振荡现象）

问题提出：考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，Ln(x)是否也更加靠近被逼近的函数。龙格（Runge）给出一个例子是极著名并富有启发性的。设区间[-1,1]上函数

f(x)1125x2in2

实验内容：考虑区间[-1,1]的一个等距划分，分点为 xi1,i0,1,2,,n

则拉格朗日插值多项式为

n Ln(x)125xi012jli(x)

其中的li(x),i0,1,2,,n是n次拉格朗日插值基函数。

实验要求：

（1） 选择不断增大的分点数目n=2,3….，画出原函数f(x)及插值多项式函数Ln(x)在[-1，1]上的图像，比较并分析实验结果。

（2）选择其他的函数，例如定义在区间[-5，5]上的函数

h(x)x1x4,g(x)arctanx

重复上述的实验看其结果如何。

（3）区间[a,b]上切比雪夫点的定义为 xkba2ba2(2k1)cos2(n1),k1,2,,n1 以x1,x2,xn1为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果。

实验2.2（样条插值的收敛性）

问题提出：多项式插值是不收敛的，即插值的节点多，效果不一定就好。对样条函数插值又如何呢？理论上证明样条插值的收敛性是比较困难的，但通过本实验可以验证这一理论结果。

实验内容：请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。考虑实验2.1中的函数或选择其他你有兴趣的函数，可以用MATLAB的函数“spline”作此函数的三次样条插值。

实验要求：

（1）随节点个数增加，比较被逼近函数和样条插值函数误差的变化情况。分析所得结果并与拉格朗日多项式插值比较。

（2）样条插值的思想是早产生于工业部门。作为工业应用的例子考虑如下问题：某汽车制造商用三次样条插值设计车门的曲线，其中一段的数据如下： xk 0 1 2 3 4 5 6 7 8 9 10 yk 0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29 yk’ 0.8 0.2

思考题一：（二维插值） 在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程数据如下。试用MATLAB的二维插值函数“interp2”进行插值，并由此找出最高点和该点的高程。 y x 100 200 300 400 100 636 697 624 478 200 698 712 630 478 300 680 674 598 412 400 662 626 552 334

相关MATLAB函数提示： plot(x,y) 作出以数据(x(i),y(i))为节点的折线图，其中x,y为同长度的向量 subplot(m,n,k) 将图形窗口分为m*n个子图，并指向第k幅图 yi=interp1(x,y,xi) 根据数据(x,y)给出在xi的分段线性插值结果yi pp=spline(x,y) 返回样条插值的分段多项式(pp)形式结构 pp=csape(x,y,‘边界类型’,‘边界值’) 生成各种边界条件的三次样条插值 yi=ppval(pp,xi) pp样条在xi的函数值 ZI=interp2(x,y,z,xi,yi) x,xi为行向量，y,yi为列向量，z为矩阵的双线性二维插值 ZI=interp2(…,'spline') 使用二元三次样条插值 ZI=griddata(x,y,z,xi,yi) x,y,z均为向量（不必单调）表示数据，xi,yi为网格向量的三角形线性插值（不规则数据的二维插值）

实验三 曲线拟合

实验3.1

n编制以函数{xk}k为基的多项式最小二乘拟合程序，并用于对下表中的数0据做3次多项式最小二乘拟合。 xi yi -1.0 -4.447 -0.5 -0.452 0.0 0.551 n*0.5 0.048 1.0 -0.447 1.5 0.549 22.0 4.552 取权数wi1,求拟合曲线*k并作离kx中的参数{k}、平方误差，k0散数据(xi,yi)的拟合函数y*(x)的图形。

实验3.2

编制正交化多项式最小二乘拟合程序，并用于求解上题中3次多项式最小二乘拟合问题，作拟合曲线的图形，计算平方误差，并于上题结果进行比较。

实验3.3（曲线逼近方法的比较）

问题提出：曲线的拟合和插值，是逼近函数的基本方法，每种方法具有各自的特点和特定的适用范围，实际工作中合理选择方法是重要的。

实验内容：考虑实验5.1中的著名问题。下面的MATLAB程序给出了该函数的二次和三次拟合多项式。

x=-1:0.2:1;

y=1/(1+25*x.*x); xx=-1:0.02:1; p2=polyfit(x,y,2); yy=polyval(p2,xx); plot(x,y,’o’,xx,yy); xlabel(‘x’); ylabel(‘y’); hold on;

p3=polyfit(x,y,3); yy=polyval(p3,xx); plot(x,y,’o’,xx,yy); hold off;

适当修改上述MATLAB程序，也可以拟合其他你有兴趣的函数。 实验要求：

（1）将拟合的结果与拉格朗日插值及样条插值的结果比较。 （2）归纳总结数值实验结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中选择方法应注意的问题。

思考题一：（病态）考虑将[0,1]30等分节点，用多项式y1xx5生成数据，再用polyfit求其3次、5次、10次、15次拟合多项式，并分析误差产生的原因。

相关MATLAB函数提示： p=polyfit(x,y,k) 用k次多项式拟合向量数据(x,y)，返回多项式的降幂系数，当k>=n-1时，实现多项式插值，这里n是向量维数

实验四 数值积分

数值实验综述：通过数值积分实验掌握数值积分的实现，理解各种数值积分公式的特性，并能用数值积分求解积分方程和微分方程。 基础实验 实验4.1

实验目的：复化求积公式计算定积分。

实验题目：计算下列各式右端定积分的近似值

(1) ln2ln3223221x1dx (2)

4111x20dx

(3) ln3103dxx (4)

e221xedxx

试验要求：

（l）若用复化梯形公式、复化Simpson公式和复化Gauss-Legendre I型公式作计算，要求绝对误差限为

12107，分别利用它们的余项对每种算法作出步长

的事前估计。

（2）分别用复化梯形公式、复化Sboson公式和复化Gauss-Legendre I型公式作计算。

（3）将计算结果与精确解作比较，并比较各种算法的计算量。

实验4.2 比较复化Simpson方法和变步长（区间逐次分半求积法）Simpson方法 试验题：计算下列定积分

(1)

20x621dxxx10xxdx (2)0 (3)

2001x5dx

试验要求：

（1）分别用复化Simpson公式、变步长Simpson公式和自适应Simpson公式计算，要求

12107绝对误差限为，输出每种方法所需的节点数和积分近似值。

（2）分析比较计算结果。

相关MATLAB函数提示： diff(x) 如果x是向量，返回向量x的差分；如果x是矩阵，则按各列作差分 diff(x,k) k阶差分 q=polyder(p) 求得由向量p表示的多项式导函数的向量表示q Fx=gradient(F,x) 返回向量F表示的一元函数沿x方向的导函数F'(x)，其中x是与F同维数的向量 z=trapz(x,y) x表示积分区间的离散化向量；y是与x同维数的向量，表示被积函数；z返回积分的近似值 z=guad(fun,a,b,tol) 自适应步长Simpson积分法求得Fun在区间[a,b]上的定积分，Fun为M文件函数句柄，tol为积分精度 z=dblquad(fun,a,b,c,d,tol,method) 求得二元函数Fun(x,y)的重积分 z=triplequad(fun,a,b,c,d,e,f,tol,method) 求得三元函数Fun(x,y,z)的重积分

实验五 解线性方程组的直接方法

实验5.1 （主元的选取与算法的稳定性）

问题提出：Gauss消去法是我们在线性代数中已经熟悉的。但由于计算机的数值运算是在一个有限的浮点数集合上进行的，如何才能确保Gauss消去法作为数值算法的稳定性呢？Gauss消去法从理论算法到数值算法，其关键是主元的选择。主元的选择从数学理论上看起来平凡，它却是数值分析中十分典型的问题。

实验内容：考虑线性方程组

Axb,ARnn,bRn

编制一个能自动选取主元，又能手动选取主元的求解线性方程组的Gauss消去过程。

实验要求：

68（1）取矩阵A161868715,b，则方程有解x*(1,1,,1)T。115614取n=10计算矩阵的条件数。让程序自动选取主元，结果如何？

（2）现选择程序中手动选取主元的功能。每步消去过程总选取按模最小或按模尽可能小的元素作为主元，观察并记录计算结果。若每步消去过程总选取按模最大的元素作为主元，结果又如何？分析实验的结果。

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

（4）将上述矩阵A中的主元改为0.00006再重新作一次数值实验看看。 （5）选取其他你感兴趣的问题或者随机生成矩阵，计算其条件数。重复上述实验，观察记录并分析实验结果。

实验5.2（线性代数方程组的性态与条件数的估计） 问题提出：理论上，线性代数方程组Axb的摄动满足

xxcond(A)1A1AbAbA

矩阵的条件数确实是对矩阵病态性的刻画，但在实际应用中直接计算它显然不现实，因为计算A1通常要比求解方程Axb还困难。

实验内容：MATLAB中提供有函数“condest”可以用来估计矩阵的条件数，

它给出的是按1-范数的条件数。首先构造非奇异矩阵A和右端，使得方程是可以精确求解的。再人为地引进系数矩阵和右端的摄动A和b，使得A和b充分小。

实验要求：

（1）假设方程Ax=b的解为x，求解方程(AA)xˆbb，以1-范数，给出

xxˆxxx的计算结果。

（2）选择一系列维数递增的矩阵（可以是随机生成的），比较函数“condest”所需机器时间的差别.考虑若干逆是已知的矩阵，借助函数“eig”很容易给出cond2(A)的数值。将它与函数“cond(A,2)”所得到的结果进行比较。

（3）利用“condest”给出矩阵A条件数的估计，针对（1）中的结果给出

xx的理论估计，并将它与（1）给出的计算结果进行比较，分析所得结果。注意，如果给出了cond(A)和A的估计，马上就可以给出A1的估计。 （4）估计著名的Hilbert矩阵的条件数。

H(hi,j)nn,hi,j1ij1,i,j1,2,,n

思考题一：（Vadermonde矩阵）设

22211 A11x0x1x2xnx0x1x2xn2nix0nin0x0inx1x1i0nx2，bni，

x2i0nxnnixni0其中，xk10.1k,k0,1,,n，

（1）对n=2,5,8，计算A的条件数；随n增大，矩阵性态如何变化？

（2）对n=5，解方程组Ax=b；设A的最后一个元素有扰动10-4，再求解Ax=b （3）计算（2）扰动相对误差与解的相对偏差，分析它们与条件数的关系。 （4）你能由此解释为什么不用插值函数存在定理直接求插值函数而要用拉格朗

日或牛顿插值法的原因吗？

相关MATLAB函数提示： 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矩阵

实验六 解线性方程组的迭代法

实验6.1 迭代法收敛速度试验 试验目的：认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。

试验题目：用迭代法求解方程组Axb，其中AR2020，它的每条对角钱元素是常 数，为

31/21/4A1/231/21/41/231/41/41/21/21/4(0)31/21/41/23

试验要求：（1）选取不同的初始向量x和不同的方程组右端项向量b，给定迭

代误差要求，用Jacobi迭代法和Gauss-Seidel迭代法计算，观测得到的迭代向量序列是否均收敛？若收敛；记录迭代次数，分析计算结果并得出你的结论； （2）取定右端问量b和初始向量x(0)，将A的主对角线元素成倍增长若干次，

非主对角线元素不变，每次用Jacobo迭代法计算，要求迭代误差满足

||x(k1)x(k)||105比较收敛速度，分析现象并得出你的结论。

实验6.2（病态的线性方程组的求解）

问题提出：理论的分析表明，求解病态的线性方程组是困难的，实际情况是否如此，会出现怎样的现象呢？

实验内容：考虑方程组Hxb的求解，其中系数矩阵H为Hilbert矩阵，

H(hij)nn,hij1ij1,i,j1,2,...,n

这是一个著名的病态问题．通过首先给定解（例如取为各分量均为1）再计算出右端的办法给出确定的问题． 实验要求：

ibocaJ（1）选择问题的维数为6，分别用Gauss消去法（即LU分解），

迭

代方法、GS迭代和 SOR迭代求解方程组，其各自的的结果如何？将计算结果

与问题的解比较，结论如何。

（2）逐步增大问题的维数，仍然用上述的方法来解它们，计算的结果如何？计算的结果说明了什么？

（3）讨论病态问题求解的算法．

实验七 非线性方程求根

数值实验综述：非线性方程（组）的数值解法是计算数学的一个重要课题，在实际问题中有广泛的应用，特别在各种非线性问题的科学计算中更显出它的重要性。本章讨论非线性方程和非线性方程组的数值解法，主要讲授了解非线性方程的二分法、迭代法、迭代加速、Newton法割线法和抛物线法等；同时介绍了解非线性方程组的Newton法及其变形和拟Newton法。本章数值实验的中心就是通过具体的数值实验来研究算法的可行性，如算法的收敛性，初值的选择，收敛速度以及算法的复杂性与稳定性等。 基础实验

实验7.1 迭代法、初始值与收敛性

实验目的：初步认识非线性问题的迭代法与线性问题迭代法的差别，探讨迭代法及初始值与迭代收敛性的关系。

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

实验内容：考虑一个简单的代数方程

xx10，

xn111xn22针对上述方程，可以构造多种迭代法，如xn1xn1,,xn1xn1等。在实轴上取初始值x0，分别用以上迭代作实验，记录各算法的迭代过程。 实验要求：

（1）取定某个初始值，按如上迭代格式进行计算，它们的收敛性如何？重复选取不同的初始值，反复实验。请读者自行设计一种比较形象的记录方式（如利用Matlab的图形功能），分析三种迭代法的收敛性与初值选取的关系。

（2）对三个迭代法中的某一个，取不同的初始值进行迭代，结果如何？试分析迭代法对不同的初值是否有差异？

（3）线性方程组迭代法的收敛性是不依赖初始值选取的。比较线性与非线性问题迭代的差异，有何结论和问题。 实验7.2 算法设计与比较

问题提出：非线性方程组的数值求解方法很多，基本的思想是线性化。不同的方法效果如何，要靠计算的实践来分析、比较。 实验内容：考虑算法：（l）牛顿法；（2）简化牛顿法； 分别编写它们的Matlab程序。 实验要求：（1）用上述各种方法，分别计算下面的两个例子。在达到精度相同的前提下，比较其迭代次数、浮点运算次数和CPU时间等。

12xx2124x370x2110x2x3110,取x(0)(0,0,0)T（I）x3210x380

3x1cos(x12x3)02x2x2181(20.1)sinx31.060取x(0)(0,0,0)Tex1x220x1(（II）33103)0

（2）取其它的初值x(0)，结果有如何？反复选取不同的初值，比较其结果。（3）总结归纳你的实验结果，试说明各种方法适用的问题。