数值分析课程课程设计Word格式文档下载.docx

1、解该问题主要使用递推算法，关于椰子数目的变化规律 可以设起初的椰子数为po,第一至五次猴子在夜里藏椰子后，椰子的数 目分别为po,p1,p2,P3,P4再设最后每个人分得x个椰子，由题：r 1 7（ pk 1） （k=0,1,2,3,4 ） x -（ p5 1）5 5所以P5 5x 1, p p if i=0.1&i0&=0.0001for n=30:-1:2i=（-1/5）*i+1/（5*n）i =1.1336e+005-5.6679e+0052.8339e+006-1.4170e+0077.0848e+007-3.5424e+0081.7712e+009-8.8560e+0094.4280

2、e+010-2.2140e+011同理输入积分初始值i=0时可以得i=0.0884第二种方法所得的结果相对来说比较精确一些，也比较可靠 因为第一种方法每一迭代都将最初的误差放大了五倍， 使得最终的误差越来越大；而第一种方法经过每一次迭代都将误差缩小为初始误差的五 分之一，使得最终的误差越来越小，因此相对来说比较可靠，性能较好。1.3绘制Koch分形曲线问题描述：从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有 5个结点C0SJ的国写（图1-4）;在新 的图形中，又将图中每一直线段中间的询之丁部分都用一个等边三角 形的另两条边代替，再次形成新的图形（图sit5）

3、,cos时，图形中共有17个结点。这种迭代继续进行下去可以形成 Koch分形曲线。问题分析：考虑由直线段（2个点）产生第一个图形（5个点）的过程， 设p1和R分别为原始直线段的两个端点。现在需要在直线段的中间依次 插入三个点P2，R，P4产生第一次迭代的图形（图1-4）。显然，只位于P点 右端直线段的三分之一处，已点绕p2旋转60度（逆时针方向）而得到 的，故可以处理为向量 巳已经正交变换而得到向量P2P3,形成算法如下：（1）P2P （任 P）/3.（2）RP 2（P5 R）/3.（3）P3% （E E） AT.在算法的第三步中，A为正交矩阵。这一算法将根据初始数据（P和R点的坐标），产生图

4、1-4中5个结点的坐标。这5个结点的坐标数组，组成一个 5X 2矩阵。这一矩阵的第一 行为为R的坐标，第二行为R的坐标，第二行为P2的坐标第五行为P5 的坐标。矩阵的第一列元素分别为 5个结点的x坐标，第二列元素分 别为5个结点的y坐标。问题思考与实验：（1）考虑在Koch分形曲线的形成过程中结点数目的变化规律。设第 k 次迭代产生结点数为nk,第k 1迭代产生结点数为nk 1,试写出nk和nk 1之 间的递推关系式；（2）参考问题分析中的算法，考虑图1-4到图1-5的过程，即由第一次 迭代的5个结点的结点坐标数组，产生第二次迭代的 17个结点的结点 坐标数组的算法；（3）考虑由第k次迭代的n

5、k个结点的结点坐标数组，产生第k 1次迭代 的nk 1个结点的结点坐标数组的算法；（4）设计算法用计算机绘制出如下的（1） nk1 4nk 3Koch分形曲线（图1-6）（2）（3辩法及（4）代码分析：p=0 0;10 0; %P 为初始两个点的坐标，第一列为x坐标，第二列为y坐标n=2; %n 为结点数A=cos（pi/3） -sin（pi/3）;sin（pi/3） cos（pi/3）; % 旋转矩阵for k=1:d=diff（p）/3; %diff %m=4*n-3; %q=p（1:n-1,:p（5:4:m,:）=p（2:n,: % p（2:）=q+d;p（3:）=q+d+d*A; %

6、p（4:）=q+2*d; % n=m;plot（p（:,1）,p（:,2） %axis（0 10 0 3.5）计算相邻两个点的坐标之差，得到相邻两点确定的向量则d就计算出每个向量长度的三分之一 ，与题中将线段三等分对应迭代公式以原点为起点，前n-1个点的坐标为终点形成向量迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标用向量方法计算迭代后处于4k+2位置上的点的坐标4k+3位置上的点的坐标4k位置上的点的坐标迭代后新的结点数目2.1用高斯消元法的消元过程作矩阵分解。设20 2 3A 1 8 12 3 15消元过程可将矩阵 A化为上三角矩阵U,试求出消元过程所用的乘数m21、m31、伉1并以

7、如下格式构造下三角矩阵L和上三角矩阵U1 20 2 3L m21 1 ,U a22） a23）4 （2）m31 m32 1 a33验证：矩阵A可以分解为L和U的乘积，即A=LU。矩阵LU分解MATLAB代码：function hl=zhjLU（A）n,n=size（A）;RA=rank（A）;if RA=ndisp（因为A的n阶行列式hl等于零，所以A不能进行LU分解.A的秩RA如下：）；RA,hl=det（A）;returnif RA=nfor p=1：nh（p）=det（A（1：p,1：p）;hl=h（1：n）;for i=1：if h（1,i）=0 因为A的各阶主子式等不等于零，所以A能

8、进彳f LU分解.A的秩RA和各阶顺序主子式值hl庆次如下：RA,hlif h（1,i）=0 因为A的各阶主子式都 不等于零，所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl如 下：for j=1：U（1,j）=A（1,j）;for k=2:for i=2:for j=2:L（1,1）=1;L（i,i）=1;jL（2,1）=A（2,1）/U（1,1）;L（i,1）=A（i,1）/U（1,1）;L（i,k）=（A（i,k）-L（i,1:k-1）*U（1:k-1,k）/U（k,k）;elseU（k,j）=A（k,j）-L（k,1:k-1,j）;RA,hl,U,L仅上上代码保存为M文件，并在命

9、令窗口输入A=20 2 3;1 8 1; 2 -3 15;b=0 0 0h=zhjLU（A）程序运行结果:L =U=1.00000 0 20.0000 2.0000 3.00000.05001.0000 0 0 7.9000 0.85000.1000-0.4051 1.0000 0 0 15.0443实验验证：可以直接使用MATLAB内置LU分解L U=lu（A） 输出结果与上程序输出结果一致。2. 2用矩阵分解方法求上题中 A的逆矩阵。记1b 0 ,b21 ,b31分别求解方程组 AX b1，AX b2，AX由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。对第一个方程组，由于A=LU

10、，所以先求解下三角方程组LY知再求解上三角 方程组UX Y,则可得逆矩阵的第一列列向量；类似可解第二、第三方 程组，得逆矩阵的第二列列向量的第三列列向量。由三个列向量拼装可 得逆矩阵A1。解:MATLAB代码如下:b1=1;0;0; b2=0;1; b3=0;1;A=20,2,3;1,8,1;2,-3,15;L=1,0,0;0.05,1,0;0.1,-0.4051,1;U=20 2 3;0 7.9 0.85;0 0 15.0443;Y1=Lb1X1=UY1Y2=Lb2 X2=UY2 Y3=Lb3X3=UY3Y1 =-0.0500-0.1203X1 =0.0517-0.0055-0.0080Y2

11、 =X1 X2 X3 ans =0.0517 -0.0164 -0.0257-0.0055 0.1237 0.1165-0.0080 0.0269 0.0934而：inv （A） =X1 X2 X3 得证0.4051X2 =-0.01640.12370.0269Y3 =1.4051X3 = -0.0257 0.1165 0.09342.3验证希尔伯特矩阵的病态性：对于三阶矩阵1 1/2 1/3H 1/2 1/3 1/41/3 1/4 1/5取右端向量b【11/6 13/12 47/6t,验证：（1） 向量X X1 X2 X3T 1 1 1T是方程组HX b的准确解；（2）取右端向量b的三位有效

12、数字得b 1.83 L08 0.783T,求方程组 的准确解X ,并与X的数据I1,1,1】作比较。说明矩阵的病态性。H=1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5;X=1；1；1；b=H*Xb =1.83331.08330.7833与题中相同（2）先求出解X,与数据，亍作比较b=1.83;1.08;0.783;X=HbX =1.08000.54001.4400与1,1,1相差较大，矩阵为病态矩阵3.1用泰勒级数的有限项逼近正弦函数sin x,x 0,x,x 0, /23 x x /6,x 0, /2 x x3/6 x5/120, x用计算机绘出上面四个函数的图形。MA

13、TLAB代码如下 （1） syms x; taylor（sin（x） x=0:0.01*pi:pi plot（x,sin（x）syms x;taylor（x）x=0:pi/2plot（x）taylor（x-xA3/6）fplot（x-xA3/6,0 pi/2） （4）taylor（x-xA3/6+xA5/120）x-xA3/6+xA5/120,0 pi/2）结果图形右：0.1:pi;y=sin（x）;plot（x,y,-skhold onpi/2;y=x;-b*）x-x.A3/6,0,pi/2,0,2,2e-3,-gxhold on fplot（x-x.A3/6+x.A5/120-r.） ho

14、ld offlegend（sin（x）,xx-A3/6+xA5/120,2）xlabel（ylabel（ytitle（Taylor approximation图中红色点线为正弦曲线，蓝色的星线为一阶泰勒逼近，绿色叉线为二阶豪勒逼近，黑色正方形线为三阶泰勒逼近。可见三阶泰勒逼近效果最 好，泰勒级改越高，逼近效果尬婷。3.2绘制飞机的降落曲线一架飞机飞临北京国际机场上空时，其水平速度为 540km/h,飞行高 度为1 000m。飞机从距机场指挥塔的横向距离12 000m处开始降落。根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。 建立直 角坐标系，设飞机着陆点为原点O，降落的飞机为动点P（x

15、,y），则x表 示飞机距指挥塔的距离，y表示飞机的飞行高度，降落曲线为2 3y（x） a0 a1x a2x a3x该函数满足条件：y（0） 0, y（12 000） 1000（0） 0, y（12000） 0（1）试利用y（x）满足的条件确定三次多项式中的四个系数；（2）用所求出的三次多项式函数绘制出飞机降落曲线。function S=f（x1,y1,x2,y2,x3,y3,x4,y4） format longa1=1 ,x1,x1A2,x1A3;a2=1 ,x2,x2A2,x2A3;a3=0,1,2*x3,3*x3A2;a4=0,1,2*x4,3*x4A2;a=a1; a2;a3;a4;b=

16、y1; y2;y3;y4;s=ab;x=-12000:250:y=s（3）*x.A2-s（4）*x.A3-d以上为M文件内容，在命令窗口键入f（0,0,12000,1000,0,0,12000,0）运行结桌： ans为多项土索数1.0e-004 *0.20833333333333-0.000011574074074.1曾任英特尔公司董事长的摩尔先生早在 1965年时，就观察到 一件很有趣的现象：集成电路上可容纳的零件数量，每隔一年半左右 就会增长一倍，性能也提升一倍。因而发表论文，提出了大大有名的年代19591962196319641965增加倍数356摩尔定律（Moore s LaW,并预测

17、未来这种增长仍会延续下去。下面 数据中，第二行数据为晶片上晶体数目在不同年代与 1959年时数目比较的倍数。这些数据是推出摩尔定律的依据：试从表中数据出发，推导线性拟合的函数表达式。解法一MATLAB：码：x=1959,1962,1963,1964,1965;y=1,3,4,5,6;p1=polyfit（x,y,1）y1=polyval（p1,x）plot（x,y1,-,x,y,r*）,ylabel（运行结果：p1 =1.0e+003 *0.0008 -1.6255y1 =0.8113 3.3019 4.1321 4.9623 5.7925线性拟合的函数表达式：Y=0.8302x-1.6255

18、e+003解法二：x=1959 1962 1963 1964 1965;xs =-1.6255283018912380.000830188679248y=1;3;4;5;6;for i=1:length（x）for j=1:A（i,j）=x（i）”1）;L,U=lu（A*A）;xs=U（L（A*y）从而年代y与增加倍数x之间的关系为：y=-1625.528301891238+0.830188679248x11=pi12= 313=3.133314=3.1416（1）牛顿-莱布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）复 合梯形公式。syms xi1=int（4/（1+xA2）,x,0,

19、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）M=100;h=（b-a）/M;i4=0;（M-1）x=a+h*k;i4=i4+4/（1+x2）;i4=h*（4/（1+aA2）+4/（1+bA2）/2+h*i4牛顿莱布尼兹公式得到精确结果a采用挣形公式得到的结果比采用 Simpson公式的精蒲I度鱼很，米用复祝榜形公式在步长取得越束裁小 的状态下可以提高精度。5. 5求空间曲线L:y sintz 2 cost sint弧长公式为L : .x2（t） y2（t） z2（t）

20、dtx cost6. 1用欧拉公式和四阶龙格-库塔法分别求解下列初值问题;0.9y（i）y , y（0） i;x 0,11 2x2；y 闩顶。（1） 欧拉公式：function t,x=Euler（fun,t0,tt,x0,N） h=（tt-t0）/N;t=t0+0:N*h;x（1,:）=x0Nf=feval（fun,t（k）,x（k,:）;f=f;x（k+1,:）=x（k,:）+h*f;以Euler.m保存function f=Euler_fun（t,x） f=0.9*x./（1+2*t） 以Euler_fun.m保存 function main_Euler t,x=Euler（Euier_fun,0,1,1,20）;fh=dsolve（Dx=0.9*x/（1+2*t）x（0）=18ft（k）=t（k）;fx（k）=subs（fh,ft（k）;t,x , 4main_Euler.m保存 输入：main_Euler四阶龙杯-库塔法function R=rk4（f,a,b,ya,N）%y=f（x,y）%a,b为左右端点%Nf