数值分析课程设计.doc

课程设计数值分析专业班级： 信息与计算科学08-1班学 号： 20085975姓 名： 李 才 芳实验1.1水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只给猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？算法分析：求解这一问题可以用迭代递推算法。首先分析椰子数目的变化规律，设最初的椰子数为p 0，即第一个水手所处理之前的椰子数，用p 1、p 2、p 3、p4、p 5 分别表示五个水手对椰子动了手脚以后剩余的椰子数目，则根据问题有 再用x表示最后每个水手平分得到的椰子数，于是有 所以 利用逆向递推的方法，有 有了逆向递推关系式，求解这一问题似乎很简单，但由于椰子数为一正整数，用任意的x作为初值递推出的p0数据不一定是合适的。 这里用 for 循环语句结合 break 语句来寻找合适的 x 和 p0 ，对任意的 x 递推计算出 p0 ，当计算结果为正整数时，结果正确，否则选取另外的x 再次重新递推计算，直到计算出的结果 p0 为正整数为止。MATLAB编程代码：(1) n=input(输出n的值：); for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end if p=fix(p) break endenddisp(x,p)输出结果：输出n的值：15001023 15621(2) for x=1:inf p=5*x+1; for k=1:5 p=5*p/4+1; end if p=fix(p) break end enddisp(x,p)输出结果：1023 15621结果分析：在解决本题的过程中，运用了迭代的方法，每步还要判断X是否能被4整除，从而试探出结果C语言编程代码：#include int count(int); void main() int i, n, y; printf( 输入水手数：n ); scanf( %d ,&n); y=count(n); for(i=0;i n;i+) printf( %dn ,y); y=(y-1)/5*4; int count(int a) int m,i,k=1,ok=0; for(i=1;) if(i=1) m=k; if(k*5+1)%4=0) k=(k*5+1)/4; i+; else k=+m; i=1; if(i=a&ok eval(ans)ans =0.1823则取I0为0.18syms x n; int(x30/(5+x),0,1) ans =931322574615478515625*log(2)+931322574615478515625*log(3)-931322574615478515625*log(5)-79095966183067699902965545527073/465817912560 eval(ans)ans = 0则由于MATLAB中小数点之后只保留四位,而I的积分值又不为0,故I30取值为0.00001-0.0001之间的数编程代码：:function I=f(I) if I=0.1& I0& I=0.0001 for n=30:-1:2 I= (-1/5)*I+1/(5*n) end else disp(请重新输入合适的I值) endend输出结果：shiyan1_2(0.18)I = 0.1000I = -4.4409e-016I = 0.3333I = -1.4167I = 7.2833I = -36.2500I = 181.3929I = -906.8393I = 4.5343e+003I = -2.2671e+004I = 1.1336e+005I = -5.6679e+005I = 2.8339e+006I = -1.4170e+007I = 7.0848e+007I = -3.5424e+008I = 1.7712e+009I = -8.8560e+009I = 4.4280e+010I = -2.2140e+011shiyan1_2f(0.00005)I = 0.0067I = 0.0056I = 0.0060I = 0.0062I = 0.0065I = 0.0067I = 0.0070I = 0.0073I = 0.0076I = 0.0080I = 0.0084I = 0.0088I = 0.0093I = 0.0099I = 0.0105I = 0.0112I = 0.0120I = 0.0130I = 0.0141I = 0.0154I = 0.0169I = 0.0188I = 0.0212I = 0.0243I = 0.0285I = 0.0343I = 0.0431I = 0.0580I =0.0884结果分析：第二种方法所得的结果相对来说比较精确一些，也比较可靠。因为第一种方法每经过一次迭代都将最初的误差放大了五倍，使得最终的误差越来越大；而第二种方法每经过一次迭代就将误差缩小为初始误差的五分之一，使得最终的误差越来越小，因此相对来说比较可靠，性能较好。实验13 绘制Koch分形曲线问题描述：从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形（图1）；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形（图2），这时，图形中共有17个结点。这种迭代继续进行下去可以形成Koch分形曲线。在迭代过程中，图形中的结点将越来越多，而曲线最终显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。Koch分形曲线的绘制与算法设计和计算机实现相关。问题分析：考虑由直线段（2个点）产生第一个图形（5个点）的过程，设和分别为原始直线段的两个端点。现在需要在直线段的中间依次插入三个点产生第一次迭代的图形（图1）。显然，位于点右端直线段的三分之一处，点绕旋转60度（逆时针方向）而得到的，故可以处理为向量经正交变换而得到向量，形成算法如下：（1）；（2）；（3）；在算法的第三步中，A为正交矩阵。；这一算法将根据初始数据（和点的坐标），产生图1中5个结点的坐标。这5个结点的坐标数组，组成一个52矩阵。这一矩阵的第一行为为的坐标，第二行为的坐标，第二行为的坐标第五行为的坐标。矩阵的第一列元素分别为5个结点的x坐标 ，第二列元素分别为5个结点的y坐标。问题思考与实验：（1）考虑在Koch分形曲线的形成过程中结点数目的变化规律。设第k次迭代产生结点数为，第迭代产生结点数为，试写出和之间的递推关系式；（2）参考问题分析中的算法，考虑图1到图2的过程，即由第一次迭代的5个结点的结点坐标数组，产生第二次迭代的17个结点的结点坐标数组的算法；（3）考虑由第k次迭代的个结点的结点坐标数组，产生第次迭代的个结点的结点坐标数组的算法；（4）设计算法用计算机绘制出如下的Koch分形曲线（图3）。算法分析：（1）第k次迭代产生的结点数为 ,第k+1 次迭代产生的结点数为 (2) 第一次迭代的5个结点的结点坐标数组， 在和 之间（1） = + ()/3（2） = + 2()/3（3） = + ()在和 之间（4） = + ()/3（5） = + 2()/3（6） = + ()在和 之间（7） = + ()/3（8） = + 2()/3（9） = + ()在和 之间（10） = + ()/3（11） = + 2()/3（12） = + ()（3）编程绘图编程代码：p=0 0;10 0;n=2; A=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); for k=1:5 d=diff(p)/3;m=4*n-3; q=p(1:n-1,:);p(5:4:m,:)=p(2:n,:); p(2:4:m,:)=q+d; p(3:4:m,:)=q+d+d*A; p(4:4:m,:)=q+2*d; n=m; end plot(p(:,1),p(:,2),k) axis equal axis off（4）运行结果： 实验2.1用高斯消元法的消元过程作矩阵分解。设消元过程可将矩阵A化为上三角矩阵U，试求出消元过程所用的乘数、并以如下格式构造下三角矩阵L和上三角矩阵U验证：矩阵A可以分解为L和U的乘积，即A=LU。编程代码：高斯消元过程：function x=nagauss(a,b,flag)if nargin b=11/6;13/12;47/60b = 1.8333 1.08330.7833(2)b=1.83;1.08;0.783; H=1,1/2,1/3;1/2,1/3,1/4;1/3,1/4,1/5;X=HbX = 1.0800 0.54001.4400结果跟第（1）题的结果相差比较大，则矩阵为变态矩阵实验3.1用泰勒级数的有限项逼近正弦函数用计算机绘出上面四个函数的图形。算法分析：用泰勒级数逼近正弦函数，从泰勒的一阶展开，到二阶展开，再到三阶展开，绘制图形，以查看图形的拟合程度。MATLAB中常用的绘制二维图形的函数是plot函数，其余的函数都是围绕其发展扩充形成的，但是在二阶展开以后fplot函数得拟合效果比plot函数得要好，图形圆滑度更好，更连贯，而且fplot函数代码简单，更为便捷。所以二阶以后使用了fplot函数快速绘图。程序代码：x=0:0.1:pi;y=sin(x);plot(x,y)x=0:0.1:pi; y=x;plot(x,y)x=0:0.1:pi/2; fplot(x-x.3/6,0,pi/2,2e-3)y=x-(x.3)/6;plot(x,y,k) x=0:0.1:pi/2; fplot(x-x.3/6+x.5/120,0,pi/2,2e-3)y=x-(x.3)/6+x.5/120;plot(x,y,k) x=0:0.1:pi; y=sin(x); plot(x,y,k); hold on; x=0:0.1:pi; y=x; plot(x,y,b); hold on; fplot(x-x.3/6,0,pi/2,2e-3,r); hold on; fplot(x-x.3/6+x.5/120,0,pi/2,2e-3,y); hold off;结果分析：图中黑色为正弦曲线，蓝色为一阶泰勒逼近，红色为二阶泰勒逼近，黄色为三阶泰勒逼近，可见黄色逼近效果最好，泰勒级数的阶数越高逼近效果越好。实验3.2 绘制飞机的降落曲线一架飞机飞临北京国际机场上空时，其水平速度为540km/h，飞行高度为1 000m。飞机从距机场指挥塔的横向距离12 000m处开始降落。根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点，则表示飞机距指挥塔的距离，表示飞机的飞行高度，降落曲线为该函数满足条件：（1）试利用满足的条件确定三次多项式中的四个系数；（2）用所求出的三次多项式函数绘制出飞机降落曲线。（1） p=a3,a2,a1,a0; （2） f=polyder(p);（3） p(0)=0;p(12000)=1000;f(0)=0;f(12000)=0;a0 a1 a2a3编程代码：function s=f(x1,y1,x2,y2,x3,y3,x4,y4) format longa1=1,x1,x12,x13 a2=1,x2,x22,x23a3=0,1,2*x3,3*x32a4=0,1,2*x4,3*x42 a=a1;a2;a3;a4 b=y1;y2;y3;y4 s=ab x=-12000:250:0 y=s(3)*x.2-s(4)*x.3 plot(x,y,-r*) xlabel(x/m)ylabel(y/m)运行结果： shiyan3_2(0,0,12000,1000,0,0,12000,0)0 0 0.20833333333333 -0.00001157407407结果分析：shiyan3_2 (0,0,12000,1000,0,0,12000,0)a0 = 1 0 0 0a1 = 1.0e+012 * 0.0000 0.0000 0.0001 1.7280a2 = 0 1 0 0a3 = 0 1 24000 432000000a = 1.0e+012 * 0.0000 0 0 0 0.0000 0.0000 0.0001 1.7280 0 0.0000 0 0 0 0.0000 0.0000 0.0004b = 0 1000 0 0p = 1.0e-004 * 0 -0.0000 0.2083 -0.0000ans = 1.0e-004 * 0 -0.0000 0.2083 -0.0000p = 1.0e-004 * 0 -0.00000000000000 0.20833333333333 -0.00001157407407实验 4.1曾任英特尔公司董事长的摩尔先生早在1965年时，就观察到一件很有趣的现象：集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。因而发表论文，提出了大大有名的摩尔定律（Moores Law），并预测未来这种增长仍会延续下去。下面数据中，第二行数据为晶片上晶体数目在不同年代与1959年时数目比较的倍数。这些数据是推出摩尔定律的依据：年代19591962196319641965增加倍数13456试从表中数据出发，推导线性拟合的函数表达式。算法分析：线性拟合，在数值分析中运用最小二乘法进行多点的线性拟和，带入点后最后求出多项式，此题可用线性函数进行拟合。编程代码：x=1959,1962,1963,1964,1965;y=1,3,4,5,6;p1=polyfit(x,y,1);y1=polyval(p1,x)plot(x,y,rx);hold onplot(x,y1)运行结果：y1 = 0.8113 3.3019 4.1321 4.9623 5.7925实验4.2 参考算法4.2设计绘制Bezier曲线的程序，选取四个点的坐标数据作为控制点绘制飞机机翼剖面图草图的下半部分图形；结合例4.4中上半部分图形绘出完整的机翼草图。最后写出机翼剖面图曲线上20个点处的坐标数据。实验4.3神经元模型用于蠓的分类识别（MCM1989A题）问题描述：生物学字试图对两类蠓虫（Af与Apf）进行鉴别，依据的资料是蠓虫的触角和翅膀的长度，已经测得9只Af和6只Apf的数据（触角长度用x表示，翅膀长度用y表示）Af数据x124136138138138140148154156Y127174164182190170182182208Apf数据x1.141.181.201.261.281.30Y1.781.961.862.002.001.96现需要解决三个问题：（1）如何凭借原始资料（15对数据，被称之为学习样本）制定一种方法，正确区分两类蠓虫；（2）依据确立的方法，对题目提供的三个样本：（1.24,1.80）,(1.28,1.84),(1.40,2.04)加以识别；（3）设Af是宝贵的传粉益虫，Apf是某种疾病的载体，是否应该修改分类方法。问题分析：首先画出15对数据的散点图，其中，Af用*标记，Apf用标记。观察图形，可以发现，Af的点集中在图中右下角，而Apf的点集中在图中左上角。应该存在一条直线L位于两类点之间，作为Af和Apf分界线，这条直线L的确定应依据问题所给的数据，即学习样本。设这条直线的方程为对于平面上任意一点P（x,y），如果该点在直线上，将其坐标代入直线方程则使方程成为恒等式，即使方程左端恒为零；如果点不在直线上，将其坐标代入直线方程，则方程左端不为零。由于f和pf的散点都不在所求的直线上，故将问题所提供的数据代入直线方程左端应该得到表达式的值大于零或者小于零两种不同的结果。这需要建立一个判别系统，引入判别函数,当属于f类时，否则。为了对判别系统引入学习机制，在学习过程中将两种不同的状态，以“”和“”表示。当属于f类时，否则。取于是由所给数据形成约束条件，这是关于判别函数中的三个待定系数的线性方程组：这是包括三个未知数共15个方程的超定方程组，可以求方程组的最小二乘解。（1）根据上面分析写出对应的正规方程组并求解。（2）确定分类边界直线的方程。由所给数据用判别函数判别三个新蠓虫的类属，即当时，判为Af类：当时，判为Apf类。算法分析：此为典型的最小二乘法的求解拟合函数的题型，以达到最佳逼近。记 = 求解最小二乘解时 当可逆时， 此为最小二乘解编程代码：p=1.24,1.36,1.38,1.38,1.38,1.40,1.48,1.54,1.56,1.14,1.18,1.20,1.26,1.28,1.30;q=1.27,1.74,1.64,1.82,1.90,1.70,1.82,1.82,2.08,1.78,1.96,1.86,2.00,2.00,1.96;y=1,1,1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1;p1=p(1:9);q1=q(1:9);p2=p(10:15);q2=q(10:15);A=p,q,ones(size(p);c=(A*A)(A*y) %求方程组的最小二乘解u=1.10:0.02:1.60;v=(-c(1)*u-c(3)/c(2);plot(p1,q1,*,p2,q2,X,u,v) %绘制散点图与判别直线运行结果：c =6.6455 -2.9128 -3.3851实验 5.1 用几种不同的方法求积分的值。（1）牛顿-莱布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）复合梯形公式。算法分析：(2)梯形公式 S=()(3)Simpson 公式 （4）复合梯形公式编程代码：syms xI1=int(4/(1+x2),x,0,1)a=0;b=1;h=b-a;I2=(4/(1+a2)+4/(1+b2)/2I3=h/6*(4/(1+a2)+4*4/(1+(a+b)/2)2)+4/(1+b2)a=0;b=1;M=10;h=(b-a)/M;s=0;for k=1:(M-1) x=a+h*k; s=s+4/(1+x2);ends=h*(4/(1+a2)+4/(1+b2)/2+h*s;I4=s结果输出：I1 =piI2 =3I3 = 3.1333I4 = 3.1399结果分析：由牛顿-莱布尼茨公式做出的答案是该积分式的精确解，将I2，I3，I4与精确解对比可以发现，梯形公式的计算精度最低，复化Simpson公式的计算结果精度最高。但是复化梯形公式的表达是最为复杂，在相对比较精度要求时，可选择相对精度相对低的梯形公式和辛普生公式，比如在求初值的欧拉公式选择了比梯形公式精度还低的矩形公式，在改进的欧拉公式中选择了梯形公式时精度便大打得到了提高，达到了精度要求。实验5.2 设X为标准正态随机变量，即XN（0，1）。现分别取，试设计算法计算30个不同的概率值；，并将计算结果与概率论教科书中的标准正态分布函数表作比较。（提示： 编程代码：function f=f (x)f=exp(-0.5*x.2);在Command Window中键入如下指令：u=0;for k=1:31 p(k)=quad(shiyan5_2,u,4); u=u+0.1;endp结果输出：p = 1.2532 1.1534 1.0546 0.9577 0.8637 0.7733 0.6874 0.6064 0.5310 0.4613 0.3976 0.3400 0.2884 0.2426 0.2023 0.1674 0.1373 0.1116 0.0900 0.0719 0.0569 0.0447 0.0348 0.0268 0.0205 0.0155 0.0116 0.0086 0.0063 0.00460.0033结果分析：将此结果与概率论教科书中的标准正态分布表做比较，发现基本一致。实验53 设某城市男子的身高XN（170，36）（单位：cm），应如何选择公共汽车门的高度H使男子与车门碰头的机会小于1%。问题分析：由题设男子身高XN（170，36）（单位：cm），应如何选择公共汽车门的高度H使男子与车门碰头的机会小于1%。问题分析：由题设男子身高数据服从平均值为170（cm），方差为6（cm）的正态分布，其分布密度函数为按正态分布的分布规律（原则），这个城市的男子身高超过188（cm）的人数极少。故可以对H=188，187，186，求出概率的值，观察使概率不超过1%的H，以确定公共汽车门应该取的高度。概念值的计算实际上是求定积分（1）选用一种数值求积公式或用数学软件分别计算出H=180、181、188时定积分近似值。（2）根据上面计算的积分值，按题目要求确定公共汽车门的高度取值（答案184cm）。如果将汽车门的高度取180cm，是否满足大多数市民的利益？（3）用计算机模拟的方法来检验你的结论，计算机产生10 000个正态随机数（它们服从均值为170，方差为6的正态分布）来模拟这个城市中10 000个男子的身高，然后统计出这10 000人中身高超过180（cm）的男子数量所占的百分比。编程代码：function y=f(x)y=exp(-(x-170).2/72)/(6*sqrt(2)*pi);在Command Window中键入如下指令u=180;for k=1:9 p(k)=quad(shiyan5_3(x),u,194); u=u+1;endpp = 0.0269 0.0188 0.0128 0.0085 0.0055 0.0035 0.0021 0.0013 0.0007 0.4613 0.3976 0.3400 0.2884 0.2426 0.2023 0.1674 0.1373 0.1116 0.0900 0.0719 0.0569 0.0447 0.0348 0.0268 0.0205 0.0155 0.0116 0.0086 0.0063 0.0046 0.0033结果分析：可以看出来，要是车门符合条件，满足大多数市民的要求，车门要设计为184cm。实验验证：F0.99，即 F故在指令窗口键入h=icdf(norm,0.99,170,6)求得 h =1.839580872442451e+002也就是184cm，所以车门高度定位180cm，不能满足大多数市民的利益实验54 用蒙特卡罗方法计算二重积分，其中，算法分析： 二重积分的几何意义是以为曲面顶，S为底的柱体D的体积。用一下思想求I的近似值：假设S被包括在几何体D的内部，D的体积已知，若在D内部产生1个均匀的随机数，那么 P（随机数落在S内） C的体积/D的体积现产生在D内的N个均匀分布的随机数。若其中个落在S的内部，那么 I D的体积/N编程代码：clear;N=1000;x=unifrnd(0,1,N,1);y=unifrnd(0,1,N,1);z=rand(N,1);c1=(x-0.5).2+(y-0.5).2)1;c2=(z-sin(sqrt(log(x+y+1) shiyan6_2(4,5)ans =85.76257352185962实验63 将下列高阶常微分方程初值问题化为一阶常微分方程组的初值问题，然后用MATEAB指令ode23求数值解。算法分析：将高阶微分方程化为低阶微分方程编程代码：function f=f(t,y)f=y(2);y(3);y(4);80-2*y(4)-y(2)-y(1);在指令窗口键入 y0=0,0,0,0; t,y=ode23(shiyan6_3,0,10,y0)结果输出：y = 1.0e+003 * 0 0 0 0 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.0000 0.0003 0.0000 0.0000 0.0000 0.0007 0.0000 0.0000 0.0000 0.0013 0.0000 0.0000 0.0000 0.0018 0.0000 0.0000 0.0000 0.0023 0.0000 0.0000 0.0001 0.0028 0.0000 0.0000 0.0001 0.0033 0.0000 0.0000 0.0001 0.0039 0.0000 0.0000 0.0001 0.0044 0.0000 0.0000 0.0002 0.0051 0.0000 0.0000 0.0002 0.0059 0.0000 0.0000 0.0003 0.0068 0.0000 0.0000 0.0004 0.0078 0.0000 0.0000 0.0006 0.0087 0.0000 0.0000 0.0007 0.0096 0.0000 0.0000 0.0008 0.0105 0.0000 0.0001 0.0010 0.0114 0.0000 0.0001 0.0012 0.0124 0.0000 0.0001 0.0015 0.0134 0.0000 0.0001 0.0018 0.0145 0.0000 0.0002 0.0021 0.0157 0.0000 0.0002 0.0025 0.0169 0.0000 0.0003 0.0031 0.0182 0.0000 0.0004 0.0037 0.0196 0.0001 0.0006 0.0044 0.0210 0.0001 0.0008 0.0053 0.0224