数值分析上机实验指导书

“数值计算方法”上机实验指导书实验一误差分析实验1.1（病态问题）实验目的：算法有“优”与“劣”之分，问题也有“好”与“坏”之别。对数值方法的研究而言，所谓坏问题就是问题本身对扰动敏感者，反之属于好问题。通过本实验可获得一个初步体会。数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。问题提出：考虑一个高次的代数多项式)1.1()()20()2)(1()(201kkxxxxxp显然该多项式的全部根为1,2,…,20共计20个，且每个根都是单重的。现考虑该多项式的一个扰动)2.1(0)(19xxp其中是一个非常小的数。这相当于是对（1.1）中19x的系数作一个小的扰动。我们希望比较（1.1）和（1.2）根的差别，从而分析方程（1.1）的解对扰动的敏感性。实验内容：为了实现方便，我们先介绍两个MATLAB函数：“roots”和“poly”。roots(a)u其中若变量a存储n+1维的向量，则该函数的输出u为一个n维的向量。设a的元素依次为121,,,naaa，则输出u的各分量是多项式方程01121nnnnaxaxaxa的全部根；而函数poly(v)b的输出b是一个n+1维向量，它是以n维向量v的各分量为根的多项式的系数。可见“roots”和“poly”是两个互逆的运算函数。))20:1((;)2();21,1(;000000001.0vepolyrootsessvezerosveess上述简单的MATLAB程序便得到（1.2）的全部根，程序中的“ess”即是（1.2）中的。实验要求：（1）选择充分小的ess，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数很小，我们自然感觉（1.1）和（1.2）的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？（2）将方程（1.2）中的扰动项改成18x或其它形式，实验中又有怎样的现象出现？（3）（选作部分）请从理论上分析产生这一问题的根源。注意我们可以将方程（1.2）写成展开的形式，)3.1(0),(1920xxxp同时将方程的解x看成是系数的函数，考察方程的某个解关于的扰动是否敏感，与研究它关于的导数的大小有何关系？为什么？你发现了什么现象，哪些根关于的变化更敏感？思考题一：（上述实验的改进）在上述实验中我们会发现用roots函数求解多项式方程的精度不高，为此你可以考虑用符号函数solve来提高解的精确度，这需要用到将多项式转换为符号多项式的函数poly2sym,函数的具体使用方法可参考MATLAB的帮助。思考题二：（二进制产生的误差）用MATLAB计算1001.010001i。结果居然有误差！因为从十进制数角度分析，这一计算应该是准确的。实验反映了计算机内部的二进制本质。思考题三：（一个简单公式中产生巨大舍入误差的例子）可以用下列式子计算自然对数的底数nnnee)11(lim1这个极限表明随着n的增加，计算e值的精度是不确定的。现编程计算nnnf)11()(与exp(1)值的差。n大到什么程度的时候误差最大？你能解释其中的原因吗？相关MATLAB函数提示：poly(a)求给定的根向量a生成其对应的多项式系数（降序）向量roots(p)求解以向量p为系数的多项式（降序）的所有根poly2sym(p)将多项式向量p表示成为符号多项式（降序）sym(arg)将数字、字符串或表达式arg转换为符号对象symsarg1arg2argk将字符arg1,arg2,argk定义为基本符号对象solve('eq1')求符号多项式方程eq1的符号解实验二插值法实验2.1（多项式插值的振荡现象）问题提出：考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，)(xLn是否也更加靠近被逼近的函数。龙格（Runge）给出一个例子是极著名并富有启发性的。设区间[-1,1]上函数22511)(xxf实验内容：考虑区间[-1,1]的一个等距划分，分点为ninixi,,2,1,0,21则拉格朗日插值多项式为201()()125nniiiLxlxx其中的nixli,,2,1,0),(是n次拉格朗日插值基函数。实验要求：（1）选择不断增大的分点数目n=2,3….，画出原函数f(x)及插值多项式函数)(xLn在[-1，1]上的图像，比较并分析实验结果。（2）选择其他的函数，例如定义在区间[-5，5]上的函数xxgxxxharctan)(,1)(4重复上述的实验看其结果如何。（3）区间[a,b]上切比雪夫点的定义为1,,2,1,)1(2)12(cos22nknkababxk以121,,nxxx为插值节点构造上述各函数的拉格朗日插值多项式，比较其结果,试分析原因。实验2.2（样条插值的收敛性）问题提出：多项式插值是不收敛的，即插值的节点多，效果不一定就好。对样条函数插值又如何呢？理论上证明样条插值的收敛性是比较困难的，但通过本实验可以验证这一理论结果。实验内容：请按一定的规则分别选择等距或者非等距的插值节点，并不断增加插值节点的个数。考虑实验2.1中的函数或选择其他你有兴趣的函数。实验要求：（1）随节点个数增加，比较被逼近函数和样条插值函数误差的变化情况。分析所得结果并与拉格朗日多项式插值比较(可以用MATLAB的函数“spline”作此函数的三次样条插值，取n=10、20，分别画出插值函数及原函数的图形)。（2）样条插值的思想是早产生于工业部门。作为工业应用的例子考虑如下问题：某汽车制造商用三次样条插值设计车门的曲线，其中一段的数据如下：xk012345678910yk0.00.791.532.192.713.033.272.893.063.193.29yk’0.80.2要求：i.自己编程计算（用三弯矩、三转角方程均可）ii.主函数myspline(x,y,边界类型，边界值，xi)其中：x节点y节点上的函数值xi未知节点返回：S(xi)iii.三对角方程组用追赶法求解(书P160)。实验2.3：(一维插值的应用—画图)画你自己的手的形状，在MATLAB中输入figure('position',get(0,'screensize'))axes('position',[0011])[x,y]=ginput;将你的手掌张开放在计算机屏幕上，然后使用计算机鼠标选取一系列点勾勒出手的轮廓，按回车键结束ginput过程，这样就获得了一系列你的手掌外形数据点),(yx。也可以这样获得数据点),(yx，先把手放在一张白纸上，并用笔画出它的轮廓线，然后将纸贴在计算机屏幕上，透过纸能看到平面上的鼠标，并通过ginput记录下轮廓上的点。将x和y坐标值看作是两个独立变量的函数，独立变量的取值为从1到记录的点的数目。利用MATLAB的插值函数进行插值，并画出你的手掌外形轮廓。6思考题：（二维插值）1.在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程数据如下。试用MATLAB的二维插值函数“interp2”进行插值，并由此找出最高点和该点的高程。yx1002003004001006366976244782006987126304783006806745984124006626265523342.画出山区地貌图。利用MATLAB的peaks函数生成某山区的一些地点及其高度三维数据（单位：m）。命令格式：[x,y,z]=peaks(n)，生成的n阶矩阵x,y,z为测量的山区地点三维数据（n=30）。根据peaks函数生成的数据，利用Matlab二维插值画出该山区的地貌图和等值线图(提示函数：interp2、meshgrid、plot3等)。相关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的分段线性插值结果yipp=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为网格向量的三角形线性插值（不规则数据的二维插值）7实验三函数逼近与曲线拟合实验3.1（曲线逼近方法的比较）问题提出：曲线的拟合和插值，是逼近函数的基本方法，每种方法具有各自的特点和特定的适用范围，实际工作中合理选择方法是重要的。实验内容：考虑实验2.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’);holdon;p3=polyfit(x,y,3);yy=polyval(p3,xx);plot(x,y,’o’,xx,yy);holdoff;适当修改上述MATLAB程序，也可以拟合其他你有兴趣的函数。实验要求：（1）将拟合的结果与拉格朗日插值及样条插值的结果比较。（2）归纳总结数值实验结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中选择方法应注意的问题。实验3.2：（最小二乘拟合的经验公式和模型）1.（已知经验公式）：某类疾病发病率为y‰和年龄段x(每五年为一段，例如0~5岁为第一段，6~10岁为第二段……)之间有形如bxaey的经验关系，观测得到的数据表如下x123456789y0.8982.383.071.842.021.942.222.774.02x10111213141516171819y4.765.466.5310.916.522.535.750.661.681.8实验要求：（1）用最小二乘法确定模型bxaey中的参数a和b(提示函数：lsqcurvefit，lsqnonlin)。（2）利用MATLAB画出离散数据及拟合函数bxaey图形。8（3）利用MATLAB画出离散点处的误差图，并计算相应的均方误差。2．（最小二乘拟合模型未知）某年美国轿车价格的调查资料如表，其中ix表示轿车的使用年数，iy表示相应的平均价格，实验要求：试分析用什么形式的曲线来拟合表中的数据，并预测使用4.5年后轿车的平均价格大致为多少？ix12345678910iy2615194314941087765538484290226204实验3.3（研究最佳平方逼近多项式的收敛性质）实验内容和要求：取函数xexf)(，在[-1，1]上以勒让德多项式为基函数，对于10,,1,0n构造最佳平方逼近多项式)(xpn，令)()()(xpxfxnn，将xxn~)(的曲线画在一个图上。令11max()()nnxEfxpx，画出~nEn的曲线。做出~nEn之间的最小二乘曲线，能否提出关于收敛性的猜测。思考题一：（病态）考虑将[0,1]30等分节点，用多项式51xxy生成数据，再用polyfit求其3次、5次、10次、15次拟合多项式，并分析误差产生的原因。思考题二：调研Matlab拟合工具箱的使用。相关MATLAB函数提示：p=polyfit(