四应用篇.docx
《四应用篇.docx》由会员分享，可在线阅读，更多相关《四应用篇.docx（26页珍藏版）》请在冰豆网上搜索。
四应用篇
四．应用篇
MATLAB为我们提供了高效的编程环境，在解决实际问题和求解数学模型时，不必象使用一般的计算机语言那样将算法的每一个细节都考虑进去。
用一系列的MAATLAB的命令解决问题时，实际上是将已有的程序模块进行组合拼装。
在这样的高层次上考虑算法设计时，自然就节约了工作时间并避免了低层次的重复劳动。
73.罗伦兹吸引子的空间曲线是如何绘制的？
罗伦兹微分方程组的解曲线——Lorenz吸引子是三维空间中的一条曲线，如图4-1所示这条曲线相互缠绕而互不相交。
如果将这条曲线视为某一动点的轨迹，这个动点将随自变量t的增大，在空间中的两个定点附近作环绕运动。
图4-1罗伦兹吸引子
罗伦兹常微分方程组为
给定初值条件
取
2.2204×10－16，则得微分方程组
将三个方程的右端函数写成向量形式，得
由于MATLAB中有常数eps=2.2204×10－16，初始条件可以用列向量[00eps]T表示。
首先建立描述微分方程组右端函数的函数文件（文件名为：
flo.m）：
functionz=flo（t,y）
A=[-8./30y
（2）；0-10.10.；-y
（2）28.-1]；
z=A*y；
将这一文件保存在MATLAB的工作目录下，然后在MATLAB环境中键入如下指令：
[t,y]=ode23（'flo',0,80,[00eps]'）;
u=y（:
1）;v=y（:
2）;w=y（:
3）;plot3（u,v,w）
MATLAB的图形窗口将显示Lorenz吸引子的图形如图4-1所示。
另外，还可以绘制动画，如果关闭MATLAB的图形窗口，并在命令窗口中再键入
comet3（u,v,w）
便可以观察到Lorenz吸引子的图形生成过程的动态演示。
命令comet3的使用格式和上面的plot3的使用格式相同，不同的是绘图效果增加了动感，其功能是在绘图时以动点模拟慧星运行并带有一条尾巴。
74．如何为足球队排名次？
已知我国8支足球队在1988～1989年全国足球甲级队联赛中的成绩，按国际足联公布的鼓励进攻的三分制计算比赛积分（胜—得3分；平—得1分；负—得0分）有下面的得分表
表4-1
T1
T2
T3
T4
T5
T6
T7
T8
T1
×
4
4
9
0
3
6
2
T2
4
×
3
5
2
2
4
1
T3
4
6
×
5
3
3
3
3
T4
0
2
2
×
0
3
1
1
T5
6
2
3
6
×
7
9
7
T6
3
2
3
3
1
×
3
4
T7
0
1
3
3
0
6
×
7
T8
2
4
3
4
1
4
1
×
表中第k行数据表示第k队与其它各队比赛所得分数。
为了依据现有的数据，对各足球队的实力给予科学客观的评价，可以用不同的方法排名次。
第一种方案是依照各队的总分排名，下面的程序就是根据上表中的数据建立得分矩阵，将每一行求和得各队总分，然
后根据总分排名次。
a=[04490362;40352241;46053333;02200311;
62360797;32331034;01330607;24341410];
f=sum（a'）
[ffl]=sort（-f）
运行程序后，屏幕显示
ff=-40-28-27-21-20-19-19-9
l=51327684
由于MATLAB中的排序命令sort将数组由小到大排序，而我们需要将各队总分数按递减排序，所以将所有分数取负号后再使用排序命令。
上面结果列表如下
表4-2总分排序数据结果
名次
1
2
3
4
5
6
7
8
队号
5
1
3
2
7
6
8
4
总分
40
28
27
21
20
19
19
9
第二种方案是根据得分矩阵的正特征值所对应的特征向量排名。
利用上面程序中已经建立的得分矩阵求矩阵的特征值，用命令eig（a）可得
ans=
20.3443
-2.8104+5.1609i
-2.8104-5.1609i
-1.9756+2.1810i
-1.9756-2.1810i
-3.5585+0.8359i
-3.5585-0.8359i
-3.6552
观察数据结果知，第一个特征值为正，取这一正特征值对应的特征向量用命令
[x,d]=eig（a）;g=x（:
1）'
得
-0.3636-0.3224-0.4133-0.1438-0.5680-0.2938-0.2837-0.2880
这一向量每一元素均为负数，由特征向量的性质将其反号仍为特征向量。
用排序命令
[f，l]=sort（g）
将特征向量元素排序，得
f=-0.5680-0.4133-0.3636-0.3224-0.2938-0.2880-0.2837-0.1438
l=53126874
由排序结果列表，得
表4-3特征向量方法排序结果
名次
1
2
3
4
5
6
7
8
队号
5
3
1
2
6
8
7
4
元素
0.5680
0.4133
0.3636
0.3224
0.2938
0.2880
0.2837
0.1438
比较两种排名方案结果，尽管结果不完全一致，但第一名和最后一名是一致的。
这说明，分数是反映各队实力的一个方面，但并不是全面。
最强的队与各队比赛总会得分，而弱队与其它队比赛总不会得分。
在比赛中，一个实力中等的球队遇上强队时不会得分，遇上弱队时会得分，所以总分排名次的方法不是很公平的，特征向量法则比较公平。
当特征向量中元素均为正时，各元素的大小就反映了各球队实力的强弱。
75.怎样研究概率论中的生日问题？
假设每个人的生日在一年365天中的任意一天都是等可能性的，那么随机找n个人（不超过365人）。
这n个人生日各不相同的概率（可能性）可用公式
计算出。
因而，n个人中至少有两人生日相同这一随机事件发生的概率为
首先利用上面的公式，用计算机分别计算出一个团体的人数分别为n=1,2,…,100时的100个概率值。
用下面程序段
m=365;n=m;
fork=1:
100
p（k）=1-n/m;
m=m*365;n=n*（365-k）;
end
使用指令p
（1），p
（2），p（3），…，p（100）可得出其相应的概率值。
用指令plot（p）绘制出图形（图4-2）,可以直观地了解概率值随团体人数变化的规律。
图4-2概率曲线
现在考虑几个特殊概率值，一个有三十个同学的学生班上，至少有两个同学“生日相同”的概率为p（30）=0.7063。
五十个人的团体中，至少有两人生日在同一天（只考虑月日，不考虑年的不同）的概率为p（50）=0.9704。
在七十个人的团体中，p（70）=0.9992。
考虑团体总人数对概率值的影响，为了求在某团体中保证“至少有两人生日相同”的概率大于99%时的总人数。
键入下列命令
find（p>0.99）（回车）
计算机屏幕将显示
ans=5758596061……979899100
所以，团体人数若超过57人，则这个团体中至少有两人生日相同的概率将大于99%。
最后，考虑用计算机数值模拟。
随机产生30个正整数，介于1到365之间（代表30个同学的生日），然后统计数据，观察是否有两人以上的人生日相同。
当30个人中有两人生日相同时，计算机输出为“1”，否则输出为“0”。
如此重复观察200次，计算出这一事件发生的频率f200。
下面由MATLAB程序做计算机模拟200次重复观察：
n=0;
form=1:
200%做200次随机试验
y=0;x=1+fix（365*rand（1,30））;%产生30个随机数
fori=1:
29
forj=i+1:
30
ifx（i）==x（j）,y=1;break,end%寻找30个随机数中是否有相同数
end
end
n=n+y;%累计有两人人生日相同的试验次数
end
f=n/m%计算频率
将程序重复运行三次后，数据结果为
f200=0.7350、0.7150、0.7250
比较前面所计算的p（30）=0.7063，说明概率值和频率值比较接近，但是存在着差异。
76．怎样用动态规划算法求最短路径问题？
动态规划是用于分析一类多阶段决策过程的最优化方法。
它解决问题的思路是将一个比较复杂的问题分解为多个同一类型的子问题，然后将这些子问题按整体最优的思想逐个地求最优解，最后再求出整个问题的最优解。
一个简单的最短路径问题描述如下：
图4-3中20个结点表示20台不同地理位置上的计算机，为了将数据文件从编号为20的结点通过中间结点传输到编号为1的结点。
求最短路径（数据在各相邻结点间传送的时间如图中所示）。
被传送的数据文件从编号为20的结点出发，最
终到达编号为1的结点。
无论哪一种方案都需要水平方向走4步，垂直方向走3步。
从第20号结点到1号结点的所有路径数为
用动态规划算法解决问题的理论依据是最优性原理，最优性原理可简述为：
不论前面的状态和策略如何，以后的最优策略只取决于最初策略所确定的当前状态。
具体求解问题时要用到递归方程，递归方程所起的作用是将一个子问题的解与下一个子问题相联系。
由最优性原理可知，在从编号为20的结点到编号为1的结点的最短路径上，各结点到编号为1的结点的的路径也是各结点到编号为1的结点的最短路径。
我们称编号为1的结点为目标结点，用dk表示编号为k的结点到目标结点的最短路径长度。
由最优性原理可得
d20=min{2+d18，1+d19}，d18=min{3+d15，2+d16}，d19=min{2+d16，3+d17}，
d15=min{2+d11，2+d12}，d16=min{4+d12，2+d13}，d17=min{1+d13，2+d14}，
d11=2+d7，d12=min{1+d7，3+d8}，d13=min{1+d8，1+d9}，d14=min{3+d9，3+d10}，
d7=3+d4，d8=min{2+d4，4+d5}，d9=min{2+d5，4+d6}，d10=4+d6，
d4=1+d2，d5=min{2+d2，5+d3}，d6=3+d3，d2=2+d1，d3=3+d1，d1=0。
由d2=2，d3=3，应用逆向递推方式，得
d4=1+d2=3，d5=min{2+d2，5+d3}=4，d6=3+d3=6，……
d20=min{2+d18，1+d19}=11。
最后再根据逆向递推的各式选择写出最短路径的路线。
用动态规划算法求解最短路径问题的程序如下：
x=[221;331;412;522;553;
633;734;824;845;925;
946;1046;1127;1217;1238;
1318;1319;1439;14310;15211;
15212;16412;16213;17113;17214;
18315;18216;19216;19317;20218;
20119];
x1=x;y=[110;212;313];
fork=4:
20
l=find（x1（:
1）==k）;
d=x1（l,2）;m=x1（l,3）;
mk=length（m）;
fori=1:
mk
j=m（i）;
mm（i,:
）=find（y（:
1）==j）;
end
[v,n]=min（d+y（mm,3））;
y（k,:
）=[km（n）v];
end
z
（1）=20;l=find（y（:
1）==20）;m=y（l,2）;
k=2;z（k）=m;
whilem>1
k=k+1;
l=find（y（:
1）==m）;
m=y（l,2）;z（k）=m;
end
rood=z
Length=y（20,3）
将程序运行后，得数据结果：
rood=201916138421
Length=11
结果表明，由动态规划算法求得的最短路径长度为11，图4-4中箭头标出了最优解路径上的结点编号依次为
20→19→16→13→8→4→2→1
上面程序中，使用了一个31×3的矩阵。
由于数据传输方向是由左向右，由下向上。
故可将图作为有向图，每两个相邻结点的编号以及它们之间的一段路径长度形成一个三元数，即
[父结点号路径长子结点号]
所以，图中共31条路段的数据形成了一个31×3的矩阵。
以此为依据，从第4号结点到20号结点，逐步寻找它们到1号结点的最短路长度，同时记录每一结点到1号结点的最短路径上的子结点。
形成一个三元数
[父结点子结点局部最短路径长]
最后，形成一个20×3的矩阵。
再以此为依据，按秩序找出全局最短路径上的结点编号。
77．如何用MATLAB解决蠓虫分类问题？
生物学家试图对两类蠓虫（Af与Apf）进行鉴别,依据的资料是蠓虫的触角和翅膀的长度,已经测得9只Af和6只Apf的数据,（触角长度用x表示,翅膀长度用y表示）具体数据为:
表4-4Af类触角和翅膀长度
x
1.24
1.36
1.38
1.38
1.38
1.40
1.48
1.54
1.56
y
1.27
1.74
1.64
1.82
1.90
1.70
1.82
1.82
2.08
表4-5Apf类触角和翅膀长度数据
x
1.14
1.18
1.20
1.26
1.28
1.30
y
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是某种疾病的载体,是否应该修改分类方法。
图4-5由触角和翅膀数据画的离散点图
首先画出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.241.27;1.361.74;1.381.64;1.381.82;1.381.90;
1.401.70;1.481.82;1.541.82;1.562.08;1.141.78;
1.181.96;1.201.86;1.262.00;1.282.00;1.301.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=[xyones（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）%在散点图中画分类直线
运行上面程序可求出超定方程组的最小二乘解并画出分类边界曲线。
为了由所给数据用判别函数判别三个新蠓虫的类属，即当
时，判为Af类；当
时，判为Apf类。
运行上面程序后,键入下面命令
xx=[1.241.801;1.281.841;1.402.041]；
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改为其它的负数以重新求超定方程组的最小二乘解获得分类边界直线的方程。
这样将与前面所求分类边界直线的方程不一样，当然对新给定的蠓虫的翅膀和触角长度数据来做判断其结果也是不同的。
78．如何绘制分形曲线图形？
从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形（图4-6）；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形（图4-7），这时，图形中共有17个结点。
图4-6第一次迭代
图4-7第二次迭代
这种迭代继续进行下去可以形成Koch分形曲线。
在迭代过程中，图形中的点将越来越多，而曲线最终显示细节的多少将取决于所进行的迭代次数和显示系统的分辩率。
Koch分形曲线的定义归结于算法和计算机实现。
考虑由直线段（2个点）产生第一个图形（5个点）的过程。
设P1和P2分别为原始直线段的两个端点。
现在需要在直线段的中间依次插入三个点Q1,Q2,Q3产生第一次迭代的图形。
显然，Q1位于P1点右端直线段的三分之一处，Q3位于P1点右端直线段的三分之二处；而Q2点的位置可以认为是，由Q3点绕Q1旋转60度（逆时针方向）而得到的，故可以处理为向量Q1Q3经正交变换而得到向量Q1Q2。
形成算法如下：
（1）Q1=P1+（P2－P1）/3；
（2）Q3=P1+2（P2－P1）/3；
（3）Q2=Q1+（Q3－Q1）×A’；
（4）P5=P2；P2=Q1；P3=Q2；P4=Q3．
在算法的第三步中，A为正交矩阵：
这一算法的结果将由初始数据P1和P2点的坐标，产生5个结点的坐标。
这5个结点的坐标数组，组成一个5×2矩阵。
这一矩阵的第一行为P1的坐标，第二行为P2的坐标，…，第五行为P5的坐标。
矩阵的第一列元素分别为5个结点的X坐标，第二列元素分别为5个结点的Y坐标。
首先考虑在Koch分形曲线的形成过程中结点数目的变化规律。
设第k次迭代产生结点数为nk，第k+1次迭代产生结点数为nk+1，nk和nk+1之间的递推关系式如下
由第k次迭代的nk个结点的结点坐标数组，产生第k+1次迭代的nk+1个结点的结点坐标数组的算法可参考上面两点到五点的算法进行设计。
根据算法编写下面迭代五次的程序并绘制Koch分形曲线
p=[00;100];%给出初始数据两个点的坐标
a=[cos（pi/3）-sin（pi/3）;sin（pi/3）cos（pi/3）]; %设置用于正交变化的正交矩阵
fork=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））%根据结点坐标绘图
运行上面程序后，MATLAB的图形窗口将显示图4-8中的分形图形。
图4-8Koch分形曲线
分形曲线的绘制可以有很多丰富多彩的变化。
例如将上面分形曲线中正交变换的逆时针600旋转修改为1200旋转，并将两点连线变成一个四边形，以四边形的四个顶点坐标为初始数据绘制第五次迭代产生的分形曲线图形，将大为改观。
修改上面程序，首先，将初始数据p=[00;100]修改为
p=[00;5-5;100;55;00]
并将语句a=[cos（pi/3）-sin（pi/3）;sin（pi/3）cos（pi/3）]中的旋转达角度由
/3改为
/3，即修
升级会员 