Matlab的实际应用设计（经典）

课 程 设 计

学院： 数学学院

学号： 20106496

姓名： 黄星奕

辅导老师： 陈晓红 殷明

题目 一 具体题目 实验一

1.1 1.2 1.3 二 2.1 2.3 三 3.1 3.2 3.3 四 4.1 4.2 4.3 五 5.1 5.2 5.3 六 6.1 6.2 七 7.4 7.5 八 8.1 8.4 总 20题 1．1 水手、猴子和椰子问题 一、 问题描述

1．1 水手、猴子和椰子问题：五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。由于旅途的颠簸，大家都很疲惫，很快就入睡了。第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？ 二、思考与实验

试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。 三、 问题分析

用递推算法。首先分析椰子数目的变化规律，设最初的椰子数为p 0，即第一个水手所处理之前的椰子数，用p 1、p 2、p 3、p4、p 5 分别表示五个水手对椰子动了手脚以后剩余的椰子数目，则根据问题有

再用x表示最后每个水手平分得到的椰子数，于是有

所以

p5 = 5x +1

利用逆向递推的方法，有

但由于椰子数为一正整数，用任意的x作为初值递推出的p0数据不一定是合适的。 在实验中可以用 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; end end

disp([x,p]);

五、实验结果

六、结果分析

从理论上分析，由于

所以

要使得最初的椰子数p0为整数，必须取 (x +1) 为 4 5( =1024)的倍数，一种简单的处理可取 x = 1023。

xndx 1．2 In??05?x1一、问题描述

xndx 1．2 设，In??05?x1（1）从I0尽可能精确的近似值出发，利用递推公式：

1In??5In?1?(n?1,2,n计算机从I1到I20的近似值；

（2）从I30较粗糙的估计值出发，用递推公式：

20)

11In?1??In?(n?30,29,55n计算从I1到I20的近似值；

,3,2)

（3）分析所得结果的可靠性以及出现这种现象的原因。

二、问题分析

有种方法可以尽可能地精确的计算I0的值，我们可以根据积分计算得I0=ln 1.2=0.1823，然后可以编程求解I1到I20的近似值。 三、源程序及运行结果 (1) fun=inline('1./(5+x)','x') z=quad(fun,0,1) for n=1:20 z=-5*z+1/n end

z =0.0188 z =0.0169 z =0.0155 z =0.0135

z =0.0156 z =-0.0011 z =0.0770 z =-0.3186 z =1.6554 z =-8.2179 z =41.1453 z =-205.6737

z =1.0284e+003（省略前7项）

（2）fun=inline('(x.^30)./(5+x)','x') z=quad(fun,0,1) for n=30:-1:2 z=-0.2*z+1/(5*n) end

z =0.0130 z =0.0141 z =0.0154 z =0.0169 四、结果分析

z =0.0188 z =0.0212 z =0.0243 z =0.0285 z =0.0343 z =0.0431 z =0.0580 z =0.0884

（限于篇幅省略前18项）

第二种算法较为可靠，原因是迭代时系数较小，第一种方法虽然I0较精确，但 后面的迭代系数绝对值为5，将误差逐步放大，所以最后结果反而不精确了。 1．3 绘制Koch分形曲线

一、问题描述

1．3 绘制Koch分形曲线：从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形（图1）；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形（图2），这时，图形中共有17个结点。这种迭代继续进行下去可以形成Koch分形曲线。在迭代过程中，图形中的结点将越来越多，而曲线最终显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。Koch分形曲线的绘制与算法设计和计算机实现相关。

二、思考与实验：

（1）考虑在Koch分形曲线的形成过程中结点数目的变化规律。设第k次迭代产生结点数为nk，第k?1迭代产生结点数为nk?1，试写出nk和nk?1之间的递推关系式； （2）参考问题分析中的算法，考虑图1到图2的过程，即由第一次迭代的5个结点的结点坐标数组，产生第二次迭代的17个结点的结点坐标数组的算法；

（3）考虑由第k次迭代的nk个结点的结点坐标数组，产生第k?1次迭代的nk?1个结点的结点坐标数组的算法；

（4）设计算法用计算机绘制出如下的Koch分形曲线（图3）。

图 1 图 2

图 3

三、问题分析

①考虑由直线段（2个点）产生第一个图形（5个点）的过程，设P1和P5分别为原始直线段的两个端点。现在需要在直线段的中间依次插入三个点P2,P3,P4产生第一次迭代的图形（图1）。显然，P2位于P1点右端直线段的三分之一处，P4点绕P2旋转60度（逆时针方向）而得到的，故可以处理为向量P2P4经正交变换而得到向量P2P3，形成算法如下：

（1）P2?P1?(P5?P1)/3； （2）P4?P1?2(P5?P1)/3；

T（3）P； ?P?(P?P)?A3242在算法的第三步中，A为正交矩阵。

??cos?3A???sin??3?3??； ??cos3???sin??这一算法将根据初始数据（P1和P5点的坐标），产生图1中5个结点的坐标。这5个结点的坐标数组，组成一个5×2矩阵。这一矩阵的第一行为为P1的坐标，第二行为P1的坐标，第二行为P2的坐标……第五行为P5的坐标。矩阵的第一列元素分别为5个结点的x坐标 ，第二列元素分别为5个结点的y坐标。

②在Koch分形曲线的形成过程中结点数目的变化规律。设第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; %统计前一轮迭代的结点数及形成结点向量 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)) %根据结点坐标绘图 五、实验结果

实验二

2．1 用高斯消元法的消元过程作矩阵分解。 一、問題描述

2．1 用高斯消元法的消元过程作矩阵分解。设

?2023??

A??181????2?315??消元过程可将矩阵A化为上三角矩阵U，试求出消元过程所用的乘数m21、m31、m31并以

如下格式构造下三角矩阵L和上三角矩阵U

?1??202?,U??(1)L??m1a2122??????m31m321???验证：矩阵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); end

for j=2:3

l(j,1)=a(j,1)/u(1,1); end

for 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); end disp(l); disp(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.044311

2.3验证希尔伯特矩阵的病态性 一、问题描述

2.3 验证希尔伯特矩阵的病态性：对于三阶矩阵

3?(1)?a23? (2)?a33?1/21/3??1?H??1/21/31/4????1/31/41/5??

Tb?[11/613/1247/60]取右端向量，验证：

X?[x1x2（1）向量

x3]T?[111]T是方程组HX?b的准确解；

T?b?[1.831.080.783]（2）取右端向量b的三位有效数字得，求方程组的准确解X，

T[111]并与X的数据作比较 。说明矩阵的病态性。

二、源程序及运行结果

(1) H=[1 1/2 1/3;1/2 1/3 1/4;1/3 1/4 1/5]

b=[11/6 13/12 47/60]' H\\b ans =

1.0000 1.0000 1.0000

X?[x1所以

x2x3]T?[111]T是方程组HX?b的准确解

(2) b=[1.83 1.08 0.783]' b =

1.8300 1.0800 0.7830 H\\b ans =

1.0800 0.5400 1.4400

所算得的解误差较大，所以希尔伯特矩阵呈现病态性

实验三

3．1 用泰勒级数的有限项逼近正弦函数 一、问题描述

3．1 用泰勒级数的有限项逼近正弦函数

y0(x)?sinx,x?[0,?]y1(x)?x,x?[0,?/2]y2(x)?x?x/6,x?[0,?/2]y3(x)?x?x3/6?x5/120,x?[0,?/2]用计算机绘出上面四个函数的图形。 二、源程序

图一

3

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('y=x','y=x-x^3/6','y=x-x^3/6+x^5/120')

图二

x=0:0.1:pi; y4=sin(x); plot(x,y4)；

legend('y=sin(x)')

三、实验结果 图一

图二

3．2 绘制飞机的降落曲线

一、问题描述

3．2 绘制飞机的降落曲线

一架飞机飞临北京国际机场上空时，其水平速度为540km/h，飞行高度为1 000m。飞机从距机场指挥塔的横向距离12 000m处开始降落。根据经验，一架水平飞行的飞机其降落曲线是一条三次曲线。建立直角坐标系，设飞机着陆点为原点O，降落的飞机为动点

P(x,y)，则x表示飞机距指挥塔的距离，y表示飞机的飞行高度，降落曲线为

y(x)?a0?a1x?a2x2?a3x3

该函数满足条件：

y(0)?0,y(12000)?1000

y'(0)?0,y'(12000)?0（1）试利用y(x)满足的条件确定三次多项式中的四个系数； （2）用所求出的三次多项式函数绘制出飞机降落曲线。 二、问题分析 由题设得a(0)=a(1)=0 a(2)=1/48000

a(3)=-1/(72000*12000)

三、源程序

fplot(inline('(1/4800)*x.^2-1/(72000*12000)*x^3'),[0,12000])

四、实验结果

3．3 追赶曲线的计算机模拟 一、问题描述

3．3 追赶曲线的计算机模拟：问题描述：欧洲文艺复兴时期的著名人物达?芬奇曾经提出一个有趣的“狼追兔子”问题，当一只兔子正在它的洞穴南面60码处觅食时，一只饿狼出现在兔子正东的100码处。兔子急忙奔向自己的洞穴，狼立即以快于兔子一倍的速度紧追兔子不放。兔子一旦回到洞穴便逃脱厄，问狼是否会追赶上兔子？

这一问题的研究方法可以推广到如鱼雷追击潜艇、地对空导弹击飞机等问题上去。 在对真实系统做实验时，可能时间太长、费用太高、危险太大、甚至很难进行。计算机模拟是用计算机模仿实物系统，对系统的结构和行为进行动态演示，以评价或预测系统的行为效果。根据模拟对象的不同特点，分为确定性模拟和随机性模拟两大类。模拟通常所用的是时间步长法，即按照时间流逝的顺序一步一步对所研究的系统进行动态演示，以提取所需要的数据。 二、思考与实验

（1）设兔子奔跑的速度为?0?1m/s，则狼运动的速度为?1?2?0。建立平面直角坐标系，若当t?tk时刻，兔子位于点Qk(0,?0tk)处，狼位于点P试根据Pk，Qkk(xk,yk)处。的坐标确定一个单位向量ek描述狼在[tk,tk?1]时段内的运动方向。

（2）根据狼的运动方向和速度推导P k(xk,yk)到Pk?1(xk?1,yk?1)的坐标的具体表达式；（3）用计算机绘制追赶曲线的图形（包括静态和动态的图形）。 三、问题分析

首先计算狼的初始位置到兔子洞穴的直线距离：

D?1002?602?116.6190

由于狼奔跑的速度是兔子速度的两倍，兔子跑60码的时间狼可以跑120码。如果狼沿直线奔向兔窝，应该是可以追上兔子的。但是，有人推导出狼在追赶兔子过程中的运动曲线为

1132002y(x)?x?10x2?

303根据曲线方程，当x?0时，y?200/3。也就是说，在没有兔窝的情况下兔子一直往北跑，在跑到大约66码处将被狼追上。由此可知，在有兔窝时狼是追赶不上兔子的。

用计算机模拟的方法也可以得到同样的结论。取时间步长为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码时，兔子被狼追上，结束追赶。下面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年时，就观察到一件很有趣的现象：集成电路上可容纳的零件数量，每隔一年半左右就会增长一倍，性能也提升一倍。因而发表论文，提出了大大有名的摩尔定律（Moore’s Law），并预测未来这种增长仍会延续下去。下面数据中，第二行数据为晶片上晶体数目在不同年代与1959年时数目比较的倍数。这些数据是推出摩尔定律的依据： 年代 增加倍数 1959 1 1962 3 1963 4 1964 5 1965 6 试从表中数据出发，推导线性拟合的函数表达式。 二、源程序

x=[0,3,4,5,6]; y=[1,3,4,5,6]; p=polyfit(x,y,1) p =

0.8302 0.8113 xi=0:0.2:6;

yi=polyval(p,xi); plot(x,y,'o',xi,yi); 五、实验结果

4．2

问题描述：

参考算法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,:); end

plot(p(:,1),p(:,2),'.:',b(:,1),b(:,2))

4．3 神经元模型用于蠓的分类识别 一、问题描述

问题描述：生物学字试图对两类蠓虫（Af与Apf）进行鉴别，依据的资料是蠓虫的触角和翅膀的长度，已经测得9只Af和6只Apf的数据（触角长度用x表示，翅膀长度用y表示）

Af数据 x Y 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

Apf数据 x Y 1.14 1.18 1.20 1.26 1.28 1.30 1.78 1.96 1.86 2.00 2.00 1.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的确定应依据问题所给的数据，即学习样本。设这条直线的方程为

?1x??2y??0?0

对于平面上任意一点P（x,y），如果该点在直线上，将其坐标代入直线方程则使方程成为恒等式，即使方程左端恒为零；如果点P(x,y)不在直线上，将其坐标代入直线方程，则方程左端不为零。由于Ａf和Ａpf的散点都不在所求的直线上，故将问题所提供的数据代入直线方程左端应该得到表达式的值大于零或者小于零两种不同的结果。

这需要建立一个判别系统，引入判别函数g(P(x,y)),当P(x,y)属于Ａf类时，

g(P(x,y))?0，否则g(P(x,y))?0。

为了对判别系统引入学习机制，在学习过程中将两种不同的状态，以“１”和“－１”表示。当P(x,y)属于Ａf类时，g(P(x,y))?1，否则g(P(x,y))??1。 取

g(P(x,y))??1x??2y??0

于是由所给数据形成约束条件，这是关于判别函数中的三个待定系数?0,?1,?2的线性方程组：

??1xj??2yj??0?1,(j?1,2,,9) ??x??y????1,(j?10,,15)2j0?1j这是包括三个未知数共15个方程的超定方程组，可以求方程组的最小二乘解。 三、思考与实验：

（1）根据上面分析写出对应的正规方程组并求解。

（2）确定分类边界直线的方程。由所给数据用判别函数判别三个新蠓虫的类属，即当

g(P(x,y))?0时，判为Af类：当g(P(x,y))?0时，判为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,'*',x2,y2,'x'),pause %绘制原始数据散点图 X=[x y ones(x)];

A=X'*X;B=X'*z;w=A\\B; %求解正规方程组 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) %在散点图中画分类直线 五，实验结果

六、结果分析

运行上面程序可求出超定方程组的最小二乘解并画出分类边界曲线。为了由所给数据用判别函数判别三个新蠓虫的类属，即当

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') 求得

时，判为 Af 类；当

时，判为

ans = -0.3877 -0.2384 -0.0235

这说明，所给数据反映出三个蠓虫均属于Apf类。

实验五 5.1用几种不同的方法求积分一、问题描述

5．1 用几种不同的方法求积分

4?01?x2dx的值。

14?01?x2dx的值。

1（1）牛顿-莱布尼茨公式；（2）梯形公式；（3）辛卜生公式；（4）复合梯形公式。。 二、源程序及运行结果

fun=inline('4./(1+x.^2)','x') fun =

Inline function: fun(x) = 4./(1+x.^2) （1）牛顿-莱布尼茨公式 syms x;

z1=int(4./(1+x.^2),x,0,1) z1= pi z1 = pi z1 = 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*fun(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.1399

5.2 设计算法计算30个不同的概率值

一、问题描述

5．2 设X为标准正态随机变量，即X~N（0，1）。现分别取u?0.1,0.2,0.3,,2.9,3.0，

试设计算法计算30个不同的概率值；P?X?ua?，并将计算结果与概率论教科书中的标准正态分布函数表作比较。

1（提示： P?X?ua??2?二、源程序

???uae?x221dx?2??4uae?x22dx)

fun=inline('exp(-x.^2/2)/sqrt(2*pi)'); a=0;h=0.1; for k=1:30 a=a+h;

p=quad(fun,a,4) end

三、实验结果 p = p = 0.4601 0.1356 p = p = 0.4207 0.1150 p = p = 0.3821 0.0968 p = p = 0.3445 0.0807 p = p = 0.3085 0.0668 p = p = 0.2742 0.0548 p = p = 0.2419 0.0445 p = p = 0.2118 0.0359 p = p = 0.1840 0.0287 p = p = 0.1586 0.0227

5．3 设某城市男子的身高X~N（170，36）（单位：cm），应如何选择公共汽车门的高度H使男子与车门碰头的机会小于1%。

问题分析：由题设男子身高数据服从平均值为170（cm），方差为6（cm）的正态分布，其分布密度函数为

p = 0.0178 p =

0.0139 p =

0.0107 p =

0.0082 p =

0.0062 p =

0.0046 p =

0.0034 p =

0.0025 p =

0.0018 p =

0.0013

??(x?170)2?1f(x)?exp??62??2?36?

按正态分布的分布规律（3?原则），这个城市的男子身高超过188（cm）的人数极少。故可以对H=188，187，186，…求出概率

??194P?X?H?的值，观察使概率不超过1%的H，以

确定公共汽车门应该取的高度。概念值的计算实际上是求定积分

（1）选用一种数值求积公式或用数学软件分别计算出H=180、181、…、188时定积分近似值。

HHP?X?H???f(x)dx??f(x)dx（2）根据上面计算的积分值，按题目要求确定公共汽车门的高度取值（答案184cm）。如果将汽车门的高度取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:188 p=quad(fun,k,194) end (2)

结果如下：

H=180、181、…、188时定积分近似值如下： p =0.0269 p =0.0188 p =0.0128 p =0.0085 p =0.0055 p =0.0035 p =0.0021 p =0.0013

p =7.4373e-004 p =0.0269 p =0.0188 p =0.0128 p =0.0085 p =0.0055 p =0.0035 p =0.0021 p = 0.0013

p =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; end end pp=s/n

3、结果如下：pp =0.0486

实验六

6．1 用欧拉公式和四阶龙格-库塔法分别求解下列初值问题； 一、问题描述

6．1 用欧拉公式和四阶龙格-库塔法分别求解下列初值问题；

0.9y,y(0)?1;x?[0,1]1?2x

xy(2)y'??,y(0)?2;x?[0,1]21?x(1)y'?二、源程序及运行结果

(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)) end

①function 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.0900000 ans =

t=0.2,y=1.1717500 ans =

t=0.3,y=1.2470768 ans =

t=0.4,y=1.3172249 ans =

t=0.5,y=1.3830861 ans =

t=0.6,y=1.4453250 ans =

t=0.7,y=1.5044519 ans =

t=0.8,y=1.5608688 ans =

t=0.9,y=1.6148989 ans =

t=1.0,y=1.6668064 ②

function 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.0000000 ans =

t=0.2,y=1.9801980 ans =

t=0.3,y=1.9421173 ans =

t=0.4,y=1.8886645 ans =

t=0.5,y=1.8235382 ans =

t=0.6,y=1.7505966 ans =

t=0.7,y=1.6733644 ans =

t=0.8,y=1.5947500 ans =

t=0.9,y=1.5169573 ans =

t=1.0,y=1.4415285 b. 四阶龙格-库塔法

fun=inline('0.9*y/(1+2*t)','t','y'); [t,y]=ode45(fun,[0,1],1); [t,y] ans =

0 1.0000

0.0250 1.0222 0.0500 1.0438

0.0750 0.1000 0.1250 1.0649 1.0855 1.1056

0.1500 1.1253 0.6000 1.4259 0.1750 1.1446 0.6250 1.4404 0.2000 1.1635 0.6500 1.4547 0.2250 1.1820 0.6750 1.4689 0.2500 1.2002 0.7000 1.4828 0.2750 1.2180 0.7250 1.4967 0.3000 1.2355 0.7500 1.5103 0.3250 1.2528 0.7750 1.5239 0.3500 1.2697 0.8000 1.5372 0.3750 1.2864 0.8250 1.5505 0.4000 1.3028 0.8500 1.5636 0.4250 1.3189 0.8750 1.5765 0.4500 1.3349 0.9000 1.5894 0.4750 1.3506 0.9250 1.6021 0.5000 1.3660 0.9500 1.6147 0.5250 1.3813 0.9750 1.6271 0.5500 1.3964 1.0000 1.6395 0.5750 1.4112

plot(t,y,'o-') %画图

（2）

odefun=inline('-x*y/(1+x*x)','x','y'); [x,y]=ode45(odefun,[0,1],1); [x,y]

ans =

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.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.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 0.9285 0.4250 0.9203 0.4500 0.9119 0.4750 0.9033 0.5000 0.8944

>> plot(x,y,'o-')

0.7250 0.7500 0.7750 0.8000 0.8250 0.8500 0.8750 0.9000 0.9250 0.9500 0.9750 1.0000 0.8096 0.8000 0.7904 0.7809 0.7714 0.7619 0.7526 0.7433 0.7341 0.7250 0.7160 0.7071

6．2 用求解常微分方程初值问题的方法计算积分上限函数 一、问题描述

6．2 用求解常微分方程初值问题的方法计算积分上限函数

y(x)??f(t)dt

ax的值，实际上是将上面表达式两端求导化为常微分方程形式，并用初值条件y(a)?0。

etdt的计算问题。 试用MATLAB中的指令ode23解决定积分?2?4t5二、源程序及运行结果

odefun=inline('2*pi*exp(x)/x','x','y');

ode23t(odefun,[4,5],0);

用梯形法则算得精确解

x=4:0.1:5;y= 2*pi*exp(x)./x;trapz(x,y) ans =

129.2178

实验七

7.4观察根的变化是否很显著 一、问题描述

7.4一个10次项式

P(x)?x10?a1x9?a2x8?L?a9x?a10

的系数为

[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 -8409500 12753576 -10628640 6328800]; >> x=roots(p) x =

p=[1 -55 1320 -18150 157773 -902055 21.7335 3416930 -8409500 12753576 -10628640 7.3521 + 7.8018i 6328800]; 7.3521 - 7.8018i >> x=roots(p) 3.6433 + 3.5988i x = 3.6433 - 3.5988i

4.9248 + 1.8070i 4.9248 - 1.8070i 1.1808 + 2.0164i 1.1808 - 2.0164i 0.0643

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.0127i

7.5 将三村短路问题推广为四村短路问题，即已知四个点A，B，C，D的具体位置A（0，

0，），B(0,3),C(8,1),D(10,5),求两个点H1，H2的具体位置，使AH1+BH1+H1H2+H2C+H2D为最短。

解：设所要求的点H1?(x1,y1)构造函数

H2?(x2,y2)

S?x12?y12?x12?(y1?3)2?(x2?x1)2?(y2?y1)2?(x2?8)2?(y2?1)2?(x2?10)2?(y2?5)2

要求使 S 取得最小值的 H1?(x1,y1)H2?(x2,y2)

??S??x?1??S??y?1???S??x2???S??y即必须满足方程组 ?2?2x1?2x1?2(x2?x1)?0?2y1?2(y1?3)?2y(2?y1)?0?2(x2?x1)?2(x2?8)?2x2(?10?)?2(y2?y1)?2(y2?1)?2y(5?)2?

00?3x1?x2?0?3y?y?3?0?12??3x2?x1?18?0?3y2?y1?6?0化简即得 ?

915??2721??x1,y1????,??x2,y2???,?解得

?48??48?

这个解也即使AH1+BH1+H1H2+H2C+H2D为最短的解。

实验八

8.1分别用直接法、雅可比迭代法、赛德尔迭代法求解线性方程组AX=b。

?41?14?A??O?1???（1）

1O41?2???1?????b??O?O????1??1??4??2??10?1 ?10?10，

??5?21??25?21?1?25?21A???OOOOO??1?25?2?1?25（2）

??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(a\\b) x =

0.4793 0.0829 0.1891 0.1608 0.1891 0.0829 0.4793

雅克比：

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]; 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; end

x',pause,eorr

??????1???0????1?b??0???2??M???5???0?20?20，?0??20?1 0.1678 0.1678 0.1608 end ans =

0.5000 0.1250 0.2188 0.1953 0.2012 0.1997 0.2001 0.2000 0.2000 0.4500

赛德尔：

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); end

x',pause,eorr end ans =

0.4688 0.0781 0.1816 0.1543 0.1615 0.1596 0.0975 0.4756 （2） 直接法：

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(a\\b) x =

Columns 1 through 12

0.2421 0.0944 -0.0218 -0.0314 -0.0069 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 雅克比：

e=ones(20,1);

a=spdiags([e,-2*e,5*e,-2*e,e],[-2,-1,0,1,2],20,20);

0.1601 0.0049 0.0000 0.1600 0.0036 -0.0000 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)=0;

y(i)=(b(i)-a(i,:)*x)/a(i,i); eorr=max(abs(s-x(i)),eorr); x=y; end

x',pause,eorr end

ans =

Columns 1 through 12

0.2000 0.0800 -0.0080 -0.0192 -0.0061 0.0004 -0.0002 -0.0002 -0.0000 0.0000

Columns 13 through 20

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]; for k=1:8 eorr=0; for i=1:20

s=x(i);x(i)=0;

x(i)=(b(i)-a(i,:)*x)/a(i,i); eorr=max(abs(s-x(i)),eorr); end

x',pause,eorr end

ans =

Columns 1 through 12

0.0014 0.0000 0.0018 0.0000

0.2000 0.0800 -0.0080 -0.0192 -0.0061 0.0014 0.0018 0.0004 -0.0002 -0.0002 -0.0000 0.0000

Columns 13 through 20

0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000

8.4小行星轨道问题 一、问题描述

8.4小行星轨道问题：一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在五个不同的对小行星作了五次观察，测得轨道上五个点的坐标数据（单位：万公里）如下表所示： X坐标 P1 53605 P2 58460 P3 62859 P4 66662 P5 68894 68894 Y坐标 6026 11179 16954 23492 由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为： a1x2+2a2xy+a3y2+2a4x+2a5y+1=0 现需要建立椭圆的方程以供研究。

（1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组 AX=R 以及方程组的系数矩阵和右端项b；

（2）用MARLAB求低阶方程的指令A＼b求出待定系数a1,a2,a3,a4,a5. 二、源程序及运行结果

a=[28.7350 6.4605 0.3631 10.7210 1.2052;34.1757 13.0705 1.2497 11.6920 2.2358;

39.5125 21.3142 2.8744 12.5718 3.3908;44.4382 31.3205 5.5187 13.3324 4.6984;47.4638 94.9277 47.4638 13.7788 13.7788] a =

28.7350 6.4605 0.3631 10.7210 1.2052 34.1757 13.0705 1.2497 11.6920 2.2358 39.5125 21.3142 2.8744 12.5718 3.3908 44.4382 31.3205 5.5187 13.3324 4.6984 47.4638 94.9277 47.4638 13.7788 13.7788 >> b=[-1;-1;-1;-1;-1]’ b = -1 -1 -1 -1 -1 >> a\\b

ans = 0.0008 -0.0087 -0.0235 -0.1090

0.1747

所以X=10^4*(a\\b)= 8 -87 -235 -1090 1747

a1=8,a2=-87,a3=-235,a4=-1090,a5=1747.

本文来源：https://www.bwwdw.com/article/00ro.html