Matlab的实际应用设计(经典)

1、课程设计学院： 数学学院 学号： 20106496 姓名： 黄星奕 辅导老师： 陈晓红 殷明 题目一二三四五六七八总具体题目1.11.21.32.12.33.13.23.34.14.24.35.15.25.36.16.27.47.58.18.420题实验一11 水手、猴子和椰子问题一、 问题描述11 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一

2、只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？二、思考与实验试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。三、 问题分析用递推算法。首先分析椰子数目的变化规律，设最初的椰子数为p 0，即第一个水手所处理之前的椰子数，用p 1、p 2、p 3、p4、p 5 分别表示五个水手对椰子动了手脚以后剩余的椰子数目，则根据问题有 再用x表示最后每个水手平分得到的椰子数，于是有所以 p5 = 5x +1 利用逆向递推的方法，有 但由于椰子数为一正整数，用任意的x作为初值递推出的p0数据不一定是合适的。 在实验中可以用

3、for 循环语句结合 break 语句来寻找合适的 x 和 p0 ，对任意的 x 递推计算出 p0 ，当计算结果为正整数时，结果正确，否则选取另外的 x 再次重新递推计算，直到计算出的结果 p0 为正整数为止。四、源程序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; endenddisp(x,p);五、实验结果六、结果分析从理论上分析，由于 所以 要使得最初的椰子数p0为整数，必须取 (x +1) 为 4 5( =1024)的倍数，一种简单的处理可取

4、x = 1023。 12 一、问题描述12 设，（1）从尽可能精确的近似值出发，利用递推公式：计算机从到的近似值；（2）从较粗糙的估计值出发，用递推公式：计算从到的近似值；（3）分析所得结果的可靠性以及出现这种现象的原因。二、问题分析有种方法可以尽可能地精确的计算的值，我们可以根据积分计算得=ln 1.2=0.1823，然后可以编程求解到的近似值。三、源程序及运行结果(1) fun=inline('1./(5+x)','x')z=quad(fun,0,1)for n=1:20z=-5*z+1/nendz =0.0188z =0.0169z =0.0155z =0

5、.0135z =0.0156z =-0.0011z =0.0770z =-0.3186z =1.6554z =-8.2179z =41.1453z =-205.6737z =1.0284e+003（省略前7项）（2）fun=inline('(x.30)./(5+x)','x')z=quad(fun,0,1)for n=30:-1:2z=-0.2*z+1/(5*n)endz =0.0130z =0.0141z =0.0154z =0.0169z =0.0188z =0.0212z =0.0243z =0.0285z =0.0343z =0.0431z =0.058

6、0z =0.0884（限于篇幅省略前18项）四、结果分析第二种算法较为可靠，原因是迭代时系数较小，第一种方法虽然较精确，但后面的迭代系数绝对值为5，将误差逐步放大，所以最后结果反而不精确了。 13 绘制Koch分形曲线一、问题描述13 绘制Koch分形曲线：从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形（图1）；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形（图2），这时，图形中共有17个结点。这种迭代继续进行下去可以形成Koch分形曲线。在迭代过程中，图形中的结点将越来越多，而曲线最终

7、显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。Koch分形曲线的绘制与算法设计和计算机实现相关。二、思考与实验：（1）考虑在Koch分形曲线的形成过程中结点数目的变化规律。设第k次迭代产生结点数为，第迭代产生结点数为，试写出和之间的递推关系式；（2）参考问题分析中的算法，考虑图1到图2的过程，即由第一次迭代的5个结点的结点坐标数组，产生第二次迭代的17个结点的结点坐标数组的算法；（3）考虑由第k次迭代的个结点的结点坐标数组，产生第次迭代的个结点的结点坐标数组的算法；（4）设计算法用计算机绘制出如下的Koch分形曲线（图3）。 图 1 图 2图 3三、问题分析考虑由直线段（2个点）产生

8、第一个图形（5个点）的过程，设和分别为原始直线段的两个端点。现在需要在直线段的中间依次插入三个点产生第一次迭代的图形（图1）。显然，位于点右端直线段的三分之一处，点绕旋转60度（逆时针方向）而得到的，故可以处理为向量经正交变换而得到向量，形成算法如下：（1）；（2）；（3）；在算法的第三步中，A为正交矩阵。；这一算法将根据初始数据（和点的坐标），产生图1中5个结点的坐标。这5个结点的坐标数组，组成一个5×2矩阵。这一矩阵的第一行为为的坐标，第二行为的坐标，第二行为的坐标第五行为的坐标。矩阵的第一列元素分别为5个结点的x坐标 ，第二列元素分别为5个结点的y坐标。在Koch分形曲线的形成

9、过程中结点数目的变化规律。设第k次迭代产生结点数为n k，第k+1次迭代产生结点数为n k+1，n k和n k+1之间的递推关系式如下由第k次迭代的n k个结点的结点坐标数组，产生第k+1次迭代的n k+1个结点的结点坐标数组的算法可参考上面两点到五点的算法进行设计。 四、源程序p=0 0;10 0; %给出初始数据两个点的坐标 a=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); %设置用于正交变化的正交矩阵 for k=1:5 %开始执行第一到第五次迭代 n=max(size(p);d=diff(p)/3; %统计前一轮迭代的结点数及形成结点向量

10、 q=p(1:n-1,:);p(5:4:4*n-3,:)=p(2:n,:); %保护前一轮的结点坐标数组 p(2:4:4*n-6,:)=q+d; %插入第一组新结点 p(3:4:4*n-5,:)=p(2:4:4*n-6,:)+d*a'%用正交变换计算第二组新结点 p(4:4:4*n-4,:)=q+2*d; %插入第三组新结点 end plot(p(:,1),p(:,2) %根据结点坐标绘图 五、实验结果实验二21 用高斯消元法的消元过程作矩阵分解。一、問題描述21 用高斯消元法的消元过程作矩阵分解。设消元过程可将矩阵A化为上三角矩阵U，试求出消元过程所用的乘数、并以如下格式构

11、造下三角矩阵L和上三角矩阵U验证：矩阵A可以分解为L和U的乘积，即A=LU。二、源程序及運行結果%将矩阵A分解为L和U的乘积，即A=LU。a=input('a=');u=zeros(3,3);l=eye(3,3);for i=1:3 u(1,i)=a(1,i);endfor j=2:3 l(j,1)=a(j,1)/u(1,1);endfor k=2:3 u(k,k:3)=a(k,k:3)-l(k,1:k-1)*u(1:k-1,k:3); l(k+1:3,k)=(a(k+1:3,k)-l(k+1:3,1:k-1)*u(1:k-1,k)/u(k,k);enddisp(l);disp

12、(u);实验结果a=20,2,3;1,8,1;2,-3,15 1.0000 0 0 0.0500 1.0000 0 0.1000 -0.4051 1.0000 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443112.3验证希尔伯特矩阵的病态性一、问题描述2.3 验证希尔伯特矩阵的病态性：对于三阶矩阵取右端向量，验证：（1）向量是方程组的准确解；（2）取右端向量b的三位有效数字得，求方程组的准确解，并与X的数据作比较 。说明矩阵的病态性。二、源程序及运行结果(1) H=1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5b=11/6

13、 13/12 47/60'Hbans = 1.0000 1.0000 1.0000所以是方程组的准确解(2) b=1.83 1.08 0.783'b = 1.8300 1.0800 0.7830Hbans = 1.0800 0.5400 1.4400所算得的解误差较大，所以希尔伯特矩阵呈现病态性实验三31 用泰勒级数的有限项逼近正弦函数一、问题描述31 用泰勒级数的有限项逼近正弦函数用计算机绘出上面四个函数的图形。二、源程序图一x=0:0.1:pi/2;y1=x;y2=x-x.3/6;y3=x-x.3/6+x.5/120;plot(x,y1,x,y2,x,y3)legend(&

14、#39;y=x','y=x-x3/6','y=x-x3/6+x5/120')图二x=0:0.1:pi;y4=sin(x);plot(x,y4)；legend('y=sin(x)')三、实验结果 图一图二32 绘制飞机的降落曲线一、问题描述32 绘制飞机的降落曲线一架飞机飞临北京国际机场上空时，其水平速度为540km/h，飞行高度为1 000m。飞机从距机场指挥塔的横向距离12 000m处开始降落。根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点，则表示飞机距指挥塔的距离，表示飞机

15、的飞行高度，降落曲线为该函数满足条件：（1）试利用满足的条件确定三次多项式中的四个系数；（2）用所求出的三次多项式函数绘制出飞机降落曲线。二、问题分析由题设得a(0)=a(1)=0a(2)=1/48000a(3)=-1/(72000*12000)三、源程序fplot(inline('(1/4800)*x.2-1/(72000*12000)*x3'),0,12000)四、实验结果33 追赶曲线的计算机模拟一、问题描述33 追赶曲线的计算机模拟：问题描述：欧洲文艺复兴时期的著名人物达芬奇曾经提出一个有趣的“狼追兔子”问题，当一只兔子正在它的洞穴南面60码处觅食时，一只饿狼出现在兔子

16、正东的100码处。兔子急忙奔向自己的洞穴，狼立即以快于兔子一倍的速度紧追兔子不放。兔子一旦回到洞穴便逃脱厄，问狼是否会追赶上兔子？这一问题的研究方法可以推广到如鱼雷追击潜艇、地对空导弹击飞机等问题上去。在对真实系统做实验时，可能时间太长、费用太高、危险太大、甚至很难进行。计算机模拟是用计算机模仿实物系统，对系统的结构和行为进行动态演示，以评价或预测系统的行为效果。根据模拟对象的不同特点，分为确定性模拟和随机性模拟两大类。模拟通常所用的是时间步长法，即按照时间流逝的顺序一步一步对所研究的系统进行动态演示，以提取所需要的数据。二、思考与实验（1）设兔子奔跑的速度为，则狼运动的速度为。建立平面直角坐

17、标系，若当时刻，兔子位于点处，狼位于点处。试根据，的坐标确定一个单位向量描述狼在时段内的运动方向。（2）根据狼的运动方向和速度推导到的坐标的具体表达式；（3）用计算机绘制追赶曲线的图形（包括静态和动态的图形）。三、问题分析首先计算狼的初始位置到兔子洞穴的直线距离：由于狼奔跑的速度是兔子速度的两倍，兔子跑60码的时间狼可以跑120码。如果狼沿直线奔向兔窝，应该是可以追上兔子的。但是，有人推导出狼在追赶兔子过程中的运动曲线为根据曲线方程，当时，。也就是说，在没有兔窝的情况下兔子一直往北跑，在跑到大约66码处将被狼追上。由此可知，在有兔窝时狼是追赶不上兔子的。用计算机模拟的方法也可以得到同样的结论。

18、取时间步长为1s，随时间步长的增加，考虑这一系统中的各个元素（狼和兔子）所处的位置变化规律，用计算机作出模拟。最后，根据第60s时狼所在的位置的坐标，判断狼是否能追上兔子。设兔子所在位置为动点 Q，狼所在位置为动点P。在时刻 tk ，两个动点的坐标分别为：Q(uk, vk), P(xk, yk)，动点P的轨迹就是追赶曲线。在t k 时刻到t k+1 时刻这个时段，P点的运动方向可以用单位向量描述： 显然，uk = 0，vk = tk 。四、源程序设在初始时刻兔子和狼的位置分别为 Q(0，0) ，P(100，0) 初始时刻的狼、兔距离为100码，我们不妨规定当狼、兔距离小于0.5码时，兔子被狼追

19、上，结束追赶。下面MATLAB程序可计算并绘制追赶曲线。 x(1)=100;y(1)=0;u(1)=0;v(1)=0;t=1;d=100;e=-1 0; while d>0.5 x(t+1)=x(t)+2*e(1); y(t+1)=y(t)+2*e(2); t=t+1;u(t)=0;v(t)=t; e=-x(t) t-y(t); d=sqrt(e(1)2+e(2)2);e=e/d; end plot(u,v,'o',x,y)五、实验结果实验四4.1线性拟合一、问题描述曾任英特尔公司董事长的摩尔先生早在1965年时，就观察到一件很有趣的现象：集成电路上可容纳的零件数量，每隔

20、一年半左右就会增长一倍，性能也提升一倍。因而发表论文，提出了大大有名的摩尔定律（Moores Law），并预测未来这种增长仍会延续下去。下面数据中，第二行数据为晶片上晶体数目在不同年代与1959年时数目比较的倍数。这些数据是推出摩尔定律的依据：年代19591962196319641965增加倍数13456试从表中数据出发，推导线性拟合的函数表达式。二、源程序x=0,3,4,5,6;y=1,3,4,5,6; p=polyfit(x,y,1)p = 0.8302 0.8113xi=0:0.2:6;yi=polyval(p,xi); plot(x,y,'o',xi,yi);五、实验结

21、果42 问题描述：参考算法4.2设计绘制Bezier曲线的程序，选取四个点的坐标数据作为控制点绘制飞机机翼剖面图草图的下半部分图形；结合例4.4中上半部分图形绘出完整的机翼草图。最后写出机翼剖面图曲线上20个点处的坐标数据。源代码及结果：p=50 -50;60 -60;100 -80;150 -60;200 -40;n=size(p,1);t=linspace(0,1)'b=0;for k=0:n-1 tmp=nchoosek(n-1,k)*t.k.*(1-t).(n-1-k); b=b+tmp*p(k+1,:);endplot(p(:,1),p(:,2),'.:',b

22、(:,1),b(:,2)43 神经元模型用于蠓的分类识别一、问题描述问题描述：生物学字试图对两类蠓虫（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）依据确立的方法，对题目提

23、供的三个样本：（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），如果该点在直线上，将其坐标代入直线方程则使方程成为恒等式，即使方程左端恒为零；如果点不在直线上，将其

24、坐标代入直线方程，则方程左端不为零。由于f和pf的散点都不在所求的直线上，故将问题所提供的数据代入直线方程左端应该得到表达式的值大于零或者小于零两种不同的结果。这需要建立一个判别系统，引入判别函数,当属于f类时，否则。为了对判别系统引入学习机制，在学习过程中将两种不同的状态，以“”和“”表示。当属于f类时，否则。取于是由所给数据形成约束条件，这是关于判别函数中的三个待定系数的线性方程组：这是包括三个未知数共15个方程的超定方程组，可以求方程组的最小二乘解。三、思考与实验：（1）根据上面分析写出对应的正规方程组并求解。（2）确定分类边界直线的方程。由所给数据用判别函数判别三个新蠓虫的类属，即当时

25、，判为Af类：当时，判为Apf类。四、源程序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

26、,'*',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 类；当

27、 时，判为 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类。 实验五5.1用几种不同的方法求积分的值。一、问题描述51 用几种不同的方法求积分的值。（1）牛顿-莱布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）复合梯形公式。二、源程序及运行结果fun=inline('4./(1+x.2)','x&#

28、39;)fun = Inline function: fun(x) = 4./(1+x.2)（1）牛顿-莱布尼茨公式syms x;z1=int(4./(1+x.2),x,0,1)z1= piz1 =piz1 =3.1416（2）梯形公式z2=1/2*fun(0)+fun(1)z2 = 3(3) 辛卜生公式z3=1/6*fun(0)+4*fun(0.5)+fun(1)z3 =3.1333(4) 复合梯形公式z4=0.05*fun(0)+2*fun(0.1)+2*fun(0.2)+2*fun(0.3)+2*fun(0.4)+2*fun(0.5)+2*fun(0.6)+2*fun(0.7)+2*fu

29、n(0.8)+2*fun(0.9)+fun(1)z4 =3.1399或者clear;x=0:0.1:1;y=4./(1+x.2);trapz(x,y)ans =3.13995.2 设计算法计算30个不同的概率值一、问题描述52 设X为标准正态随机变量，即XN（0，1）。现分别取，试设计算法计算30个不同的概率值；，并将计算结果与概率论教科书中的标准正态分布函数表作比较。（提示： 二、源程序fun=inline('exp(-x.2/2)/sqrt(2*pi)');a=0;h=0.1;for k=1:30a=a+h;p=quad(fun,a,4)end三、实验结果p =0.0178

30、p = 0.0139p = 0.0107p = 0.0082p = 0.0062p = 0.0046p = 0.0034p = 0.0025p = 0.0018p = 0.0013p = 0.1356p = 0.1150p = 0.0968p = 0.0807p = 0.0668p = 0.0548p = 0.0445p = 0.0359p = 0.0287p = 0.0227p = 0.4601p = 0.4207p = 0.3821p = 0.3445p = 0.3085p = 0.2742p = 0.2419p = 0.2118p = 0.1840p = 0.158653 设某城市男子的

31、身高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）。如果将汽车门的高

32、度取180cm，是否满足大多数市民的利益？（3）用计算机模拟的方法来检验你的结论，计算机产生10 000个正态随机数（它们服从均值为170，方差为6的正态分布）来模拟这个城市中10 000个男子的身高，然后统计出这10 000人中身高超过180（cm）的男子数量所占的百分比。解：程序如下：（1）fun=inline('exp(-1*(x-170).2)/(2*36)/(6*pi*sqrt(2)');for k=180:188p=quad(fun,k,194)end (2)结果如下：H=180、181、188时定积分近似值如下：p =0.0269p =0.0188p =0.012

33、8p =0.0085p =0.0055p =0.0035p =0.0021p =0.0013p =7.4373e-004p =0.0269p =0.0188p =0.0128p =0.0085p =0.0055p =0.0035p =0.0021p = 0.0013p =7.4373e-004 (3)2，程序如下：n=10000;r=170+6*randn(n,1);s=0;for i=1:n if r(i)>180 s=s+1; endendpp=s/n3、结果如下：pp =0.0486实验六61 用欧拉公式和四阶龙格-库塔法分别求解下列初值问题；一、问题描述61 用欧拉公式和四阶龙格

34、-库塔法分别求解下列初值问题；二、源程序及运行结果(1)a.欧拉法function euler(a,b,h,alpha)t=a:h:b;n=(b-a)/h;y(1)=alpha;for i=2:n+1 y(i)=y(i-1)+h*f(t(i-1),y(i-1); sprintf('t=%3.1f,y=%9.7f',t(i),y(i)endfunction r=f(t,y)r=(0.9*y)/(1+2*t);在Mathlab Command Window里输入如下代码:euler(0,1,0.1,1)得到：ans =t=0.1,y=1.0900000ans =t=0.2,y=1.

35、1717500ans =t=0.3,y=1.2470768ans =t=0.4,y=1.3172249ans =t=0.5,y=1.3830861ans =t=0.6,y=1.4453250ans =t=0.7,y=1.5044519ans =t=0.8,y=1.5608688ans =t=0.9,y=1.6148989ans =t=1.0,y=1.6668064function r=f(x,y)r=-x*y/(1+x*x);在Mathlab Command Window里输入如下代码euler(0,1,0.1,2)结果：ans =t=0.1,y=2.0000000ans =t=0.2,y=1

36、.9801980ans =t=0.3,y=1.9421173ans =t=0.4,y=1.8886645ans =t=0.5,y=1.8235382ans =t=0.6,y=1.7505966ans =t=0.7,y=1.6733644ans =t=0.8,y=1.5947500ans =t=0.9,y=1.5169573ans =t=1.0,y=1.4415285b. 四阶龙格-库塔法fun=inline('0.9*y/(1+2*t)','t','y');t,y=ode45(fun,0,1,1);t,yans = 0 1.0000 0.0250

37、 1.0222 0.0500 1.0438 0.0750 1.0649 0.1000 1.0855 0.1250 1.1056 0.1500 1.1253 0.1750 1.1446 0.2000 1.1635 0.2250 1.1820 0.2500 1.2002 0.2750 1.2180 0.3000 1.2355 0.3250 1.2528 0.3500 1.2697 0.3750 1.2864 0.4000 1.3028 0.4250 1.3189 0.4500 1.3349 0.4750 1.3506 0.5000 1.3660 0.5250 1.3813 0.5500 1.396

38、4 0.5750 1.4112 0.6000 1.4259 0.6250 1.4404 0.6500 1.4547 0.6750 1.4689 0.7000 1.4828 0.7250 1.4967 0.7500 1.5103 0.7750 1.5239 0.8000 1.5372 0.8250 1.5505 0.8500 1.5636 0.8750 1.5765 0.9000 1.5894 0.9250 1.6021 0.9500 1.61470.9750 1.62711.0000 1.6395plot(t,y,'o-') %画图（2）odefun=inline('-

39、x*y/(1+x*x)','x','y');x,y=ode45(odefun,0,1,1);x,yans = 0 1.0000 0.0250 0.9997 0.0500 0.9988 0.0750 0.9972 0.1000 0.9950 0.1250 0.9923 0.1500 0.9889 0.1750 0.9850 0.2000 0.9806 0.2250 0.9756 0.2500 0.9701 0.2750 0.9642 0.3000 0.9578 0.3250 0.9510 0.3500 0.9439 0.3750 0.9363 0.4000

40、 0.9285 0.4250 0.9203 0.4500 0.9119 0.4750 0.9033 0.5000 0.8944 0.5250 0.8854 0.5500 0.8762 0.5750 0.8669 0.6000 0.8575 0.6250 0.8480 0.6500 0.8384 0.6750 0.8288 0.7000 0.8192 0.7250 0.8096 0.7500 0.8000 0.7750 0.7904 0.8000 0.7809 0.8250 0.7714 0.8500 0.7619 0.8750 0.7526 0.9000 0.7433 0.9250 0.734

41、1 0.9500 0.7250 0.9750 0.7160 1.0000 0.7071>> plot(x,y,'o-')62 用求解常微分方程初值问题的方法计算积分上限函数一、问题描述62 用求解常微分方程初值问题的方法计算积分上限函数的值，实际上是将上面表达式两端求导化为常微分方程形式，并用初值条件。试用MATLAB中的指令ode23解决定积分的计算问题。二、源程序及运行结果odefun=inline('2*pi*exp(x)/x','x','y');ode23t(odefun,4,5,0);用梯形法则算得精确解x=

42、4:0.1:5;y= 2*pi*exp(x)./x;trapz(x,y)ans = 129.2178实验七7.4观察根的变化是否很显著一、问题描述7.4一个10次项式的系数为1 a1 a2 a9 a10=1 55 1320 18150 157773 902055 3416930 -8409500 12753576 -10628640 6328800试用多项式的求根指令roots求出该10次方程的10个根，然后修改9次项的系数-55为-56，得新的10次方程，求解新的方程，观察根的变化是否很显著。p=1 -56 1320 -18150 157773 -902055 3416930 -840950

43、0 12753576 -10628640 6328800;>> x=roots(p)x = 21.7335 7.3521 + 7.8018i 7.3521 - 7.8018i 3.6433 + 3.5988i 3.6433 - 3.5988i 4.9248 + 1.8070i 4.9248 - 1.8070i 1.1808 + 2.0164i 1.1808 - 2.0164i 0.0643 二、源程序及运行结果p=1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;>>

44、 x=roots(p) x = 10.6051 + 1.0127i 10.6051 - 1.0127i 8.5850 + 2.7898i 8.5850 - 2.7898i 5.5000 + 3.5058i 5.5000 - 3.5058i 2.4150 + 2.7898i 2.4150 - 2.7898i 0.3949 + 1.0127i 0.3949 - 1.0127i7.5 将三村短路问题推广为四村短路问题，即已知四个点A，B，C，D的具体位置A（0，0，），B(0,3),C(8,1),D(10,5),求两个点H1，H2的具体位置，使AH1+BH1+H1H2+H2C+H2D为最短。解：设所

45、要求的点构造函数 要求使 S 取得最小值的 即必须满足方程组 化简即得 解得 这个解也即使AH1+BH1+H1H2+H2C+H2D为最短的解。实验八8.1分别用直接法、雅可比迭代法、赛德尔迭代法求解线性方程组AX=b。（1），（2），源代码及结果：（1）直接法：e=ones(10,1);a=spdiags(e,4*e,e,-1,0,1,10,10);b=2;1;1;1;1;1;1;1;1;2;x=transpose(ab)x = 0.4793 0.0829 0.1891 0.1608 0.1678 0.1678 0.1608 0.1891 0.0829 0.4793雅克比：e=ones(10,

46、1);a=spdiags(e,4*e,e,-1,0,1,10,10);b=2;1;1;1;1;1;1;1;1;2;x=0;0;0;0;0;0;0;0;0;0;y=0;0;0;0;0;0;0;0;0;0;for k=1:8 eorr=0; for i=1:10 s=x(i);x(i)=0; y(i)=(b(i)-a(i,:)*x)/a(i,i); eorr=max(abs(s-x(i),eorr); x=y; endx',pause,eorrendans =0.5000 0.1250 0.2188 0.1953 0.2012 0.1997 0.2001 0.2000 0.2000 0.4

47、500赛德尔：e=ones(10,1);a=spdiags(e,4*e,e,-1,0,1,10,10);b=2;1;1;1;1;1;1;1;1;2;x=0;0;0;0;0;0;0;0;0;0;for k=1:8 eorr=0; for i=1:10 s=x(i);x(i)=0; x(i)=(b(i)-a(i,:)*x)/a(i,i); eorr=max(abs(s-x(i),eorr); endx',pause,eorrendans =0.4688 0.0781 0.1816 0.1543 0.1615 0.1596 0.1601 0.1600 0.0975 0.4756（2）直接法：

48、e=ones(20,1);a=spdiags(e,-2*e,5*e,-2*e,e,-2,-1,0,1,2,20,20);b=1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;x=transpose(ab)x = Columns 1 through 12 0.2421 0.0944 -0.0218 -0.0314 -0.0069 0.0049 0.0036 0.0002 -0.0008 -0.0004 0.0001 0.0001 Columns 13 through 20 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000雅克比：e=ones(20,1);a=spdiags(e,-2*e,5*e,-2*e,e,-2,-1,0,1,2,20,20);b=1;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;x=0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;y=0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;for k=1:8 eorr=0; for i=1:20 s=x(i);x(i)