《数值计算方法》

1《数值计算方法》上机实验指导书邹昌文编2009年10月2“数值计算方法”上机实验指导书实验一误差分析实验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”是两个互逆的运算函数。3))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大到什么程度的时候误差最大？你能解释其中的原因吗？4相关MATLAB函数提示：poly(a)求给定的根向量a生成其对应的多项式系数（降序）向量roots(p)求解以向量p为系数的多项式（降序）的所有根poly2sym(p)将多项式向量p表示成为符号多项式（降序）sym(arg)将数字、字符串或表达式arg转换为符号对象symsarg1arg2argk将字符arg1,arg2,argk定义为基本符号对象solve('eq1')求符号多项式方程eq1的符号解5实验二插值法实验2.1（多项式插值的振荡现象）问题提出：考虑一个固定的区间上用插值逼近一个函数。显然拉格朗日插值中使用的节点越多，插值多项式的次数就越高。我们自然关心插值多项式的次数增加时，)(xLn是否也更加靠近被逼近的函数。龙格（Runge）给出一个例子是极著名并富有启发性的。设区间[-1,1]上函数22511)(xxf实验内容：考虑区间[-1,1]的一个等距划分，分点为ninixi,,2,1,0,21则拉格朗日插值多项式为niijnxlxxL02)(2511)(其中的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中的函数或选择其他你有兴趣的函数，可以用MATLAB的函数“spline”作此函数的三次样条插值。实验要求：（1）随节点个数增加，比较被逼近函数和样条插值函数误差的变化情况。分6析所得结果并与拉格朗日多项式插值比较。（2）样条插值的思想是早产生于工业部门。作为工业应用的例子考虑如下问题：某汽车制造商用三次样条插值设计车门的曲线，其中一段的数据如下：xk012345678910yk0.00.791.532.192.713.033.272.893.063.193.29yk’0.80.2思考题一：（二维插值）在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程数据如下。试用MATLAB的二维插值函数“interp2”进行插值，并由此找出最高点和该点的高程。yx100200300400100636697624478200698712630478300680674598412400662626552334相关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）归纳总结数值实验结果，试定性地说明函数逼近各种方法的适用范围，及实际应用中选择方法应注意的问题。思考题一：（病态）考虑将[0,1]30等分节点，用多项式51xxy生成数据，再用polyfit求其3次、5次、10次、15次拟合多项式，并分析误差产生的原因。相关MATLAB函数提示：p=polyfit(x,y,k)用k次多项式拟合向量数据(x,y)，返回多项式的降幂系数，当k=n-1时，实现多项式插值，这里n是向量维数8实验四数值积分与数值微分实验4.1（高斯数值积分方法用于积分方程求解）问题提出：线性的积分方程的数值求解，可以被转化为线性代数方程组的求解问题。而线性代数方程组所含未知数的个数，与用来离散积分的数值方法的节点个数相同。在节点数相同的前提下，高斯数值积分方法有较高的代数精度，用它通常会得到较好的结果。实验内容：求解第二类Fredholm积分方程btatfdssystktyba),()(),()(首先将积分区间[a,b]等分成n份，在每个子区间上离散方程中的积分就得到线性代数方程组。实验要求：分别使用如下方法，离散积分方程中的积分1.复化梯形方法；2.复化辛甫森方法；3.复化高斯方法。求解如下的积分方程。1）ttedssyeety10)(12)(，方程的准确解为te；2）ttsteedssyety210211)()(，方程的准确解为te；比较各算法的计算量和误差以分析它们的优劣。实验4.2（高维积分数值计算的蒙特卡罗方法）问题提出：高维空间中的积分，如果维数不很高且积分区域是规则的或者能等价地写成多重积分的形式，可以用一元函数积分的数值方法来计算高维空间的积分。蒙特卡罗方法对计算复杂区域甚至不连通的区域上的积分并没有特殊的困难。实验内容：对于一般的区域，计算其测度（只要理解为平面上的面积或空间中的体积）的一般方法是：先找一个规则的区域A包含，且A的测度是已知的。生成区域A中m个均匀分布的随机点mipi,,2,1,，如果其中有n个落在区域中，则区域的测度m()为n/m。函数f(x)在区域上的积分可以近似为：区域的测度与函数f(x)在中n个随机点上平均值的乘积，即kpkpfnmdxxf)(1)()(实验要求：假设冰琪淋的下部为一锥体而上面为一半球，考虑冰琪淋体积问题：计算锥面222yxz上方和球面1)1(222zyx内部区域的体积。如果使用球面坐标，该区域可以表示为如下的积分：4/0cos20202sinddd用蒙特卡罗方法可以计算该积分。另一方面，显然这样的冰琪淋可以装在如下立方体的盒子里20,11,11zyx9而该立方体的体积为8。只要生成这个盒子里均匀分布的随机点，落入冰琪淋锥点的个数与总点数之比再乘以8就是冰琪淋锥的体积。比较两种方法所得到的结果。类似的办法可以计算复杂区域的测度（面积或体积）。试求由下列关系所界定区域的测度：0)sin(31,21,10)1(yzyezyxx22941,31)2(33xeyyxyx1sin10,10,10)3(2xezxzyxzyx相关MATLAB函数提示：diff(x)如果x是向量，返回向量x的差分；如果x是矩阵，则按各列作差分diff(x,k)k阶差分q=polyder(p)求得由向量p表示的多项式导函数的向量表示qFx=gradient(F,x)返回向量F表示的