数值分析课程设计讲解.docx

《数值分析课程设计讲解.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计讲解.docx（31页珍藏版）》请在三一文库上搜索。

1、实验1.1水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个 水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又 去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰 子分成五堆，恰多一只给猴子，私藏一堆，再去入睡，天亮以后，大家把余下的 椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰 子？算法分析：求解这一问题可以用迭代递推算法。首先分析椰子数目的变化规律， 设最初的椰子数为P 0,即第一个水手所处理之前的椰子数， 用P 1、P 2、p 3、

2、P4、 P5分别表示五个水手对椰子动了手脚以后剩余的椰子数目，则根据问题有4Pk 1 = 5(Pk一1)，(k=0,1,2,3,4)再用x表示最后每个水手平分得到的椰子数，于是有x = 1(P5-1)5所以P5 = 5x 1利用逆向递推的方法，有5PkPk1 1， (k=4,3,2,1,0)4有了逆向递推关系式，求解这一问题似乎很简单，但由于椰子数为一正整数， 用任意的x作为初值递推出的P0数据不一定是合适的。这里用for 循环语句结合break 语句来寻找合适的x和p0 ,对任意的x 递推计算出Po ,当计算结果为正整数时，结果正确，否则选取另外的 x再次 重新递推计算，直到计算出的结果 P

3、o为正整数为止。MATLA射程代码：(1) n=inPut(输出 n 的值：);for x=1:np=5*x+1;for k=1:5 p=5*p/4+1;endif p=fix(p)breakendenddisp(x,p)输出结果：输出 n 的值： 1500102315621(2) for x=1:infp=5*x+1;for k=1:5p=5*p/4+1;endif p=fix(p) breakendenddisp(x,p)输出结果：102315621C语言编程代码：#include int count(int);void main()int i, n, y;printf( 输入水手数： n

4、 );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 4)ok+; k=+m; i=1;if(i=a&ok=4) break;return(k*5+1);实验1 . 3绘制Koch分形曲线问题描述：从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5 个结

5、点的新的图形（图 1） ；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替， 再次形成新的图形（图 2） ，这时，图形中共有17 个结点。这种迭代继续进行下去可以形成 Koch 分形曲线。在迭代过程中，图形中的结点将越来越多，而曲线最终显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。 Koch 分形曲线的绘制与算法设计和计算机实现相关。问题分析：考虑由直线段（ 2 个点） 产生第一个图形（ 5 个点）的过程， 设 P1和 P5 分别为原始直线段的两个端点。现在需要在直线段的中间依次插入三个点P2,P3,P4产生第一次迭代的图形（图1）。显然，P2位于P

6、i点右端直线段的三分之一处，R点绕P2旋转60度(逆时针方向)而得到的，故可以处理为向量 P2 P4经正交变换而得到向量P2P3,形成算法如下：(1)P2 =P +(F5 P)/3;(2) R+2(F5P)/3;(3)2=P2 +(P4-P2)X AT ；在算法的第三步中,A为正交矩阵。冗icos3冗 ,sin -3ji sin3in cos3这一算法将根据初始数据(Pi和P5点的坐标)，产生图1中5个结点的坐标。 这5个结点的坐标数组，组成一个5X2矩阵。这一矩阵的第一行为为Pi的坐标， 第二行为P的坐标，第二行为PJ勺坐标第五行为P5的坐标。矩阵的第一列 元素分别为5个结点的x坐标，第二列

7、元素分别为5个结点的y坐标。问题思考与实验：(1)考虑在Koch分形曲线的形成过程中结点数目的变化规律。设第 k次迭 代产生结点数为nk,第k+1迭代产生结点数为nk书，试写出nk和nk书之间的递推 关系式；(2)参考问题分析中的算法，考虑图1到图2的过程，即由第一次迭代的5 个结点的结点坐标数组，产生第二次迭代的17个结点的结点坐标数组的算法；(3)考虑由第k次迭代的nk个结点的结点坐标数组，产生第 k+1次迭代的 nk ,个结点的结点坐标数组的算法；(4)设计算法用计算机绘制出如下的 Koch分形曲线(图3)。算法分析：(1)第k次迭代产生的结点数 为 错误！未找到引用源。，第k+1次迭

8、代产生的结点数为错误！未找到引用源。错误！未找 到引用源。（2） 第一次迭代的 5 个结点的结点坐标数组， 错误！未找到引用源。 错误！未找到引用源。 错误！未找到引用源。 错误！未找到引用源。在 错误！未找到引用源。 和 错误！未找到引用源。 之间（ 1） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + （ 错误！ 未找到引用源。 ）/3（ 2） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + 2（错误！未找到引用源。 ）/3（ 3） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + （错误！未找到引用源。 ）X错误！未找到引用源。在 错误！未找到引用源。 和 错误！未

9、找到引用源。 之间（ 4） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + （ 错误！ 未找到引用源。 ）/3（ 5） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + 2（错误！未找到引用源。 ）/3（ 6） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + （错误！未找到引用源。 ）X错误！未找到引用源。在 错误！未找到引用源。 和 错误！未找到引用源。 之间（ 7） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + （错误！ 未找到引用源。 ）/3（ 8） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + 2（错误！未找到引用源。 ）/3（ 9） 错误！

10、未找到引用源。 = 错误！ 未找到引用源。 + （错误！未找到引用源。 ）X错误！未找到引用源。在 错误！未找到引用源。 和 错误！未找到引用源。 之间（ 10） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + （错误！未找到引用源。 ）/3（ 11） 错误！未找到引用源。 = 错误！未找到引用源。 + 2（错误！未找到引用源。 ）/3（ 12） 错误！ 未找到引用源。 = 错误！ 未找到引用源。 + （错误！未找到引用源。 ）X错误！未找到引用源。（ 13） 绘图编程代码：p=0 0;10 0;n=2;A=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3

11、); for k=1:5d=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;endplot(p(:,1),p(:,2),k)axis equalaxis off4）运行结果：用高斯消元法的消元过程作矩阵分解。设20 23A= 1812-3 15_消元过程可将矩阵A化为上三角矩阵U,试求出消元过程所用的乘数m2i、m3i、m3i并以如下格式构造下三角矩阵L和上三角矩阵U120 23L = m21 1 ,U =a22) a23(2)口3

12、1 m)32 1 Ja33验证：矩阵A可以分解为L和U的乘积，即A=LU。编程代码：高斯消元过程：function x=nagauss(a,b,flag)if nargin b=11/6;13/12;47/60 b =1.83331.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.08000.54001.4400结果跟第（1）题的结果相差比较大，则矩阵为变态矩阵实验3.1用泰勒级数的有限项逼近正弦函数y0(x) =sinx, x 0,二yi(x) = x, x 0，二 /23y2(x) =

13、x-x /6,x 0,二 /235_y3(x) =x-x /6 x /120, x 0,二 /2用计算机绘出上面四个函数的图形 算法分析：用泰勒级数逼近正弦函数，从泰勒的一阶展开，到二阶展开，再到三阶展开， 绘制图形，以查看图形的拟合程度。MATLAB中常用的绘制二维图形的函数是 plot函数，其余的函数都是围绕 其发展扩充形成的，但是在二阶展开以后fplot函数得拟合效果比plot函数得要好，图形圆滑度更好，更连贯，而且 fplot函数代码简单，更为便捷。所以二阶 以后使用了 fplot函数快速绘图。程序代码: x=0:0.1:pi; y=sin(x);plot(x,y)x=0:0.1:pi

14、 y=x； piot(x,y)fplot(x-xA3/6,0,pi/2,2e-3)x=0:0.1:pi/2;y=x-(x.A3)/6; plot(x,y,k) |x=0:0.1:pi/2;fplot(x-xA3/6+xA5/120,0,pi/2,2e-3)y=x-(xA3)/6+xA5/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-xA3/6,0,pi/2,2e-3,r);hold on;fplotCx-x.A3/6+x.A5/120,0,p

15、i/2,2e-3,y);hold off;结果分析：图中黑色为正弦曲线，蓝色为一阶泰勒逼近，红色为二阶泰勒逼近， 黄色为三阶泰勒逼近，可见黄色逼近效果最好，泰勒级数的阶数越高逼近效果越 好。实验3.2绘制飞机的降落曲线一架飞机飞临北京国际机场上空时，其水平速度为540km/h,飞行高度为1000m。飞机从距机场指挥塔的横向距离 12 000m处开始降落。根据经验，一架 水平飞行的飞机其降落曲线是一条三次曲线。建立直角坐标系，设飞机着陆点为 原点O,降落的飞机为动点P(x,y),则x表示飞机距指挥塔的距离，y表示飞机 的飞行高度，降落曲线为y(x); a0 a1x a2x2 a3x3该函数满足条

16、件：y(0) =0, y(12000) =1 000y(0) =0, y(12000) =0(1)试利用y(x)满足的条件确定三次多项式中的四个系数；(2)用所求出的三次多项式函数绘制出飞机降落曲线。(1) p=a3,a2,a1,a0;(2) f=polyder(p);(3) p(0)=0；p(12000)=1000;f(0)=0;f(12000)=0;a0a1a2a3编程代码：function S=f(x1,y1,x2,y2,x3,y3,x4,y4)format longa1=1,x1,x1A2,x1A3a2=1,x2,x2A2,x2A3a3=0,1,2*x3,3*x3A2a4=0,1,2*

17、x4,3*x4A2a=a1;a2;a3;a4b=y1;y2;y3;y4s=abx=-12000:250:0y=s(3)*x.A2-s(4)*x.A3plot(x,y,-r*)xlabel(x/m)ylabel(y/m)运行结果： shiyan3_2(0,0,12000,1000,0,0,12000,0)000.20833333333333-0.000011574074071000900800700600300200100m 500400a00, a1 0,a2 0.20833333333333,0-12000-10000-8000-6000-4000-2000x/m结果分析:330.00001

18、15740740700.00011.72800shiyan3_2 (0,0,12000,1000,0,0,12000,0)a0 = 100a1 = 1.0e+012 *0.00000.0000a2 =010a3 =240004320000001.0e+012 *0.00000000.00000.00000.00011.728000.00000000.00000.00000.0004b =01000 0 0P =1.0e-004 * 0-0.00000.2083-0.0000ans =1.0e-004 * 0-0.00000.2083-0.0000P =1.0e-004 *0-0.0000000

19、00000000.20833333333333-0.00001157407407实验4.1曾任英特尔公司董事长的摩尔先生早在 1965年时，就观察到一件很有趣的现象： 集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。 因而发表论文，提出了大大有名的摩尔定律（Moore s LaW,并预测未来这种增长仍会延续下去。下面数据中，第二行数据为晶片上晶体数目在不同年代与 1959 年时数目比较的倍数。这些数据是推出摩尔定律的依据：年代19591962196319641965增加倍数13456试从表中数据出发，推导线性拟合的函数表达式算法分析：线性拟合，在数值分析中运用最小二乘法

20、进行多点的线性拟和，带入点后最后求出多项式，此题可用线性函数进行拟合。编程代码：x=1959,1962,1963,1964,1965;y=1,3,4,5,6;p1=poly攸x,y,1);y1=polyval(p1,x)plot(x,y,rx);hold onplot(x,y1)运行结果：y1 =0.81133.30194.13214.96235.79251959196019611962196319641965验4.2参考算法4.2设计绘制Bezier曲线的程序，选取四个点的坐标数据作为控制点绘 制飞机机翼剖面图草图的下半部分图形； 结合例4.4中上半部分图形绘出完整的 机翼草图。最后写出机翼

21、剖面图曲线上 20个点处的坐标数据。实验5.1用几种不同的方法求积分 Hdx的值01 x2(1)牛顿-莱布尼茨公式；(2)梯形公式；(3)辛卜生公式；(4)复合梯形公式。算法分析：(2)梯形公式S=错误！未找到引用源。(错误！未找到引用源。)(3)Simpson公式 错误！未找到引用源。h门(4)复合梯形公式 Tn(f (a) f (b) h f(Xk)2k编程代码：syms xI1=int(4/(1+xA2),x,0,1)a=0;b=1;h=b-a;I2=(4/(1+aA2)+4/(1+bA2)/2I3=h/6*(4/(1+aA2)+4*4/(1+(a+b)/2)A2)+4/(1+bA2)a

22、0;b=1;M=10;h=(b-a)/M;s=0;for k=1:(M-1)x=a+h*k;s=s+4/(1+xA2);ends=h*(4/(1+aA2)+4/(1+bA2)/2+h*s;I4=s结果输出：11 =pi12 =313 =3.133314 =3.1399结果分析：由牛顿-莱布尼茨公式做出的答案是该积分式的精确解，将I2 , I3 , I4与精确解对比可以发现，梯形公式的计算精度最低，复化Simpson公式的计算结果精度最 高。但是复化梯形公式的表达是最为复杂，在相对比较精度要求时，可选择相对精度相对低的梯形公式和辛普生公式，比如在求初值的欧拉公式选择了比梯形公 式精度还低的矩形

23、公式，在改进的欧拉公式中选择了梯形公式时精度便大打得到 了提高，达到了精度要求。实验5.2设X为标准正态随机变量，即XN (0, 1)。现分别取u = 0.1,0.2,0.3,|,2.9,3.0 ,试设计算法计算30个不同的概率值；PX Ua,并将计算结果与概率论教科 书中的标准正态分布函数表作比较。22,g一.1 二:14:(提小：P IX uaJ= e 2 dx = e 2dx).2 二,2 二 Ua编程代码：function f=f (x)f=exp(-0.5*x.A2);在Command Window中键入如下指令:u=0;for k=1:31p(k)=quad(shiyan5_2,u

24、4);u=u+0.1;endP结果输出：p =1.25321.15341.05460.95770.86370.77330.68740.60640.53100.46130.39760.34000.28840.24260.20230.16740.13730.11160.09000.07190.05690.04470.03480.02680.02050.01550.01160.00860.00630.00460.0033结果分析：将此结果与概率论教科书中的标准正态分布表做比较，发现基本一致实验5. 4用蒙特卡罗方法计算二重积分 fsin Jln(x +y +1)dxdt,其中，D尸9911 )2(

25、1 )21d =i(x,y)| x-+ y-I2)I2J4算法分析：二重积分的几何意义是以 错误！未找到引用源。为曲面顶，S为底的柱体 D的体积。用一下思想求I的近似值：假设S被包括在几何体D的内部，D的体 积已知，若在D内部产生1个均匀的随机数，那么P (随机数落在S内)错误！未找到引用源。C的体积/D的体积现产生在D内的N个均匀分布的随机数。若其中 错误！未找到引用源。个落在 S的内部，那么 I错误！未找到引用源。D的体积X Ns/N编程代码：clear;N=1000;x=unifrnd(0,1,N,1);y=unifrnd(0,1,N,1);z=rand(N,1);c1=(x-0.5)A

26、2+(y-0.5)A2)1;c2=(z-sin(sqrt(log(x+y+1) y0=0,0,0,0; t,y=ode23(shiyan6_3,0,10,y0)结果输出：y =1.0e+003 *00.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00010.00000.00000.00000.00030.00000.00000.00000.00070.00000.00000.00000.00130.00000.00000.0

27、0000.00180.00000.00000.00000.00230.00000.00000.00010.00280.00000.00000.00010.00330.00000.00000.00010.00390.00000.00000.00010.00440.00000.00000.00020.00510.00000.00000.00020.00590.00000.00000.00030.00680.00000.00000.00040.00780.00000.00000.00060.00870.00000.00000.00070.00960.00000.00000.00080.01050.0

28、0000.00010.00100.01140.00000.00010.00120.01240.00000.00010.00150.01340.00000.00010.00180.01450.00000.00020.00210.01570.00000.00020.00250.01690.00000.00030.00310.01820.00000.00040.00370.01960.00010.00060.00440.02100.00010.00080.00530.02240.00010.00100.00630.02390.00020.00140.00760.02540.00030.00190.0

29、0900.02680.00040.00250.01080.02830.00060.00330.01280.02960.00090.00440.01530.03080.00140.00590.01810.03180.00210.00800.02140.03250.00320.01070.02520.03280.00480.01430.02950.03260.00740.01920.03450.03150.01140.02590.03990.02920.01780.03500.04570.02510.02840.04750.05130.01810.04700.06520.05530.00560.0

30、6590.07970.0552-0.00680.08850.09370.0515-0.02140.11760.10740.0428-0.03950.15270.11840.0273-0.06060.19170.12370.0043-0.08310.23190.1202-0.0266-0.10480.26430.1084-0.0584-0.12090.29400.0848-0.0973-0.13360.31530.0470-0.1409-0.13940.3223-0.0062-0.1858-0.13480.3093-0.0748-0.2278-0.11640.2774-0.1457-0.2578

31、0.08730.2196-0.2290-0.2781-0.04160.1318-0.3192-0.28170.02190.0289-0.3981-0.26560.0921-0.0983-0.4707-0.22690.1748-0.2187-0.5198-0.17550.2500-0.3697-0.5579-0.09310.3404-0.5393-0.56910.02340.4370-0.7133-0.53900.17440.5295-0.8604-0.46480.33910.5998-0.9850-0.33010.53320.6479-1.0599-0.12560.74350.6571-1.

32、06570.08760.90970.6262-1.01110.34221.06290.5550-0.88800.62991.19060.4402最后得到的数值解为0.4402实验6. 4列出函数f (x) =1 -lnu(x)u(x)在区间0, e上的函数值表并作出它的图象。其中，u(x)是初值问题dln (x) =2xdx(0) = 0的解。编程代码：v=dsolve(Dv*log(x)=2*x,v(0)=0,x)f=(1-log(v)*vezplot(f)结果输出：v =-2*Ei(1,-2*log(x)f =-2*(1-log(-2*Ei(1,-2*log(x)*Ei(1,-2*log(x)f6420-2-4-6-6-4-20246f实验6. 5鱼雷击舰问题。一敌舰在某海域内沿正北方向航行时,的正西方1 n mile处。我舰向敌舰发射鱼雷，敌舰速度为 速度为敌股速度的两部。试问敌舰航行多远时将被击中？我方战舰恰位于敌舰0.42n mile/min ,鱼雷提示：本问题归结为常微分方程初值问题. 1 y2y=-JL-,x (0,1)2(1-x)y |x=0 = 0, y |x=