数值分析课程设计.docx

数值 课分 程析 设计专业班级：信息与计算科学09-1班姓 名：陈育伟 学 号：20096361 实验一11 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题（15621）。【问题分析】首先分析问题的解答方法，应该采用逆推的方式解决问题。每个水手起来的椰子数量等于前一个水手醒来时椰子数量少一的五分之四。最后每个水手得到的椰子数量会等于最后一堆椰子少一的五分之一。根据这个逆推，由于椰子的数量是整数，所以利用循环语句知道整数的解答为止。【程序实现】n=input(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，end end disp(x，p)运行这段程序后，屏幕出现要求从键盘输入 x 数据的信息input n，输入1200后，MATLAB计算出合适的 x 和 p0 的值为 1023 15621 【结果分析】要使得最初的椰子数p0为整数，必须取 (x +1) 为 4 5( =1024)的倍数，一种简单的处理可取 x = 1023。12设，（1）从尽可能精确的近似值出发，利用递推公式：计算机从到的近似值；（2）从较粗糙的估计值出发，用递推公式：计算从到的近似值；（3）分析所得结果的可靠性以及出现这种现象的原因【算法分析】 syms x n; int(1/(x+5),0,1) ans log(2)+log(3)-log(5)eval(ans)ans = 0.1823Now we use to compute In where n rang from 1 to 20.s=0.1823for n=1:20s=-5*s+1/n;disp(s)ends = 0.1823 0.0885 0.0575 0.0458 0.0208 0.0958 -0.3125 1.7054 -8.4018 42.1200 -210.5002 1.0526e+003 -5.2629e+003 2.6314e+004 -1.3157e+005 6.5786e+005 -3.2893e+006 1.6447e+007 -8.2233e+007 4.1116e+008 -2.0558e+009（2）首先计算I（30）的近似值 syms x n; int(x30/(x+5),0,1)eval (ans)0 /即I（30）0s=0;for n=30:-1:2s=-s/5+1/5*n;disp(s)ends = 0 6 4.6000 4.6800 4.4640 4.3072 4.1386 3.9723 3.8055 3.6389 3.4722 3.3056 3.1389 2.9722 2.8056 2.6389 2.4722 2.3056 2.1389 1.9722 1.8056 1.6389 1.4722 1.3056 1.1389 0.9722 0.8056 0.6389 0.4722 0.3056【结果分析】第一小题当中首先算出较为精确的I0，之后经过代入的方法计算出 I（20），结果是比较准确的。但是第二小题当中先是算出必是很精确的I（30），本来所得误差比第一步就稍微大一点了，然后再用这个误差稍微大的I（30）回代计算I（20）I（1），而从I（30）算到I（20）的时候误差这时已经够大了，采用这个误差更大的结果去计算我们想要的值，当然误差大得惊人啦，显然比第一种误差大得多。所以我们使用第一种方法得到的结果较接近准确值。实验二21用高斯消元法的消元过程作矩阵分解。设消元过程可将矩阵A化为上三角矩阵U，试求出消元过程所用的乘数、并以如下格式构造下三角矩阵L和上三角矩阵U验证：矩阵A可以分解为L和U的乘积，即A=LU。【算法分析】(1) The process of Gaussian elimination A= 20, 2,3;1,8,1;2,-3,15;B(1)=0,0,0;-1,-2/20,-3/20;0,0,0;C(1)=A+B(1)得到结果为：C(1) = 20.0000 2.0000 3.0000 0 7.9000 0.8500 2.0000 -3.0000 15.0000B(2)=0,0,0; 0,0,0; -20/10,-2/10,-3/10;C(2)=C(1)+B(2)得到结果为：C(2) = 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 -3.2000 14.7000B(3)=0,0,0; 0,0,0; 0,0.32,0.85*3.2/7.9;C(3)=B(3)+C(2)得到结果为：C(3)= 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443So U= 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443And from the Gaussion elimination we have thatL= 1 0 0 1/20 1 0 1/10 -3.2/7.9 1 验证：L=1，0，0；0.05，1，0；0.1，-0.405063，1; U=20，2，3；0，7.9，0.85；0，0，15.0443; A*Bans = 20.0000 2.0000 3.0000 1.0000 8.0000 1.0000 2.0000 -2.9995 14.9998 所以LU = 1 0 0 20.0000 2.0000 3.0000 1/20 1 0 0 7.9000 0.8500 1/10 -3.2/7.9 1 0 0 15.0443 L*U 20 2 3 1 8 1 2 -3 15【结果分析】因为在运算过程中遇到除不尽的分数采用的是四舍五入的方法来记录数据，所以所得结果是一个近似值，但是结果非常的接近准确值，所以说这种方法是正确的。22 用矩阵分解方法求上题中A的逆矩阵。记分别求解方程组由于三个方程组系数矩阵相同，可以将分解后的矩阵重复使用。对第一个方程组，由于A=LU，所以先求解下三角方程组,再求解上三角方程组，则可得逆矩阵的第一列列向量；类似可解第二、第三方程组，得逆矩阵的第二列列向量的第三列列向量。由三个列向量拼装可得逆矩阵。【程序实现】A = 20 2 3;1 8 1;2 -3 15;b1 = 1;0;0;b2 = 0;1;0;b3 = 0;0;1;L,U = lu(A);Y = inv(L)*b1;a1 = inv(U)*Y;Y = inv(L)*b2;a2 = inv(U)*Y;Y = inv(L)*b3;a3 = inv(U)*Y;a1 = 0.0517 -0.0055 -0.0080a2 = -0.0164 0.12370.0269a3 = -0.0093 -0.00720.0665【结果分析】所以A的逆矩阵是0.0517 -0.0164 -0.0093-0.0055 0.1237 -0.0072-0.0080 0.0269 0.0665同时，a1、a2、a3分别是三个方程组的解。实验三31 用泰勒级数的有限项逼近正弦函数用计算机绘出上面四个函数的图形。【算法分析】(1) y1=x,x 0,pi/2则算法如下： x=0:.1:pi/2; y1=sin(x); y2=x; plot(x,y1,x,y2)得到的图形如右图：（2）y(x)=sin x ,x0,；则y(x)x- x/3！具体算法如下： x=0:.1:pi;y1=sin(x);y2=x-x.3/6;plot(x,y1,x,y2)得到的图形如右图：（3）y3=x-x3/6+x5/120,x0，pi/2具体算法如下：x=0:.1:pi/2;y1=sin(x);y2=x-x.3/6+x.5/120;plot(x,y1,x,y2)得到的图形如右图：【结果分析】函数采用taloy展开项来逼近原函数，依次采用一阶、二阶、三阶的泰勒公式，显然有限项越多，逼近效果越好，从图形就可以很明显的看出来了。33追赶曲线的计算机模拟问题描述：欧洲文艺复兴时期的著名人物达芬奇曾经提出一个有趣的“狼追兔子”问题，当一只兔子正在它的洞穴南面60码处觅食时，一只饿狼出现在兔子正东的100码处。兔子急忙奔向自己的洞穴，狼立即以快于兔子一倍的速度紧追兔子不放。兔子一旦回到洞穴便逃脱厄，问狼是否会追赶上兔子？这一问题的研究方法可以推广到如鱼雷追击潜艇、地对空导弹击飞机等问题上去。在对真实系统做实验时，可能时间太长、费用太高、危险太大、甚至很难进行。计算机模拟是用计算机模仿实物系统，对系统的结构和行为进行动态演示，以评价或预测系统的行为效果。根据模拟对象的不同特点，分为确定性模拟和随机性模拟两大类。模拟通常所用的是时间步长法，即按照时间流逝的顺序一步一步对所研究的系统进行动态演示，以提取所需要的数据。问题分析：首先计算狼的初始位置到兔子洞穴的直线距离：由于狼奔跑的速度是兔子速度的两倍，兔子跑60码的时间狼可以跑120码。如果狼沿直线奔向兔窝，应该是可以追上兔子的。但是，有人推导出狼在追赶兔子过程中的运动曲线为根据曲线方程，当时，。也就是说，在没有兔窝的情况下兔子一直往北跑，在跑到大约66码处将被狼追上。由此可知，在有兔窝时狼是追赶不上兔子的。用计算机模拟的方法也可以得到同样的结论。取时间步长为1s，随时间步长的增加，考虑这一系统中的各个元素（狼和兔子）所处的位置变化规律，用计算机作出模拟。最后，根据第60s时狼所在的位置的坐标，判断狼是否能追上兔子。问题思考与实验：（1）设兔子奔跑的速度为，则狼运动的速度为。建立平面直角坐标系，若当时刻，兔子位于点处，狼位于点处。试根据，的坐标确定一个单位向量描述狼在时段内的运动方向。（2）根据狼的运动方向和速度推导到的坐标的具体表达式；（3）用计算机绘制追赶曲线的图形（包括静态和动态的图形）。【问题求解】 由于P点的运动方向始终指向Q,设在t=tk时刻的位置为Pk(xk,yk) ,则P点的运动方向可以由下面的单位向量表示Ek=PkQkPkQk = 0-xk(0-xk)2+(tk-yk)2，tk-yk(0-xk)2+(tk-yk)2P点的速度是2m/s，取时间步长为t ，设在t=tk+1时刻，P点的位置为Pk+1(xk+1,yk+1)，于是P点的位置变化规律为xk+1=xk+2t（0-xk）(0-xk)2+(tk-yk)2yk+1=yk+2t(tk-yk)(0-xk)2+(tk-yk)2【编程实现】建立M脚本文件 m3_3.mx2(1)=100;y2(1)=0;t(1)=0;tstep=1;for k=1:60 t(k+1)=t(k)+tstep; x2(k+1)=x2(k)+2*(0-x2(k)/sqrt(0-x2(k)2+(t(k)-y2(k)2); y2(k+1)=y2(k)+2*(t(k)-y2(k)/sqrt(0-x2(k)2+(t(k)-y2(k)2); x1(k)=0; y1(k)=t(k);endaxis(0,100,0,60)hold onfor k=1:61 plot(x2(k),y2(k),*,x1(k),y1(k),o) pause(0.5)endhold offd=sqrt(0-x(61)2+(t(61)-y(61)2);fprintf(狼与兔子的距离为%5.4n,d)运行程序可以出现动态的图形去掉语句pause(0.5)后，可以显示静态的图形如下图：解决本题时，使用了计算机模拟的方法，动态演示了系统的过程，使得问题的解决大大地简化了。实验四43 神经元模型用于蠓的分类识别（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 用 标记。观察图4-5 可以发现，代表 Af 的点位于图中偏右，而代表 Apf 的点位于图中偏左。应该存在一条位于两类点之间的直线 L， 作为 Af 和 Apf 分 界线，这条直线 L 应依据问题所给的数据，即学习样本来确定。设这条直线的方程为。对于平面上任意一点 P(x，y)，如果该点在直线上，将其坐标代入直线方程则使方程成为恒等式，即使方程左端恒为零；如果点 P(x，y) 不在直线上，将其坐标代入直线方程，则方程左端不为零。由于 Af 和 Apf 的散点都不在所求的直线上，故将问题所提供的数据代入直线方程左端所得到的表达式的值应有大于零或者小于零两种不同的结果。【算法分析】建立一个判别系统，引入判别函数 g(x，y)，当 (x，y) 代表 Af 类时，令g(x，y)0，否则 g(x，y)0。 为了对判别系统引入学习机制，在学习过程中将两种不同的状态，以“1”和“1”表示。取 则由所给数据形成约束条件，这是关于判别函数中的三个待定系数值 的线性方程组： 此为包括三个未知数共15个方程的超定方程组，可以求方程组的最小二乘解。 下列程序用于求超定方程组的最小二乘解，并绘制出分类边界曲线的图解。xy=1.24 1.27;1.36 1.74;1.38 1.64;1.38 1.82;1.38 1.90; 1.40 1.70;1.48 1.82;1.54 1.82;1.56 2.08;1.14 1.78; 1.18 1.96;1.20 1.86;1.26 2.00;1.28 2.00;1.30 1.96; %学习样本数据 z=1;1;1;1;1;1;1;1;1;-1;-1;-1;-1;-1;-1; x=xy(:,1);y=xy(:,2);x1=x(1:9);y1=y(1:9);x2=x(10:15);y2=y(10:15); plot(x1,y1,*,x2,y2,x),pause %绘制原始数据散点图 X=x y ones(x); A=X*X;B=X*z;w=AB; %求解正规方程组 a=w(1)/w(2);b=w(3)/w(2); t=1.10:0.02:1.60;u=a*x+b ; %确定分类直线数据 plot(x1,y1,*,x2,y2,x,t,u) %在散点图中画分类直线 运行上面程序可求出超定方程组的最小二乘解并画出分类边界曲线。为了由所给数据用判别函数判别三个新蠓虫的类属，即当 时，判为 Af 类；当 时，判为 Apf 类。运行上面程序后,键入下面命令 xx=1.24 1.80 1;1.28 1.84 1;1.40 2.04 1； xx*w plot(t，u，xx(:,1)，xx(:,2)，o) 求得 ans = -0.3877 -0.2384 -0.0235 这说明，所给数据反映出三个蠓虫均属于Apf类。【结果分析】上面的解决问题方案是在学习过程中用 +1 和 1 分别代表正数和负数来完成的。这只是一种人为的规定，并不是一成不变的。当Apf是害虫时，可以修改超定方程组的右端项中“-1”为“-0.6”，或者将 1 改为其它的负数以重新求超定方程组的最小二乘解获得分类边界直线的方程。这样将与前面所求分类边界直线的方程不一样，当然对新给定的蠓虫的翅膀和触角长度数据来做判断其结果也是不同的。实验五51 用几种不同的方法求积分的值。（1）牛顿-莱布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）复合梯形公式。【程序实现】function t=natrapz(fname,a,b,n)h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0.001*h);t=h*(0.5*(fa+fb)+sum(f) 由定步长梯形法求得值为 3function q,srm,err = naadapt( fname,a,b,e)srm=zeros(30,6);it=0;done=1;m=1;sr=simpson(fname,a,b,e);srm(1,1:6)=sr;state=it;while (state=it) n=m; for j=n:-1:1 p=j; sr0=srm(p,:); err=sr0(5);e=sr0(6); if e=err state=done; a=sr0(1);b=sr0(2);err=sr0(5);e=sr0(6); c=0.5*(a+b);e2=0.5*e; sr1=simpson(fname,a,c,e2); sr2=simpson(fname,c,b,e2); err=abs(sr0(3)-sr1(3)-sr2(3)/10;if err x,y=ode23(lore,4,5,0)x = 4.00000000000000 4.00000093280783 4.00000559684697 4.00002891704269 4.00014551802130 4.00072852291430 4.00364354737933 4.01821866970447 4.09109428133017 4.19109428133017 4.29109428133017 4.39109428133017 4.49109428133017 4.59109428133017 4.69109428133017 4.79109428133017 4.89109428133017 5.00000000000000y = 1.0e+002 * 0 0.00000080000028 0.00000480001007 0.00002480026893 0.00012480681052 0.00062497072747 0.00312907383321 0.01573209090479 0.08086238316915 0.17628097287779 0.27924859887657 0.39042430504800 0.51052612335396 0.64033652423016 0.78070837535945 0.93257145727939 1.09693958889058 1.29146855032127【结果分析】第二个积分的求解问题是通过转化为常微分初值问题来解决的，所得结果不是一个近似值，而是一个精确值，由此可知，积分和微分的求解问题是可以相互转化的，数学个问题之间是相互贯通的。64 列出函数在区间0，e上的函数值表并作出它的图象。其中，是初值问题的解。【问题分析】从题目中看出，In(v)在零点无定义，所以计算f：dfdx=dfdvdvdx=-Invdvdx=-2xf0=limv0(1-In(v)v=0【程序实现】 clear;simplify(dsolve(Dy=2*x/log(y),x) ans = (x2+2*C1)/lambertw(x2+2*C1)*exp(-1) fun=inline(2*x/log(y),x,y) ode23(fun,1,2.7,1)fun = Inline function: fun(x,y) = 2*x/log(y)Warning: Divide by zero.作图得：实验七7.4一个10次项式的系数为1 a1 a2 a9 a10=1 55 1320 181 50 157 773 902 055 341 693 0 -840 950 0 127 535 76 -106 286 40 632 880 0试用多项式的求根指令roots求出该10次方程的10个根，然后修改9次项的系数-55为-56，得新的10次方程，求解新的方程，观察根的变化是否很显著。【问题解答】此题应用roots(p)函数。p=1 -551320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;所以roots(p)= 551320.032920469 + 0.00000000000000i-2.24160771698781 + 0.00000000000000i-0.626975644192937 + 1.55805673614711i-0.626975644192937 - 1.55805673614711i1.21700668114928 + 0.00000000000000i0.820510882955133 + 0.840751943017165i0.820510882955133 - 0.840751943017165i0.302305044859346 + 0.994768021247990i0.302305044859346 - 0.994768021247990ip=1 -561320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;所以roots(p)= 561320.032333995 + 0.00