详情
本书以最新推出的MATLAB R2013a软件为基础，详细讲解了遗传算法、免疫算法、退火算法、粒子群算法、鱼群算法、蚁群算法和神经网络算法等最常用的智能算法的MATLAB实现。书中给出的每个案例都是一个使用智能算法解决问题的具体实例，所有案例均由理论讲解、案例背景、MATLAB程序实现和扩展阅读几个部分组成，并配有完整的原创程序，使读者在掌握算法的同时更能快速提高使用算法求解实际问题的能力。
工程软件应用精解
北京
图书在版编目（CIP）数据
MATLAB智能算法超级学习手册/高飞编著.--北京：人民邮电出版社，2014.5
ISBN 978-7-115-34879-1
Ⅰ.①M… Ⅱ.①高… Ⅲ.①Matlab软件—手册 Ⅳ.①TP317-62
中国版本图书馆CIP数据核字（2014）第046848号
内容提要
MATLAB为广大科研工作者的必备工具之一，智能算法在工程实际上得到较广泛的应用。本书基于MATLAB R2013a 软件，全面地介绍和举例验证智能算法的有效性。
智能算法种类较多，本书的内容主要包括马尔科夫链模型、层次分析法、粒子群算法、遗传算法、蚁群算法、鱼群算法、PID控制算法、神经网络算法等。智能算法对于很多初学者而言，有一定的困难，很难理解程序流程、数据的运算过程，因此给实际应用带来困难。本书将围绕智能算法展开综述，深入浅出地介绍和分析各类智能算法，用智能算法解决工程应用问题。
本书以工程应用为目标，深入浅出，实例引导，讲解详实，适合作为理工科高等院校研究生、本科生的教学用书，也可作为广大科研和工程技术人员的参考用书。
◆编著 MATLAB技术联盟 高飞
责任编辑 王峰松
责任印制 彭志环 杨林杰
◆人民邮电出版社出版发行 北京市丰台区成寿寺路11号
邮编 100164 电子邮件 315@ptpress.com.cn
网址 http://www.ptpress.com.cn
北京艺辉印刷有限公司印刷
◆开本：787×1092 1/16
印张：31.25
字数：738千字 2014年5月第1版
印数：1-3000册 2014年5月北京第1次印刷
定价：69.00元
读者服务热线：(010)81055410 印装质量热线：(010)81055316
反盗版热线：(010)81055315
广告经营许可证：京崇工商广字第0021号
在科学研究和工程计算领域经常会遇到一些非常复杂的计算问题。这些问题利用计算器或手工计算无法完成，只能借助计算机完成。MATLAB在数值计算方面表现卓越，同时，MATLAB语言具有编程效率高、图形界面友好、扩充能力强、交互性好、可移植性强、全方位的帮助系统等特点。因此，MATLAB广泛应用于各行各业。
目前，MATLAB已成为信号处理、通信原理、自动控制等专业的重要基础课程的首选实验平台。对于学生而言，最有效的学习途径是结合某一专业课程的学习掌握该软件的使用与编程。
1．本书特点
由浅入深，循序渐进：本书以MATLAB爱好者为对象，首先从MATLAB使用基础讲起，再辅以MATLAB智能算法在工程中的应用案例，帮助读者尽快掌握MATLAB进行智能算法的学习和开发。
步骤详尽、内容新颖：作者结合多年 MATLAB 使用经验与实际工程应用案例，将MATLAB 软件的使用方法与技巧详细地讲解给读者。本书在智能算法的讲解过程中，步骤详尽，并辅以相应的流程图，使读者在阅读时一目了然，从而快速掌握书中所讲内容。
实例典型，轻松易学：通过学习实际工程应用案例，运用智能算法进行操作求解，是掌握MATLAB智能算法编程应用最好的方式。本书通过综合应用案例，透彻详尽地讲解了MATLAB智能算法的应用研究。
2．本书内容
本书基于MATLAB R2013a 版，讲解了MATLAB 的基础知识和核心内容。本书主要围绕智能算法在工程问题中的应用进行算法验证。全书包括23章。
第1章，MATLAB基础知识：对MATLAB的基本界面和功能进行介绍，并围绕矩阵的应用以及简单的工程问题求解进行阐述。
第2章，种群竞争微分方程的求解：以种群竞争模型为例，利用MATLAB对该微分模型进行分析求解。
第3章，基于Markov的食品物价趋势预测：主要围绕Markov模型进行食品价格的预测分析。
第4章，基于时间序列的物价预测算法：围绕常用的时间序列模型，进行食品价格数据的分析。
第5章，基于层次分析法的食堂服务质量评价算法：基于食堂服务质量模型，运用层次分析方法，给出合理全面的评价。
第6章，MATLAB优化工具箱的使用：MATLAB优化工具箱有很多。针对极值寻优计算，MATLAB 优化工具箱提供了大量的函数以供调用。本章讲述了优化工具箱函数的使用方法。
第7章，基于RBF 网络的优化逼近：围绕径向基函数（RBF，Radial Basis Function）神经网络，讲述了使用智能算法优化的RBF网络的优化逼近算法分析。
第8章，自适应模糊控制算法：自适应模糊控制是指具有自适应学习算法的模糊逻辑系统，其学习算法是依靠数据信息来调整模糊逻辑系统的参数。本章围绕常规的二阶控制系统，进行自适应模糊控制算法的研究分析。
第9章，基于PID的控制算法：PID是现今工业控制广泛采用的控制算法。本章围绕PID控制算法，研究了不同的PID算法的控制效果。
第10章，基于LQR+PID的倒立摆控制算法：针对倒立摆系统，分别采用PID和LQR控制算法对倒立摆系统进行控制，达到稳定控制的目的。
第11章，基于粒子群算法的寻优计算：粒子群算法（PSO）是一种基于群体的随机优化技术，根据迭代寻优算法找出目标函数的极值。本章主要使用粒子群算法对有约束和无约束的函数方程进行粒子群算法验证。
第12章，基本粒子群改进算法分析：粒子群算法寻优迭代过程中，常常出现粒子早熟等现象，陷入局部最优解。因此，本章基于基本的粒子群算法讨论了粒子群算法的改进策略。
第13章，基于免疫算法的物流中心选址：物流中心选址模型是非凸和非光滑的带有复杂约束的非线性规划模型，属于 NP-hard 问题。本章采用免疫算法，一方面对算法进行验证，另一方面拟解决物流中心选址等问题。
第14章，基于人工免疫粒子群聚类算法：聚类算法多根据数据特征进行相关性近亲聚类，基于人工免疫粒子群算法能够较快、较准确地实现图像数据的聚类分析。
第15章，基于ART的植物种类自动分类：目前植物识别和分类主要由人工完成，工作量较大，效率较低。因此，本章采用ART神经网络对植物种类数据进行自动分类分析，简化实际工程应用工作量。
第16章，基于贝叶斯网络的数据预测：贝叶斯预测不同于传统预测方法。贝叶斯预测在预测过程中应用了决策者的主观信息。本章基于贝叶斯网络，应用贝叶斯网络算法进行数据的预测分析。
第17章，基于遗传算法的寻优计算：遗传算法具有较好的收敛性和健壮性。本章采用遗传算法对函数进行寻优计算，验证算法的可行性。
第18章，基于遗传算法的TSP求解：TSP也称旅行商问题，属于NP问题。采用遗传算法能较快地实现TSP算法的求解。
第19章，基于蚁群算法的路径规划计算：蚁群具有较好的寻优能力，本章主要采用蚁群算法模拟了二维和三维路径规划优化计算问题。
第20章，基于蚁群算法的TSP求解：本章主要借助蚁群算法实现TSP问题的求解，验证算法的可行性。
第21章，基于模拟退火的粒子群算法：模拟退火算法基于对固体退火过程的模拟。基于模拟退火的粒子群优化算法以基本粒子群优化算法运算流程作为主体流程，为把模拟退火机制引入其中，采用杂交粒子群优化算法中的杂交运算和带高斯变异的粒子群优化算法中的变异运算，以进一步调整优化群体。
第22章，基于人群搜索算法的函数优化：人群搜索算法是对人的随机搜索行为进行分析。本章采用人群搜索算法实现函数的寻优计算，验证了算法的可行性。
第23章，数控机床进给伺服系统的SOA-PID参数整定：基于SOA对于函数优化特性的分析，本章主要采用人群搜索算法实现对实际数控机床进给伺服系统进行PID参数的整定。
注：本书中用到的所有程序代码和数据，请到作者的博客下载。
3．读者对象
本书适合于MATLAB初学者和研究算法提高并解决工程应用能力的读者，包括：
★相关从业人员 ★初学MATLAB的技术人员
★大中专院校的教师和在校生 ★相关培训机构的教师和学员
★参加工作实习的“菜鸟” ★MATLAB爱好者
★广大科研工作人员 ★初中级MATLAB从业人员
4．本书作者
本书由MATLAB技术联盟高飞编著，另外，余胜威、李昕、刘成柱、史洁玉、孙国强、代晶、贺碧蛟、石良臣、孔玲军、柯维娜等人为本书的编写提供了大量的帮助，在此一并表示感谢。
虽然作者在本书的编写过程中力求叙述准确、完善，但由于水平有限，书中欠妥之处在所难免，希望读者能够及时指出，共同促进本书质量的提高。
5．读者服务
读者如在学习过程中遇到与本书有关的技术问题，可以发邮件到邮箱 book_hai@126.com，或者访问博客http://blog.sina.com.cn/tecbook，编者会尽快给予解答，我们将竭诚为您服务。
编者
MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似。用MATLAB解决问题要比用C、FORTRAN等语言简捷得多，并且MATLAB吸收了Maple 等软件的优点，从而成为一个强大的数学软件。本章从最基本的运算单元出发，讲述了 MATLAB 矩阵的表示方法、符号变量的应用、线性方程组的求解，并着重讲解了MATLAB在工程上的简单应用研究。
学习目标：
（1）熟练掌握MATLAB矩阵的表示方法；
（2）熟练运用符号变量求解实际物理模型；
（3）熟练掌握线性方程组的求解方法，包括齐次线性方程组和非齐次线性方程组的求解；
（4）熟练使用MATLAB工具解决简单工程问题。
本书基于MATLAB R2013a 版本进行程序设计，涉及程序在2009版本以及以后版本均可以运行。在集成开发环境下，MATLAB集成了管理文件、变量和应用程序的许多编程工具。在MATLAB桌面上可以得到和访问的窗口主要有：
命令窗口（The Command Window）
命令历史窗口（The Command History Window）
启动平台（Launch Pad）
编辑调试窗口（The Edit/Debug Window）
工作台窗口和数组编辑器（Workspace Browser and Array Editor）
帮助空间窗口（Help Browser）
当前路径窗口（Current Directory Browser）
单击HOME页，在界面下的布局（Layout）中可选择性选择显示的窗口。例如，在图1-1 中，界面只显示了 Current Folder、Workspace、Command Window 窗口，默认打开的Command History窗口关掉。
MATLAB支持程序的开发，并且内部函数的代码也是开源的，用户可以调用自己设计的程序文件。例如，图1-2为MATLAB程序脚本文件，用户可在里面书写代码并修改和调试，很方便。直接在HOME页单击New，默认的文件名为untitled.m文件。
利用MATLAB进行2D图形的绘制。MATLAB默认的图像，曲线颜色是自动编号的，曲线颜色呈现是不同的。MATLAB R2013a 版本提供了快捷的画图界面，在Plots 界面有很多不同类型的输出图形，简单的2D曲线如图1-3所示。
对于3D图形的绘制，MATLAB也能较好、较迅速地表征出来。例如，图1-4为MATLAB Logo的3D曲面造型，采用mesh()函数可快速地表征MATLAB曲面。
对于图1-3、图1-4所示的图形，MATLAB数据都存在矩阵数组中，用whos命令可产生一个在当前工作区内的所有变量和数组状况表。在 Command Window 直接输入 whos，显示的为Workspace所有变量的属性，如图1-5所示。
MATLAB函数查询如图1-6所示，单击输入命令行左侧的 fx，上拉一个对话框，在对话框输入函数名称，即可进行相关查询。对于查询到的函数，单击鼠标，在右侧出现该函数的用法对话框，可以进行该函数的用法预览。
MATLAB提供了help查询功能。单击HOME页上的help，会弹出如图1-7所示的查询界面。该界面列出了MATLAB工具箱，用户可以进行相关函数和demo的查询。
MATLAB集成了很多工具箱，不同版本MATLAB工具箱的更新程度不同。对于查询工具箱种类以及查询工具箱版本，可直接在MATLAB Command Window窗口输入ver命令，按回车键可得到MATLAB版本信息，如图1-8所示。
MATLAB功能相当强大，几乎可以解决所有的工程分析问题，而且MATLAB计算精度较高，拥有强大的工具箱和矩阵处理能力，被广大学术界的研究人员所认可。因此，MATLAB是一款高效的科学计算软件。
矩阵和向量是一样的，用来描述某一个问题的方程组的系数、由方程组的系数和常数构成的方阵。矩阵包括数值矩阵、符号矩阵、特殊矩阵等3种基本样式。
1．实数矩阵的输入
MATLAB的强大功能之一体现在能直接处理向量或矩阵。前提是用户根据具体的问题输入待处理的向量或矩阵。
一般简单的定义矩阵，可以直接按行方式输入每个元素：同一行中的元素用逗号（，）或者用空格符来分隔，且空格个数不限；不同的行用分号（；）分隔。
所有元素处于一个方括号（[ ]）内。当矩阵是多维（三维以上）的，且方括号内的元素是维数较低的矩阵时，会有多重方括号。
【例1-1】实数矩阵输入实例。
>> T=[11 12 1 2 3 4 5 6 7 8 9 10]
T=
Columns 1 through 11
11 12 1 2 3 4 5 6 7 8 9
Column 12
10
>>X=[2.32 3.43;4.37 5.98]
X=
2.3200 3.4300
4.3700 5.9800
>> va=[1 2 3 4 5]
va=
1 2 3 4 5
>>MB=[1 2 4;2 3 3;5 4 5]
MB=
1 2 4
2 3 3
5 4 5
>>Null=[ ] %生成一个空矩阵
Null=
[]
2．复数矩阵的输入
复数在现行的控制工程以及复平面计算中应用较多。复数矩阵是指带有虚数的数值矩阵。复数矩阵的生成方式如例1-2和例1-3所示。
【例1-2】复数矩阵输入方式一实例。
>> a=1.7;b=3/25;
C=[1,3*a+i*b,b*sqrt(a); sin(pi/5),a+7*b,3.9+1]
C=
1.0000 5.1000+ 0.1200i 0.1565
0.5878 2.5400 4.9000
【例1-3】复数矩阵输入方式二实例。
>> R=[1 2 3;4 5 6],M=[11 12 13;14 15 16]
R=
1 2 3
4 5 6
M=
11 12 13
14 15 16
>> RM=R+i*M
RM=
1.0000 +11.0000i 2.0000 +12.0000i 3.0000 +13.0000i
4.0000 +14.0000i 5.0000 +15.0000i 6.0000 +16.0000i
在 MATLAB 中输入符号向量或者矩阵的方法和输入数值向量或者矩阵在形式上很相似，只不过要用到符号矩阵定义函数sym，或者是用到符号定义函数syms。先定义一些必要的符号变量，再像定义普通矩阵一样输入符号矩阵。
1．用命令sym定义矩阵
这时的函数sym实际上在定义一个符号表达式，符号矩阵中的元素可以是任何符号或者表达式，而且长度没有限制，只是将方括号置于用于创建符号表达式的单引号中。
【例1-4】用命令sym定义矩阵实例。
>> sym_m=sym('[a b c;Jack,Help Me!,NO WAY!]')
sym_m=
[ a, b, c, 0, 0]
[ Jack,Help,factorial(Me),NO,factorial(WAY)]
>> sym_d=sym('[1 2 3;a b c;sin(x) cos(y) tan(z)]')
sym_d=
[ 1, 2, 3]
[ a, b, c]
[ sin(x),cos(y),tan(z)]
2．用命令syms定义矩阵
先定义矩阵中的每一个元素为一个符号变量，然后如数值矩阵操作那样输入符号矩阵。
【例1-5】用命令syms定义矩阵实例。
>> syms a b c
>> M1=sym('Classical');
>> M2=sym('Claysw');
>> M3=sym('yellow');
>> yswM123=[a,b,c;M1,M2,M3;2,3,5;5,4,6]
yswM123=
[ a, b, c]
[Classical,Claysw,yellow]
[ 2, 3, 5]
[ 5, 4, 6]
3．把数值矩阵转化成相应的符号矩阵
数值型和符号型在 MATLAB 中是不相同的，它们之间不能直接进行转化。MATLAB提供了一个将数值型转化成符号型的命令，即sym。
【例1-6】数值型转化成符号型实例。
>> Digit_Ma=[1/3 sqrt(3) 3.1;exp(0.3) log(10) 23^.5]
Syms_Ma=sym(Digit_Ma)
Digit_Ma=
0.3333 1.7321 3.1000
1.3499 2.3026 4.7958
Syms_Ma=
[ 1/3, 3^(1/2), 31/10]
[ 3039611811401035/2251799813685248,2592480341699211/1125899906842624,23^(1/2)]
注意：思考矩阵是用分数形式还是浮点形式表示的。一般情况下，矩阵是以浮点型变量保存的。针对本例，矩阵转化成符号矩阵后以最接近原值的有理数形式或者函数形式表示。
（1）全零阵
函数 zeros
格式 B=zeros(n) %生成n×n全零阵
B=zeros(m,n) %生成m×n全零阵
B=zeros([mn]) %生成m×n全零阵
B=zeros(d1,d2,d3…) %生成d1×d2×d3×…全零阵或数组
B=zeros([d1 d2 d3…]) %生成d1×d2×d3×…全零阵或数组
B=zeros(size(A)) %生成与矩阵A大小相同的全零阵
（2）单位阵
函数 eye
格式 Y = eye(n) %生成n×n单位阵
Y = eye(m,n) %生成m×n单位阵
Y = eye(size(A)) %生成与矩阵A大小相同的单位阵
（3）全1阵
函数 ones
格式 Y=ones(n) %生成n×n全1阵
Y=ones(m,n) %生成m×n全1阵
Y=ones([m n]) %生成m×n全1阵
Y=ones(d1,d2,d3…) %生成d1×d2×d3×…全1阵或数组
Y=ones([d1 d2 d3…]) %生成d1×d2×d3×…全1阵或数组
Y=ones(size(A)) %生成与矩阵A大小相同的全1阵
（4）均匀分布随机矩阵
函数 rand
格式 Y=rand(n) %生成n×n随机矩阵，其元素在（0，1）内
Y=rand(m,n) %生成m×n随机矩阵
Y=rand([m n]) %生成m×n随机矩阵
Y=rand(m,n,p,…) %生成m×n× p×…随机矩阵或数组
Y=rand([mn p…]) %生成m×n× p×…随机矩阵或数组
Y=rand(size(A)) %生成与矩阵A大小相同的随机矩阵
rand %无变量输入时只产生一个随机数
s=rand('state') %产生包括均匀发生器当前状态的35个元素的向量
rand('state',s) %使状态重置为 s
rand('state',0) %重置发生器到初始状态
rand('state',j) %对整数 j重置发生器到第 j个状态
rand('state',sum (100*clock)) %每次重置到不同状态
（5）正态分布随机矩阵
函数 randn
格式 Y=randn(n) %生成n×n正态分布随机矩阵
Y=randn(m,n) %生成m×n正态分布随机矩阵
Y=randn([m n]) %生成m×n正态分布随机矩阵
Y=randn(m,n,p,…) %生成m×n× p×…正态分布随机矩阵或数组
Y=randn([mn p…]) %生成m×n× p×…正态分布随机矩阵或数组
Y=randn(size(A)) %生成与矩阵A大小相同的正态分布随机矩阵
randn %无变量输入时只产生一个正态分布随机数
s=randn('state') %产生包括正态发生器当前状态的2个元素的向量
s=randn('state',s) %重置状态为 s
s=randn('state',0) %重置发生器为初始状态
s=randn('state',j) %对于整数 j重置状态到第 j个状态
s=randn('state',sum(100*clock)) %每次重置到不同状态
（6）产生随机排列
函数 randperm
格式 p=randperm(n) %产生1～n之间整数的随机排列
（7）产生线性等分向量函数linspace
格式 y=linspace(a,b) %在 (a,b)上产生100个线性等分点
y=linspace(a,b,n) %在 (a,b)上产生n个线性等分点
（8）产生对数等分向量
函数 logspace
格式 y=logspace(a,b) %在 (10a,10b)之间产生50个对数等分向量
y=logspace(a,b,n) %在 (10a,10b)上产生n个对数等分向量
y=logspace(a,pi) %在 (10a,π)上产生50个对数等分向量
（9）计算矩阵中元素个数
函数 numel
格式 n=numel(a) %返回矩阵A中元素的个数
（10）产生以输入元素为对角线元素的矩阵
函数 blkdiag
格式 out=blkdiag(a,b,c,d,…) %产生以a,b,c,d,…为对角线元素的矩阵
（11）友矩阵
函数 compan
格式 A=compan(u) %u为多项式系统向量，A为友矩阵，A的第1行元素为−u(2 :n) /u(1)，其中u(2:n)为u的第2到n个元素，A的特征值就是多项式的特征根
（12）Hadamard矩阵
函数 hadamard
格式 H=hadamard(n) %返回n 阶Hadamard矩阵
（13）hankel方阵
函数 hankel
格式 H=hankel(c) %第1列元素为c，反三角以下元素为0
一个元素不同，交叉位置的元素取为c的最后一个元素
H=hankel(c,r) %第 1列元素为c，最后一行元素为 r。如果c的最后一个元素与 r的第（14）Hilbert矩阵
函数 hilb
格式 H=hilb(n) %返回n阶Hilbert矩阵，其元素为H (i,j)=1/ (i+ j−1)。
（15）逆Hilbert矩阵
函数 invhilb
格式 H=invhilb(n) %产生n阶逆Hilbert矩阵
（16）Magic（魔方）矩阵
函数 magic
格式 M=magic(n) %产生n阶魔方矩阵
（17）Pascal矩阵
函数 pascal
格式 A=pascal(n) %产生n阶Pascal矩阵。它是对称正定矩阵，元素由Pascal三角组成。它的逆矩阵的所有元素都是整数
A=pascal(n,1) %返回由下三角的Cholesky系数组成的Pascal矩阵
A=pascal(n,2) %返回pascal(n,1)的转置和交换的形式
（18）托普利兹矩阵
函数 toeplitz
格式 T=toeplitz(c,r) %生成一个非对称的托普利兹矩阵，将c作为第1列，将 r作为第1行，其余元素与左上角相邻元素相等
T=toeplitz(r) %用向量 r生成一个对称的托普利兹矩阵
（19）Wilkinson特征值测试阵
函数 wilkinson
格式 W=wilkinson(n) %返回n阶Wilkinson特征值测试阵
符号变量在解决工程问题中的应用较多。对于一个工程问题而言，一般首先从变量出发，把问题用符号变量表示出来（得到符号矩阵），然后通过符号变量求解得到一般表达式，再根据该表达式，代入相应的初始条件，即可得到问题的具体的解。
本节主要从符号变量与实际生活实证分析出发，应用符号变量求解质点系转动惯量、油罐剩余油量体积和光的反射定理等问题。
已知在平面上的n个质点P1( x1，y1)，P2( x2，y2)，…，Pn( xn，yn)，其质量分别为m1，m2，…，mn。请确定一个点P(x,y)，使得质点系关于此点的转动惯量最小。
解：设质点系关于此点的转动惯量为J，由转动惯量定义可知：
要满足质点系的转动惯量最小，即
最小，这是一个二元一次极值问题。
由式（1.1）可知：
式（1.2）满足最小值条件时，
，可得：
由式（1.3）可得：
则P点在
处有极小值，此时质点系关于P的转动惯量最小。
油罐在一般的加油站均有应用。如何求油罐剩余油量的体积是一个亟待解决的问题。例如，如图1-9（a）所示，一平放的椭圆柱体形状的油罐，长度为L，椭圆的长半轴为a，短半轴为b，油的密度为ρ，问：当油罐中油的高度为h时，油量是多少？
解：由题意可知，该柱体在长度方向上是均匀的，故在此取该椭圆柱体的一横截面进行分析探讨，如图1-9（b）所示。
假设椭圆柱体的横截面为标准的椭圆形，且椭圆柱体完好无损，放置平稳，外界干扰可忽略不计；椭圆柱体里油高为h的油面与所建坐标系上的椭圆柱体横截面相交，且一交点为A点，A点坐标为(x,y)。
图1-9（b）所示的阴影部分面积为椭圆柱体中的油所占空间，积分可得：
其中，S为椭圆柱体的横截面面积。对式（1.4）积分可得：
则可求出：
从而可求出体积V：
椭圆柱体油罐中油的高度为h，与模型中的y存在如下关系：
将式（1.8）代入式（1.7）可得：
其中，V为椭圆柱体里油高为h的油的体积。
则高度为h的油量为：
其中m为椭圆柱体里油高为h的油的质量（油量）。由式（1.10）得：
该模型结果符合题目要求。程序如下：
>> syms a b h y L
m=sqrt(b^2-y^2);
m1=int(m);
m2=int(m,'-b','n')
m2=
(b^2*asin(b/(b^2)^(1/2)))/2 + (b^2*asin(n/(b^2)^(1/2)))/2 + (n*(b^2 - n^2)^(1/2))/2
>> m3=subs(m2,'n','y');
>> S=2*a/b*m3;
>> simplify(S)
ans=
a*b*(asin(b/(b^2)^(1/2)) + asin(y/(b^2)^(1/2))) + (a*y*(b^2 - y^2)^(1/2))/b
>> V=S*L
V=
(2*L*a*((b^2*asin(b/(b^2)^(1/2)))/2 + (b^2*asin(y/(b^2)^(1/2)))/2 + (y*(b^2 - y^2)^(1/2))/2))/b
>> V=2*a/b*(1/2*y*(b^2-y^2)^(1/2)+1/2*b^2*atan(y/(b^2-y^2)^(1/2)))*L;
>> y=h-b;
>> V1=subs(V,'y','h-b')
V1=
-(2*L*a*((b^2*atan((b - h)/(b^2 - (b - h)^2)^(1/2)))/2 + ((b^2 - (b - h)^2)^(1/2)*(b - h))/2))/b
>> simplify(V1)
ans=
- L*a*b*atan((b - h)/(b^2 - (b - h)^2)^(1/2)) - (L*a*(b^2 - (b - h)^2)^(1/2)*(b - h))/b
>>
光的发射定理最早由费马提出（费马原理）光总是沿用时最短的光程传播。试根据这一原理利用极值的有关知识证明光的反射定律：入射角等于反射角。下文将借助于符号变量证明入射角等于反射角。
根据题意，光线的入射、反射过程可由图1-10直观地表示出来。在图1-10中，光线从1入射，反射到2点。
针对图1-10 所示的光线反射路径图，假设一束自然光线沿路径L10照射到x轴，y轴设为实物体表面，且为理想状态，光线传播过程中无阻碍，与法线y轴的夹角为θ1；光线经实物体表面x轴反射后，沿路径L02反射，与法线y轴的夹角为θ2。由费马原理可得，路径L10、路径L02为直线；光线从1点到2点在坐标轴上的竖直方向上的投影相等，且为H；光在空气中传播的速度为光线在真空中传播的速度C；1点与2点之间的距离为定值L；光线从L10到L02所需时间为T。
则光线从L10到L02所需时间为：
其中，θ1、θ2∈(0°，90°)。
又1点与2点之间的距离为定值L，则可得：
设，则
由三角代换变形得：
代入式（1.12）得：
对式（1.15）求一阶导数得：
因函数时间T有极小值，令T'=0，整理结果可得：
此时，将式（1.17）代入式（1.14）得：
从而有：
由式（1.19）可得入射角等于反射角，即结果成立。故由费马原理（光总是沿用时最短的光程传播）可知，实际光线从1点到达2点，经过路径L10、L02，与法线y轴所成夹角均为θ。
程序代码如下。
>> syms H C K x
>> T=(H/C)*((1/cos(x))+[1+(K-tan(x))^2]^(1/2));
>> dfdx=diff(T,x)
dfdx=
H/C*(1/cos(x)^2*sin(x)+1/(1+(K-tan(x))^2)^(1/2)*(K-tan(x))*(-1-tan(x)^2))
>> a=solve(dfdx,'x');
>> tan(a)
ans=
1/2*K
1/2*K
线性方程组的求解在日常生活中的应用较多，特别是解决企业规划、任务分配等问题。线性方程组的求解一般分为两类：一类是求唯一解或求特解，另一类是求通解。可以通过由MATLAB求解线性方程组系数矩阵的秩来判断：
若系数矩阵的秩r=n（n为方程组中未知变量的个数），则有唯一解；
若系数矩阵的秩r<n，则可能有无穷解；
线性方程组的通解（无穷解）=对应齐次方程组的通解 + 非齐次方程组的一个特解，其特解的求法属于解的第一类问题，通解部分属第二类问题。
在MATLAB中，函数null( )用来求解零空间，即满足A·X=0的解空间，实际上是求出解空间的一组基。
格式 z=null % z的列向量为方程组的正交规范基，满足Z′×Z=I
z=null(A,’r’) % z的列向量是方程A·X=0的有理基
【例1-7】求解方程组
的通解。
解：MATLAB求解程序代码如下。
>>A=[1 2 2 1;
2 1 -2 -2;
1 -1 -4 -3]; %原始系数矩阵
format rat %指定有理式格式
B=null(A,'r') %求解空间的有理基
B=
2 5/3
-2 -4/3
1 0
0 1
或通过最简行得到基：
>> B=rref(A)
B=
1 0 -2 -5/3
0 1 2 4/3
0 0 0 0
则相应地写出线性方程组的通解：
syms k1 k2 %定义符号变量
X=k1*B(:,1)+k2*B(:,2) %写出方程组的通解
%运行结果显示
X=
2*k1 + (5*k2)/3
- 2*k1 - (4*k2)/3
k1
k2
pretty(X) %让通解表达式更精美
+- -+
| 5 k2 |
| 2 k1+ ---- |
| 3 |
| |
| 4 k2 |
| - 2 k1 - ---- |
| 3 |
| |
| k1 |
| |
| k2 |
+- -+
需要先判断非齐次线性方程组是否有解，若有解，然后求通解，步骤如下。
Step1：判断A·X=b是否有解，若有解，则进行第二步，否则终止求解；
Step2：求A·X=b的一个特解；
Step3：求A·X=0的通解；
Step4：A·X=b的通解等于 A·X=b的通解加上A·X=b的一个特解。
【例1-8】求解方程组 
解：在MATLAB中建立脚本M文件：
A=[1 -2 3 -1;3 -1 5 -3;2 1 2 -2];
b=[1 2 3]';
B=[A b];
n=4;
RA=rank(A)
RB=rank(B)
format rat
if RA==RB&RA==n %判断是否有唯一解
X=A\b
elseif RA==RB&RA<n %判断是否有无穷解
X=A\b %求特解
C=null(A,'r') %求AX=0的基础解系
else X='equition no solve' %判断无解
end
运行后结果显示为：
RA=
2
RB=
3
X=
equition no solve
【例1-9】求解方程组： 
解法一：在MATLAB编辑器中建立M文件：
A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
n=4;
R_A=rank(A)
R_B=rank(B)
format rat
if R_A==R_B&R_A==n
X=A\b
elseif R_A==R_B&R_A<n
X=A\b
C=null(A,'r')
else X='Equation has no solves'
end
运行后结果显示为：
R_A=
2
R_B=
2
Warning: Rank deficient,rank=2,tol=3.826647e-15.
X=
0
0
-8/15
3/5
C=
3/2 -3/4
3/2 7/4
1 0
0 1
所以，原方程组的通解为 
解法二：用rref( )求解：
>> A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8];
b=[1 4 0]';
B=[A b];
C=rref(B) %求增广矩阵的行最简
运行后结果显示为：
C=
Columns 1 through 5
1 0 -3/2 3/4 5/4
0 1 -3/2 -7/4 -1/4
0 0 0 0 0
对应齐次方程组的基础解系为：
；非齐次方程组的特解为： 
所以，原方程组的通解为X=k1×ξ1+k2×ξ2+η*。
函数symmlq的格式如下：
x=symmlq(A,b) %求线性方程组A·X=b的解X。A必须为n阶对称方阵，b为n元列向量。A可以是由afun定义并返回A × X 的函数。如果收敛，将显示结果信息；如果收敛失败，将给出警告信息并显示相对残差norm(b−A·X)/norm(b)和计算终止的迭代次数
symmlq(A,b,tol) %指定误差 tol，默认值是1e−6
symmlq(A,b,tol,maxit) %maxit指定最大迭代次数
symmlq(A,b,tol,maxit,M) %M为用于对称正定矩阵的预处理因子
symmlq(A,b,tol,maxit,M1,M2) %M=M1×M2
symmlq(A,b,tol,maxit,M1,M2,x0) % x0为初始估计值，默认值为0
[x,flag]=symmlq(A,b,…) % flag的取值为：0表示在指定迭代次数内按要求精度收敛；1表示在指定迭代次数内不收敛；2表示M 为坏条件的预处理因子；3表示两次连续迭代完全相同；4表示标量参数太小或太大；5表示预处理因子不是对称正定的
[x,flag,relres]=symmlq(A,b,…) % relres表示相对误差norm(b−A·x) /norm(b)
[x,flag,relres,iter]=symmlq(A,b,…) % iter表示计算 x的迭代次数
[[x,flag,relres,iter,resvec]=symmlq(A,b,…) % resvec表示每次迭代的残差：norm(b−A·x0)
[x,flag,relres,iter,resvec,resveccg]=symmlq(A,b,…) %resveccg 表示每次迭代共轭梯度残差的范数
MATLAB在工程上的应用较多，例如机械机构优化分析、机械控制、通信领域、数值计算等方面。MATLAB因其强大的数据处理能力，逐渐成为工程应用领域主导辅助工具。本节主要应用MATLAB解决简单的工程问题，例如机械中的内燃机转角与升程之间的插值拟合关系、海上航行区域警示线问题、物种竞争模型的求解等。这类问题的优化求解对解决复杂的工程问题起着关键的作用。
内燃机的气门开启状态是由与发动机转轴连接的凸轮来控制的，气门开启的大小决定了发动机的性能。因此，驱动气门的凸轮加工曲线是非常重要的。表1-1 给出了某凸轮转角与升程的几组实测数据。为了机床加工的需要，试确定ϕ从91°到130°每隔1°的升程。
解：这是个一元函数插值问题。运用 MATLAB 软件，采用双三次插值法“cubic”，易求得ϕ从91°到130°每隔1°的升程，如表1-2所示。
根据表1-2，ϕ从91°到130°每隔1°的升程图像如图1-11所示。
从图1-11 中可看出，插入值与原函数值组成的图像能很好地吻合，故双三次插值法“cubic”较精确。
程序如下。
clc，clear，close all
x=[91，105，110，115，120，124，128];
h=[0，1.0869，1.9710，2.7555，3.3986，4.9073，8.3409];
j=1;
for i=91:130
y(1，j)=interp1(x，h，i，'cubic');
j=j+1;
end
x1=91:130;
plot(x，h，'r')
hold on
plot(x1，y，'o')
C=contour(x2，y2，z2); %绘制等位线
海上运输涉及世界上各个水域。由于船舶航行的航线、区域、季节的不同，其海上风险也不一样。因此，为了避免特殊风险，船舶不能驶入某些区域或一定时间内的某些地区，相关部门也会根据一定时间内该区域的水位情况予以警示。本小节针对航行区域警示线模型进行分析求解。
某海域上各种吨位的船只频繁地经过。为了保证船只的航行安全，有关机构在低潮时对水深进行了测量。表1-3给出了测量数据。
表1-3中，(x，y)为测量点，z为(x，y)处的水深（单位：米）。船的吨位可以用其吃水深度来反映，分为4吨、4.5吨、5吨和5.5吨4档。
现在，航运部门要在矩形海域（75，200）×（−50，150）上为不同吨位的航船设置警示标志。根据测量数据描述该海域的地貌，并绘制不同吨位的警示线，供航运部门使用。相邻监测点之间地势没有剧烈变化；测量时该海域潮高没有明显变化；忽略大气压引起的水位变化，以及测量仪器引起的测量误差；题中数据基本符合实际情况。
航行区域的警示线关系船只航行安全，即由于不同船只载重的不同，吃水深度也不同，因此行船需要选择合适的航道以满足其航行需要。由题意，把水深z作为x，y的函数，需要根据已知数据运用插值函数的求解算法，通过求解，作地貌和等深线图，由图能清楚地了解该海域的复杂地貌。
解：由题意知，船吃水深度划分为4档。为了使过往船只能够清楚该海域的水深情况，需要根据等深线的不同，标注不同吃水深度船只的警示线。使用MATLAB操作如下。
Step1：输入表1-3中的数据，对数据进行meshgrid ( )操作，使数组x和y产生网格，以绘制3D曲面；
Step2：x和y平面设置好后，用griddata ( )命令对矩形海域内的格点样条函数内插，得到相应的水深z；
Step3：运用mesh ( )作图命令，绘出该海域的地貌图，如图1-12（a）所示；
Step4：等高线易于观察地形地貌，使用contour ( )命令绘制等位线，并用clabel ( )命令标识等位线所处位置的深度值，如图1-12（b）所示；
Step5：对于特殊等位线在航海警示线上的应用，应特殊标记，再次运用contour ( )命令绘制警示线，运用gtext ( )命令标识不同吨位的警示线，如图1-12（c）所示，航行船只应该予以避让。
从图1-12中可看出，该地貌高低起伏，最大水位为10米，最低水位为4米。根据警示线标识情况，航运部门可对船只进行宏观调控，并对过往船只给以提示，避免发生安全事故。
程序如下。
clc,clear,close all
x=[129，140，103.5，88，185.5，195，105.5，157.5，107.5，77，81，162，162，117.5];
y=[7.5，141.5，23，147，22.5，137.5，85.5，-6.5，-81，3，56.5，-66.5，84，-33.5];
z=[4，8，6，8，6，8，8，9，9，8，8，9，4，9];
x1=75:0.5:200;
y1=-50:0.5:150;
[x2，y2]=meshgrid(x1，y1);
z2=griddata(x，y，z，x2，y2，'v4'); % 'v4'MATLAB 4格点样条函数内插
subplot(1，3，1)
mesh(x2，y2，z2);
hold on
clabel(C);
subplot(1，3，2)
C=contour(x2，y2，z2); %绘制等位线
clabel(C); %等位线标识
[C，d]=contour(x2，y2，z2);
clabel(C，d，'manual');
grid on
hold on
subplot(1，3，3)
contour(x2，y2，z2，[4，4，4.5，4.5，5，5，5.5，5.5]); %绘制不同吨位的警示线
>> title('警示线图');
>> hold on
xlabel('X')，ylabel('Y');
grid on
gtext('4吨的警戒线');
gtext('4.5吨的警戒线');
gtext('5吨的警戒线');
gtext('5.5吨的警戒线');
MATLAB以矩阵形式进行数值运算，运行效率较高，在工程上得到广泛应用。本章主要阐述了MATLAB基础知识，借助于实际应用讲述了符号变量在实际问题中的应用论证、插值模型的求解、等高线的绘制等问题。
自然界中常见到竞争种共存于同一环境，例如，水中浮游植物大都用类似的水体营养物进行光合作用，并能长期共存。有人认为这是某些物种共生的表现；也有人认为是由于优势种的增长受到捕食者的限制，因而劣势种未被排除。然而，这些解释仍不足以说明大量浮游植物长期共存的现象。有人则认为：种群竞争公式的推导是假定竞争种处于平衡状态，因而竞争排除也要求达到平衡状态，但水体环境在季节变化条件下总达不到平衡状态，因此物种得以长期共存。
种群（population）指在一定时间内占据一定空间的同种生物的所有个体。种群中的个体并不是机械地集合在一起的，而是彼此可以交配，并通过繁殖将各自的基因传给后代。种群是进化的基本单位。种群竞争模型，对于分析种群之间的竞争关系尤为重要。
学习目标：
（1）熟练掌握MATLAB编程表示方法；
（2）熟练运用MATLAB求解实际物理模型；
（3）熟练运用MATLAB求解微分方程；
（4）熟练使用MATLAB工具解决简单工程问题。
种群竞争模型是一个动态的过程，种群生存期间有着出生、死亡、迁入、迁出等问题。因此，种群数量较难确定，其种群竞争的数学模型只能通过反复的修正不断地完善，从而更加接近实际。本节就种群竞争微分方程模型的求解展开讨论。
设有甲、乙两种群，当它们独自生存时，数量演变服从Logistic规律：
其中，x(t)、y(t)分别为甲、乙两种群的数量，r1、r2为它们的固有增长率，n1、n2为它们的最大容量。
当两种群在同一环境中生存时，它们之间的一种关系是为了争夺同一资源而进行竞争。考察乙消耗有限的资源对甲的增长产生的影响，可以合理地将种群甲的方程修改为：
式（2.2）中s1的含义：对于供养甲的资源而言，单位数量乙（相对于n2）的消耗为单位数量甲（相对n1）消耗的s1倍。类似地，如果甲的存在也影响了乙的增长，则乙的方程应改为：
式（2.2）中s2的含义：对于供养乙的资源而言，单位数量甲（相对于n1）的消耗为单位数量乙（相对n2）消耗的s1倍。给定种群的初始值：
并给定参数r1、r2、s1、s2、n1、n2后，由式（2.1）～式（2.4）可确定两种群数量的变化规律。
（1）设r1=r2=1，n1=n2=100，s1=0.5，s2=2，x0=y0=10，计算x(t)、y(t)，画出它们的图形及相图x(t)、y(t)，说明时间t充分大时x(t)、y(t)的变化趋势。
解：对于微分方程的求解，首先建立微分方程函数多数情况下，用数值解代替代数解进行方程的模拟，自定义种群函数程序zhongqun( )：
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=1;r2=1;
n1=100;n2=100;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
%在此函数中，改变r1、r2、n1、n2、s1、s2的值，达到相关要求
针对题目中已知的初始条件，编写相应的MATLAB脚本文件程序，运行得图2-1、图2-2。
由图2-1、图2-2可知：在t=10时，x达到稳定值100，y达到稳定值0。
结论：时间t充分大时，x(t)、y(t)的值稳定在x=100,y=0。
图2-1的程序如下。
%绘制当r1=1，r2=1，n1=100，n2=100，s1=0.5，s2=2时的函数图像
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1;r2=1;s1=0.5;s2=2;n1=100;n2=100;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
图2-2的程序如下。
%绘制曲线向量解曲线
>>syms r1 r2 s1 s2 n1
r1=1;r2=1;s1=0.5;s2=2;n1=100;n2=100;
Xmin=0;
Xmax=140;
Ymin=0;
Ymax=100;
n=50;
% 计算切线矢量
>> [X，Y]=meshgrid(linspace(Xmin，Xmax，n)，linspace(Ymin，Ymax，n));
>> Fx=r1.*X.*(1-X./n1-s1.*Y./n2);
Fy=r2.*Y.*(1- s2.*X./n1-Y./n2);
Fx=Fx./(sqrt(Fx.^2+Fy.^2+1));
Fy=Fy./(sqrt(Fx.^2+Fy.^2+1));
%求解微分方程
>> options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
>> [T1，Y1]=ode45(@zhongqun，[0 50]，[10 10]，options);
>>% 绘制斜率场
hold on
grid on
box on
axis([Xmin，Xmax，Ymin，Ymax])
quiver(X，Y，Fx，Fy，0.5);
>>% 绘制解曲线
plot(Y1(:，1)，Y1(:，2)，'g'，'LineWidth'，2)
（2）改变r1、r2、n1、n2、x0、y0，维持s1、s2不变，绘出r1=1.2，r2=1.1，n1=200，n2=120，x0=y0=10时的函数图像及r1=0.9,r2=1.5,n1=500,n2=800，x0=y0=10时的函数图像。
解：同第（1）问，改变初始值，程序运行后依次如图2-3、图2-4所示。图2-3为r1=1.2，r2=1 . 1,n1=200，n2=120，x0=y0=10时的函数图像；图2-4 为r1=0.9，r2=1.5，n1=500，n2=800，x0=y0=10时的函数图像。
图2-3的程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=1.2;r2=1.1;
n1=200;n2=120;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
% 保持s1、s2不变，改变其他变量
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1.2;r2=1.1;s1=0.5;s2=2;n1=200;n2=120;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
图2-4的程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=0.9=1.5
n1=500;n2=800;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=0.9;r2=1.5;s1=0.5;s2=2;n1=500;n2=800;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
改变r1、r2、n1、n2、x0、y0，维持s1、s2不变，种群甲将占优势，而种群乙变为0。
综合结论：改变r、n的初始值，甲、乙种群的最终稳定状态不会改变，都是种群甲达到环境最大承载值，而种群乙变为0。参数r、n的初始值的改变仅会影响达到稳定的速度，不会改变优势种群甲的优势地位，即最终的稳定状态情况。
若s1=1.5，s2=0.7，绘出函数图像，如图2-5所示。
在图2-5中，当t=20时，y达到最大容量稳定值100，种群x变为零。当s1小而s2大时（s1与s2相差较大时，s1<1，s2>1），乙消耗甲的资源少，乙对甲影响小，同时甲消耗乙的资源多，甲对乙影响大，从而乙处于不利地位，甲处于有利地位，所以最后种群甲达到环境最大承载量，而种群乙则变为0。
相反，当s1大而s2小时（s1与s2相差较大时，s1>1，s2<1），乙消耗甲的资源多，乙对甲影响大，同时甲消耗乙的资源少，甲对乙影响小，乙处于有利地位，甲处于不利地位，所以最后种群乙达到环境最大承载量，而种群甲则变为0。
程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=1.5;r2=0.7
n1=100;n2=100;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1;r2=1;s1=1.5;s2=0.7;n1=100;n2=100;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
（3）试验当s1=0.8，s2=0.7时会有什么结果；当s1=1.5，s2=1.7时又会有什么结果。你能解释这些结果吗？
解：当s1=0.8，s2=0.7时，绘制图像如图2-6所示；当s1=1.5，s2=1.7时，绘制图像，如图2-7所示。
由这两个图可知：
①在s1<1，s2<1，且s1、s2相近时，由于双方的抑制作用数值较小，所以任何一方都无法占据绝对优势，所以双方均无法消灭对方，只能在最大承载量以下达到稳定状态，且受到影响小的（图2-6、图2-7中乙受到的影响小）稳定值大。
②当s1>1，s2>1，且s1、s2相近时，由于s1、s2均比较大，对于种群的抑制力非常大，以至于一旦竞争处于下风，就会受到很大的抑制作用而无法生存。s1、s2相近，仍然有一方会灭绝。但这是一个较长的过程，一方无法短时间内完全抑制对方，另一方因受到抑制，故增长也会减缓。虽然r、n初值均相等，但由于s1<s2，所以最后仍是甲受到的影响小，最终胜出。第（1）问中t=10时达到稳定，而此时t=20时才达到稳定，足足晚了10秒。
图2-6的程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=0.8;r2=0.7
n1=100;n2=100;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1;r2=1;s1=0.8;s2=0.7;n1=100;n2=100;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
图2-7的程序如下。
%自定义种群函数
function dy=zhongqun(t，y)
syms r1 r2 s1 s2 n1
%r、n赋予不同的参数时，有不同的解
r1=1.5;r2=1.7
n1=100;n2=100;
s1=0.5;s2=2;
dy=zeros(2，1);
dy(1)=r1*y(1)*(1-y(1)/n1-s1*y(2)/n2);
dy(2)=r2*y(2)*(1-s2*y(1)/n1- y(2)/n2);
运行的脚本文件程序如下。
>> x0=10;y0=10;
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun'，[0 50]，[x0 y0]，options);
grid on
axis equal
plot(T，Y(:，1)，'b-'，T，Y(:，2)，'r-')
title('r1=1;r2=1;s1=1.5;s2=1.7;n1=100;n2=100;x0=10;y0=10;')
h=legend('x(t)'，'y(t)'，2);
蓝鲸和长须鲸是两个生活在同一海域的相似的种群，因此认为它们之间存在竞争。估计蓝鲸的固有增长率每年为5%，长须鲸为每年8%。估计蓝鲸环境承载力（环境能够支持的鲸鱼的最大数量）为150 000 条，长须鲸为400 000 条。鲸鱼竞争的程度是未知的。过去约100 年剧烈的捕捞已经使鲸鱼数量减少，蓝鲸大约5 000 条，长须鲸大约70 000 条。蓝鲸是否会灭绝？
解：该问题是对2.1节第（1）、（2）、（3）问的具体化。根据实际统计数据，进行蓝鲸和长须鲸的种群竞争模型模拟计算。假设蓝鲸和长须鲸的增长情况仅与两者之间的竞争有关，与其他动物无关；不考虑环境改变带来的影响，环境承载力是稳定的；人类停止对鲸鱼的捕杀，鲸鱼按自然条件繁衍。
建立相应的数学模型：
由于蓝鲸和长须鲸的竞争是未知的，根据实际情况，当s1、s2取2时，种群间的相互影响已经非常大，所以估计s1、s2的区间均为（0，2）时，就可以很好地模拟实际情况。以 0.1 为步长，采用穷列举法，代入不同的s1、s2，求出最后蓝鲸与长须鲸的稳定状态数值。求解的图形如图2-8、图2-9所示。
不同情况下s1、s2的具体数值列表如下。
（1）蓝鲸与长须鲸的最后稳定状态都不为0时的s1、s2值如表2-1所示。
（2）当蓝鲸具有优势，最后稳态值为x=150 000，y=0时的s1、s2值如表2-2所示。
（3）当长须鲸具有优势，最后稳态值为x=0，y=400 000时的s1、s2值如表2-3所示。
由以上三表所示结果可知，由于蓝鲸的固有增长率每年为5%，长须鲸为每年8%，蓝鲸环境承载力为150 000条，长须鲸为400 000 条，所以长须鲸在相同的竞争程度条件下明显具有一定的优势。因而，当蓝鲸具有优势，最后稳态值为x=150 000，y=0时的s1、s2必定是s2较大、s1较小，而表2-2 中s1基本小于 1、s2基本大于 1，正好体现了在蓝鲸对长须鲸影响小、蓝鲸对长须鲸影响大的条件下，蓝鲸具有优势的实际情况。
同样，当长须鲸具有优势，最后稳态值为x=0，y=400 000时，根据实际，必然是s1较大（s1>1），而s2较小（s2<1）；同时，由于长须鲸本身对环境的适应力较强（体现在增长率和环境承载力较大），所以即使在相互抑制作用均较大（s1，s2>1）时，仍然是长须鲸处于优势地位。综合所述，s1>1时，无论s2取何值，长须鲸都具有优势。表2-3中，s1在1.1以上时，x=0，说明理论计算结果很好地符合实际情况。
通过以上分析可知，蓝鲸不一定毁灭，s1、s2的取值不同时，最后的稳定状态不同。
种群函数：
function dy=zhongqun4(t，y)
global E F B K
dy=zeros(2，1);
dy(1)=0.05*y(1)*(1-y(1)/150000-E(B)*y(2)/400000);
dy(2)=0.08*y(2)*(1-F(K)*y(1)/150000-y(2)/400000);
主程序：
global E F B K
E=0:0.1:2;
E=E';
F=0:0.1:2;
F=F';
B=1;K=1;a=1;b=1;c=1;d=1;
S=zeros(140，2); %记录当蓝鲸与长须鲸的最后稳定状态都不为0时的 s1，s2值
H=zeros(250，2); %记录蓝鲸具有优势，最后稳态值为x=150 000，y=0时的 s1，s2值
U=zeros(250，2); %记录长须鲸具有优势，最后稳态值为x=0，y=400000时的 s1，s2值
Num=zeros(441，4); %记录蓝鲸与长须鲸的最后稳态值
while B<22
K=1;
while K<22
options=odeset('RelTol'，1e-4，'AbsTol'，[1e-4 1e-5]);
[T，Y]=ode45('zhongqun4'，[0 2000]，[5000 70000]，options);
[m，n]=size(Y);
Num(a，1)=Y(m，1);
if Num(a，1)<1
Num(a，1)=0;
end
Num(a，2)=Y(m，2);
if Num(a，2)<1
Num(a，2)=0;
end
Num(a，3)=E(B);
Num(a，4)=F(K);
if ((Y(m，1)-1>0)&(Y(m，2)-1>0))==1
S(b，1)=E(B);
S(b，2)=F(K);
b=b+1;
end
if (Y(m，2)-1)<0
H(c，1)=E(B);
H(c，2)=F(K);
c=c+1;
end
if (Y(m，1)-1)<0
U(d，1)=E(B);
U(d，2)=F(K);
d=d+1;
end
a=a+1;
K=K+1;
end
B=B+1;
end
当s1、s2改变时，蓝鲸与长须鲸的量数：
>> [s1，s2]=meshgrid(Num(:，3)，Num(:，4));
>> LANJING=griddata(Num(:，3)，Num(:，4)，Num(:，1)，s1，s2，'v4');
CHANGXUJING=griddata(Num(:，3)，Num(:，4)，Num(:，2)，s1，s2，'v4');
>> mesh(s1，s2，LANJING)
>> title('s1和s2改变时，蓝鲸的数量变化图像');
mesh(s1，s2，CHANGXUJING);
title('s1和s2改变时，长须鲸的数量变化图像')
MATLAB是以矩阵形式进行数值运算，运行效率较高，在工程上得到广泛应用。本章主要阐述了微分方程在种群竞争模型中的应用，并利用MATLAB对种群竞争微分方程的应用求解。
食品是城市居民生活的重要元素，食品的零售价格与居民的生活是息息相关的。近年来，我国食品的零售价格指数受自然因素、市场因素的影响逐年增大；在国家宏观调控和市场的自动调控下，其也在一定范围内稳定变化。但是由于市场的随机变化性，商品的零售价格指数变化受市场上众多不确定因素共同影响，因而其变化也有不确定性，所以简单地根据其价格变化趋势对其价格进行预测显得不合理。
马尔可夫链（Markov）是数学中具有马尔可夫性质的离散时间随机过程。该过程中，在给定当前知识或信息的情况下，过去（即当期以前的历史状态）对于预测将来（即当期以后的未来状态）是无关的。因此马尔可夫链有众多的生物学应用，特别是人口过程，可以帮助模拟生物人口过程的建模，还可用来建模排队理论和统计学中的建模。隐蔽马尔可夫链模型还被用于生物信息学，用以编码区域或基因预测。本章将着重讲解 Markov 在食品价格趋势上的应用研究。
学习目标：
（1）熟练掌握分析和解决工程实际问题的思路；
（2）熟练掌握MATLAB编程的方法以及思路；
（3）熟练掌握矩阵在求解工程问题上的应用；
（4）熟练运用Markov求解实际物理模型；
（5）熟练使用MATLAB工具解决复杂的工程项目问题。
消费者物价指数（Consumer Price Index，CPI），也称消费价格指数，是反映与居民生活有关的产品及劳务价格统计出来的物价变动指标，通常作为观察通货膨胀水平的重要指标，是与人们生活密切相关的参考指标。
城市居民食品零售价格是消费者物价指数的重要组成部分。权威机构研究认为粮食生产、流通成本上涨一定会带动农产品价格总体上涨，特别是2011年异常的气候情况，导致生产成本大量增加，国际粮价对国内供需的影响，食品价格未来可能发生上涨。
2011年3月，CPI增幅达5.4%，创32个月来的新高，这使得年内的通货膨胀压力增大。表3-1是2010年3月5日—2011年3月25日城市居民食品零售价格数据。
注：其中42种食品数据来自中国统计年鉴网站。
试解决下面两个问题：
（1）根据表3-1建立数学模型，将所涉及食品适当分类，并分析每类食品的特点；
（2）根据表3-1建立数学模型，预测2011年4月、5月城市居民食品零售价格的走势。
针对食品价格每隔十天一个价格值，价格波动变化范围较小，实际价格图形为无规律的曲线，随机性很强，在此对两问题做出如下假设。
假设1：表3-1中的城市居民食品零售价格数据基本属实。
假设2：增长率在(≥0.03)的区间内，食品价格快速上升；增长率在（0，0.03）区间内，价格缓慢上升；增长率为0，食品价格相对不变；增长率在（−0.03，0）区间内，价格缓慢下降；增长率在(≤−0.03)的区间内，食品价格快速下降。
根据食品价格变化规律来看，通常有以下两种方法对食品价格变化的特性进行描述。
（1）数值描述法：直接用价格增长为描述对象，绘制价格—时间曲线。数值描述法可以直观地描述价格的变动趋势，但是由于各种食品的价格相差较大，且数据量很大，无法体现食品价格变化的整体规律。
（2）增长率描述法：计算同种食品相邻时间点的价格增长率，绘制不同食品的增长率—时间曲线。增长率描述法依次计算相邻时期各食品价格的增长率，虽不能直观地反映价格值的变化，但可以直观地反映价格水平的变化。
对序列x1,x2,x3,…,xn进行变换：
对问题中的42种食品进行分类，xi（i=1，…，42）依次表示题中相应的食品，选择了两种方法。
MATLAB聚类分析的主要步骤：
Step1：首先对数据进行分析，由于42种食品有3种量纲，在此，先不考虑花生油、大豆调和油和酱油、醋，从而横向保证了单位的统一，无须无量纲标准化，可直接对样品进行聚类分析。
数据标准化方法主要有以下3种。
（1）规范化方法
则新序列y1,y2,…,yn∈[ 0,1 ]且无量纲。一般的数据需要时都可以考虑先进行规范化处理。
（2）正规化方法
对序列x1,x2,x3,…,xn进行变换：
其中，
。则新序列y1,y2,…,yn的均值为0，而方差为1，且无量纲。
（3）归一化方法
对正项序列x1,x2,x3,…,xn进行变换： 
则新序列y1,y2,…,yn∈[ 0,1 ]且无量纲，并且显然有
。
归一化方法在确定权重时经常用到，再次选用规范化方法进行数据的量纲归一化处理。
Step2：计算不同变量之间的距离dij，根据需要选用欧氏（Euclidean）距离公式：
得到的距离矩阵是实对称矩阵，其主对角线元素为零。
Step3：根据需要选用Ward最小方差法的逐步归类方法进行聚类。
经过MATLAB聚类分析后得到图3-1。
clc,clear,close all
load('jlx.mat')
x=jlx;
BX=zscore(x); %标准化数据矩阵
Y=pdist(x); %用欧氏距离计算两两之间的距离
D=squareform(Y); %欧氏距离矩阵
Z=linkage(Y); %最短距离法
% T=cluster(Z,4);
[H,T]=dendrogram(Z,'colorthreshold','default');
从图3-1可知，MATLAB聚类分析的38种食品可分为四类，再加上未聚类的两类，总共可分为如下六类，具体的食品种类如下。
第一类共17种：芹菜、大白菜、油菜、黄瓜、萝卜、茄子、圆白菜、西红柿、土豆、胡萝卜、青椒、韭菜、芦柑、西瓜、香蕉、豆腐、食用盐；
第二类共12种：白砂糖、菜籽油、大豆油、鸡蛋、草鱼、鲤鱼、尖椒、豆角、蒜薹、苹果、红糖、鲜牛奶；
第三类共3种：鲜羊肉1、鲜羊肉2、鲜牛肉；
第四类共6种：鲜猪肉1、鲜猪肉2、活鸡、鸡肉、带鱼、绵白糖；
第五类共2种：花生油、大豆调和油；
第六类共2种：酱油、醋。
end
运用MATLAB对食品价格进行聚类分析，只是对各类食品的价格进行欧氏（Euclidean）距离的聚类，即各类食品价格之间的波动距离的最小值的聚合。例如，第一类食品的价格变动范围为1—3元，第二类食品的价格变动范围为 3—6 元，第三类食品的价格变动范围为15—17元，第四类食品的价格变动范围为 7—12元，第五类食品的价格变化范围为50—120元，第六类食品的价格变动范围为3.5—5.5元。这六类食品的价格变化范围虽然在一定的范围内，但是仍不能形象直观地反映各类食品的价格和价格增长率的变化情况，需对这六类进行图示分析。
1．第一类食品分析
第一类食品价格变化程序：
%x=[ 42种食品各个时间的价格];在MATLAB中的Workspace中载入
%x1(时间),y1、y2、y3、y4、y5、y6（分类后的食品的物价）均在MATLAB中的Workspace中载入
%画出已分类的食品的价格—时间曲线
clc,clear,close all
load('x.mat')
load('y1.mat')
%%
for i=1:17
plot(x1,y4(:,i));
hold on
end
xlabel('月日');ylabel('价格')
第一类食品价格增长率变化程序：
%提取经过聚类后的已分类的食品,画出增长率—时间曲线
%y111、y222、y333、y444、y555、y666为分类后的各类食品的相邻时间的增长率
%第一类的增长率
for i=2:39
y44((i-1),:)=y4(i,:)-y4((i-1),:);
end
for i=1:38
for j=1:17
y444(i,j)=y44(i,j)/y4(i,j);
end
end
%第一类食品的增长率—时间曲线
for i=1:17
subplot(4,5,i);
plot(x1,y444(:,i));
end
该类食品以西红柿为代表，即该类其他食品价格紧靠西红柿的价格变化所得的欧氏距离总和最小，西红柿的价格变动情况如图3-2 中的粗线所示，价格按先升高后降低的趋势变化。
从图3-3中可看出，芹菜、大白菜、油菜、黄瓜、萝卜、茄子、圆白菜、西红柿、土豆、胡萝卜、青椒、韭菜等的价格增幅变化同步，随时间的变化呈现明显的变化；芦柑、西瓜、豆腐、食用盐的增幅变化不明显，先保持稳定后增长，特别是食盐基本无任何变化。
2．第二类食品分析
第二类食品价格变化程序：
clc,clear,close all
load('x.mat')
load('y1.mat')
%%
for i=1:12
plot(x1,y3(:,i));
hold on
end
xlabel('月日');ylabel('价格')
第二类食品价格增长率变化程序：
%第二类的增长率
for i=2:39
y33((i-1),:)=y3(i,:)-y3((i-1),:);
end
for i=1:38
for j=1:12
y333(i,j)=y33(i,j)/y3(i,j);
end
end
%第二类食品的增长率—时间曲线
for i=1:12
subplot(3,4,i);
plot(x1,y333(:,i));
end
该类食品以红糖为代表，如图3-4 中的粗实线所示，各食品的价格变化仍无规律。但从图3-5 中看出，菜籽油、大豆油价格增幅情况相当，红糖、白砂糖价格变化仍相似。从该大类的价格增幅图得出，各种食品的价格增长率变化迥异。
3．第三类食品分析
第三类食品价格变化程序：
clc,clear,close all
load('x.mat')
load('y1.mat')
%%
for i=1:3
plot(x1,y2(:,i));
hold on
end
xlabel('月日');ylabel('价格')
第三类食品价格增长率变化程序：
%第三类的增长率
for i=2:39
y22((i-1),:)=y2(i,:)-y2((i-1),:);
end
for i=1:38
for j=1:3
y222(i,j)=y22(i,j)/y2(i,j);
end
end
%第三类食品的增长率—时间曲线
for i=1:3
subplot(1,3,i);
plot(x1,y222(:,i));
end
该类食品为鲜羊肉去骨、鲜牛肉、鲜羊肉带骨，价格变化情况大致一致，在2011年1到3月肉价格上涨加快，且价格变化以鲜牛肉为中心所得的欧氏距离值最小。在图3-7中，价格增幅也较明显，先保持稳定后快速上涨再下降保持稳定，这也符合市场肉类价格的变化规律。
4．第四类食品分析
第四类食品价格变化程序：
clc,clear,close all
load('x.mat')
load('y1.mat')
%%
for i=1:6
plot(x1,y1(:,i));
hold on
end
xlabel('月日');ylabel('价格')
第四类食品价格增长率变化程序：
load('x1-38.mat')
%第四类的增长率
for i=1:38
x1(i)=i;
end
for i=2:39
y11((i-1),:)=y1(i,:)-y1((i-1),:);
end
for i=1:38
for j=1:6
y111(i,j)=y11(i,j)/y1(i,j);
end
%第四类食品的增长率—时间曲线
for i=1:6
subplot(3,2,i);
plot(x1,y111(:,i));
end
从图3-8 中可看出，该类食品为鲜猪肉精瘦肉、鲜猪肉肋条肉、活鸡、鸡肉、带鱼、绵白糖，其中欧式距离值最小以鲜猪肉肋条肉为中心。从图3-9 中可看出，绵白糖的价格变化规律与鲜猪肉精瘦肉、鲜猪肉肋条肉、活鸡、鸡肉、带鱼的变化步调有明显区别，先保持稳定后快速增长，且保持稳定的时间很长。
5．第五类食品分析
第五类食品价格变化程序：
clc,clear,close all
load('x.mat')
load('y1.mat')
%%
for i=1:2
plot(x1,y5(:,i));
hold on
end
xlabel('月日');ylabel('价格')
第五类食品价格增长率变化程序：
%第五类的增长率
for i=2:39
y55((i-1),:)=y5(i,:)-y5((i-1),:);
end
for i=1:38
for j=1:2
y555(i,j)=y55(i,j)/y5(i,j);
end
end
%第五类食品的增长率—时间曲线
for i=1:2
subplot(1,2,i);
plot(x1,y555(:,i));
end
从图3-10、图3-11 中可看出，该类食品以花生油、大豆调和油为一类，其中花生油的价格变化范围为95—120元，居高不下，价格增长幅度先保持稳定后上涨；大豆调和油的价格在46—65元之间变化，价格走势与花生油大致相当，增长率变化先保持稳定后快速增长。
6．第六类食品分析
第六类食品价格变化程序：
clc,clear,close all
load('x.mat')
load('y1.mat')
%%
for i=1:2
plot(x1,y6(:,i));
hold on
end
xlabel('月日');ylabel('价格')
第六类食品价格增长率变化程序：
%第六类的增长率
for i=2:39
y66((i-1),:)=y6(i,:)-y6((i-1),:);
end
for i=1:38
for j=1:2
y666(i,j)=y66(i,j)/y6(i,j);
end
end
%第六类食品的增长率—时间曲线
for i=1:2
subplot(1,2,i);
plot(x1,y666(:,i));
end
从图3-12、图3-13中可看出，该类食品由酱油、醋组成，酱油、醋的价格变化趋势相当，先保持稳定后有所下降，再上升保持稳定；在价格增长率图中，两者价格增幅在时间序列上呈现延迟性的变化相当，此类聚合在一起符合实际情况。
综合上述分析可知，采用简单的按照欧氏（Euclidean）距离对各食品进行价格上的聚合分类，从而划分食品种类，在一定程度上不能反映客观规律。国家对食品价格进行宏观上的调控，主要针对食品大类进行调控，如肉类及其制品、鲜菜类、鲜果类、鱼类等等，而不是按照笼统的价格分类。
食品价格受国家宏观调控，国家根据市场食品种类进行价格调控，因此根据“国家食品分类系统”、表3-1中的各食品数据（来自中国国家统计年鉴数据）以及中国国家统计年鉴数据查询结果等资料，把42种食品分为食用油类、肉类及其制品、鱼类、蔬菜类、水果类、调味品类、奶类七大类来进行样本属性分析。
根据食品属性，即食用油类、肉类及其制品、鱼类、蔬菜类、水果类、调味品类、奶类等分类，具体如下所示。
第一类 油类共4 种：菜籽油、大豆油、花生油、大豆调和油；
第二类 肉类及蛋类共8种：鲜猪肉1、鲜猪肉2、活鸡、鸡肉、鲜羊肉1、鲜羊肉2、鲜牛肉、鸡蛋；
第三类 水产品类共3 种：草鱼、鲤鱼、带鱼；
第四类 鲜菜类共15 种：芹菜、大白菜、油菜、黄瓜、萝卜、茄子、圆白菜、西红柿、土豆、胡萝卜、青椒、尖椒、韭菜、豆角、蒜薹；
第五类 鲜果类共4 种：香蕉、苹果、西瓜、芦柑；
第六类 调味品类共7 种：食用盐、红糖、白砂糖、绵白糖、酱油、醋、豆腐；
第七类 奶类1 种：鲜牛奶。
根据食品属性分类的食品，绘制每类食品的价格增长率曲线图（图3-14、图3-15、图3-16、图3-17、图3-18、图3-19、图3-20中对称的四条虚线的价格增长率为±0.03和±0.05）。
1．第一类食品价格增长率分析
%z1、z2、z3、z4、z5、z6为分类后的各食品的物价，且在Workspace中直接载入
%z11、z22、z33、z44、z55、z66为分类的各食品的相邻价格差值，且在Workspace中直接载入
%z111、z222、z333、z444、z555、z666为各种食品的相邻时间价格的增长率，且在Workspace中直接载入
%第一类食品的增长率图
for i=1:39
x2(i)=i;
end
for i=1:38
for j=1:4
z111(i,j)=z11(i,j)/z1(i,j);
end
end
for i=1:4
subplot(2,2,i);
plot(x1,z111(:,i),x2,0.03*ones(1,39),'r--',x2,-0.03*ones(1,39),'r--',x2,0.05*ones(1,39),'r--',x2,-0.05*ones(1,39),'r--' );hold on
end
从图3-14中可知，这四种油类的变化趋势相当，都是先保持不变后上涨。这四种食品价格都是平稳或者缓慢增长的，且菜籽油、大豆油、花生油、大豆调和油价格都是在2010年10月25日至11月25日左右有个增幅较大的变化。
2．第二类食品价格增长率分析
%第二类食品的增长率图
for i=1:38
for j=1:8
z222 (i,j)=z22(i,j)/z2(i,j);
end
end
for i=1:8
subplot(3,3,i);
plot(x1,z222(:,i),x2,0.03*ones(1,39),'r--',x2,-0.03*ones(1,39),'r--',x2,0.05*ones(1,39),'r--',x2,-0.05*ones(1,39),'r-');
hold on
end
从图3-15 中可知，猪肉和牛肉的价格是上下波动变化的；而羊肉的价格只是在 2010年11月、12月有较大幅度的增加；鸡肉和鸡蛋的价格基本是在平稳变化的，在2010年，还是有显著的价格上升，但持续的时间不长。
3．第三类食品价格增长率分析
%第三类食品的增长率图
for i=1:38
for j=1:3
z333 (i,j)=z33(i,j)/z3(i,j);
end
end
for i=1:3
subplot(1,3,i);
plot(x1,z333(:,i),x2,0.03*ones(1,39),'r--',x2,-0.03*ones(1,39),'r--',x2,0.05*ones(1,39),'r--',x2,-0.05*ones(1,39),'r--');hold on
end
从图3-16中可知，对于水产品而言，鱼的价格没有方向性，下降和上升变化幅度大，不同鱼的变化步调不一致。鱼类产品的价格随季节性变化很大，例如夏天鱼的价格相对便宜，而冬季就相对贵了。它还和该地方的经济水平、人们的喜好和其生产成本及流通成本有很明显的联系。
4．第四类食品价格增长率分析
%第四类食品的增长率图
for i=1:38
for j=1:15
z444 (i,j)=z44(i,j)/z4(i,j);
end
end
for i=1:15
subplot(4,4,i);
plot(x1,z444(:,i),x2,0.03*ones(1,39),'r--',x2,-0.03*ones(1,39),'r--',x2,0.05*ones(1,39),'r--',x2,-0.05*ones(1,39),'r--');hold on
end
从图3-17中可知，鲜菜类食品的价格快速增长和快速下降的概率是很大的，说明蔬菜的价格上下波动的幅度特别大。因为蔬菜类食品受天气、自然灾害、季节、生产成本、流通成本、国家物价部门的宏观调控、全国的蔬菜价格、国际蔬菜价格的影响非常大，所以才会造成蔬菜价格的这种变化。
5．第五类食品价格增长率分析
%第五类食品的增长率图
>> for i=1:38
for j=1:4
z555 (i,j)=z55(i,j)/z5(i,j);
end
end
for i=1:4
subplot(2,2,i);
plot(x1,z555(:,i),x2,0.03*ones(1,39),'r--',x2,-0.03*ones(1,39),'r--',x2,0.05*ones(1,39),'r--',x2,-0.05*ones(1,39),'r--');hold on
end
从图3-18中可知，鲜果类食品价格增长率的变化也非常大，说明瓜果类食品的价格上下波动的幅度也较大。鲜果类食品受天气、季节、生产成本、流通成本、国家物价部门的宏观调控、全国的鲜果类价格、国际鲜果类价格的影响非常大，所以才会造成鲜果价格的快速下降和快速增长。鲜果类的食品不宜长时间储藏，所以有关生产部门要控制水果的生产量，广泛推广生态瓜果，以减少蔬菜对气候和季节的依赖，而且国家物价部门也在适当控制鲜果类价格，使之相对稳定地变化。
6．第六类食品价格增长率分析
%第六类食品的增长率图
>> for i=1:38
for j=1:7
z666(i,j)=z66 (i,j)/z6(i,j);
end
end
for i=1:7
subplot(2,4,i);
plot(x1,z666(:,i),x2,0.03*ones(1,39),'r--',x2,-0.03*ones(1,39),'r--',x2,0.05*ones(1,39),'r--',x2,-0.05*ones(1,39),'r--');
hold on
end
从图3-19中可知，生活调味品的增长率变化是相对稳定的，增长幅度不大，如豆腐价格上涨后，后期价格一直保持不变；其他几种调味品价格变化也是先保持相对不变后有一定的增幅。
7．第七类食品价格增长率分析
%第七类食品的增长率图
for i=1:38
for j=1
z777(i,j)=z77 (i,j)/z7(i,j);
end
end
for i=1
plot(x1,z777(:,i),x2,0.03*ones(1,39),'r--',x2,-0.03*ones(1,39),'r--',x2,0.05*ones(1,39),'r--',x2,-0.05*ones(1,39),'r-');hold on
end
从图3-20中可知，奶类食品的价格增长率先上涨后保持稳定，价格变化不很明显。鲜牛奶价格主要依赖奶牛的产奶量，而国内牛奶的价格在很大程度上依赖于国际牛奶价格，由于国外金融危机的影响，国外农业生产加强，国外牛奶的供给量保持了国内市场的稳定，故鲜牛奶价格表现为先上涨后保持稳定，总体表现为价格是相对稳定的。
综合上述分析，汇总各类食品价格特点，如表3-2所示。
注：各类食品的特点由根据2010年3月5日—2011年3月25日相邻增长率出现情况一致的次数总计得到。总计和越大，越能代表该种食品的特点。
近年来，我国食品的零售价格指数受自然因素、市场因素的影响逐年变大。但是在国家宏观调控和市场的自动调控下，其价格指数也在一定范围内稳定变化。由于市场的随机变化性，商品的零售价格指数变化受市场上众多不确定因素共同影响，因而其变化也有不确定性，所以简单地根据其价格变化趋势对其价格进行拟合或者对其增长率进行拟合，以此来预测后面两个月的价格变化趋势，不能很好地反映市场随机性和影响因素的不确定性。本节采用的是马尔可夫链模型，用MATLAB软件进行数据分析，预测42种食品2011年4月5日—2011年5月25日的价格增幅。
食品价格预测模型基本假设如下：
（1）假设表3-1所提供数据基本属实可靠；
（2）假设食品价格的变化客观反映了 CPI、食品消费价格指数产品、流通成本、生产资料、通货膨胀的变化；
（3）假设2011年4、5月无自然灾害等突发事件的影响，环境因素的影响基本和前期一致；
（4）在预测的6组数据中，前一组数据对其后的数据预测没有影响。
Bi：价格增长率划分为5 段，分别用Bi表示(i=1,2,3,4,5)，即快速上升(≥0.03)、缓慢上升(0,0.03)、相对不变为0、缓慢下降(−0.03,0)、快速下降(≤−0.03)；具体表示方法如下：
B1：快速上升，在程序中用“2”表示；
B2：缓慢上升，在程序中用“1”表示；
B3：相对不变，在程序中用“0”表示；
B4：缓慢下降，在程序中用“-1”表示；
B5：快速下降，在程序中用“-2”表示。
(E1,E2,…,Ei)：状态向量，用来表示状态变量(i=1,2,3,4,5)，即对应问题一中的快速上升、缓慢上升、相对不变、缓慢下降、快速下降五种状态；
Ei(0)：初始状态向量；
E：一步转移概率矩阵；
Eij：从状态i到状态 j的一步转移概率；
Cij：42 种食品空间中从状态i一步转移到状态 j的食品个数；
Di：各种食品价格系统中处于状态i时的样本个数。
问题分析部分对问题二的分析显示，由于市场的随机变化性，食品的零售价格指数变化受市场上众多不确定因素共同影响，因而其变化也有不确定性。
针对未来一段时间内的价格，既对历史数据有一定的参考性，又具有一定的随机变化性的特点，在此，引入动态经济学模型中的马尔可夫链的概念来对4、5月的价格增长趋势进行预测。
本模型中马尔可夫链的数据是42种食品在2010年3月5日—2011月3月25日之间每隔10天的价格增长率。在第一问中已算出，每种产品都有38组数据，如表3-3所示。
本书中将以第一种食品菜籽油为例，取2010年3月5日—2011年3月15日之间37组增长率进行计算，将第38组增长率与模型的预测数据进行比对，从而验证模型的合理性。
在MATLAB程序中，主要参见以下步骤。
Step1：求出各食品数据的增长率。
Step2：编程，使快速上升(≥0.03)，在程序中用“2”表示；缓慢上升( 0,0.03)，在程序中用“1”表示；相对不变0，在程序中用“0”表示；缓慢下降(−0.03,0)，在程序中用“-1”表示；快速下降(≤−0.03)，在程序中用“-2”表示。得到各类食品在不同时期的Bi。
Step3：统计相邻增长率变化一致的次数Cij。就菜籽油而言，在2010年3月5日—2011年3月15日，其增长率为先快速增长后相对不变出现的次数为‘20’=1。
Step4：对每一种食品，分五类进行叠加求和，分别为“先快速上升”、“先缓慢上升”、“先相对不变”、“先缓慢下降”、“先快速下降”，得到Di。
Step5：求出一步转移概率矩阵E，即Cij每行的每五个数据依次与对应的Di相除，得到概率矩阵E的一行。
Step6：使用求出的一步转移概率矩阵E，求出不同时期的状态概率。
针对Di，由MATLAB程序得到：
clc,clear,close all
load('x.mat')
for j=1:42
for i=1:38
a(j,i)=(x(j,i+1)-x(j,i))/x(j,i); %求增长率
end
end
b=zeros(42,38);
c=zeros(42,25);
% 分别为不同增长率赋值
for j=1:42
for i=1:38
if a(j,i)>=0.03
b(j,i)=2;
elseif (a(j,i)<0.03)&&(a(j,i)>0)
b(j,i)=1;
elseif a(j,i)==0
b(j,i)=0;
elseif a(j,i)>-0.03 && a(j,i)<0
b(j,i)=-1;
elseif a(j,i)<-0.03
b(j,i)=-2;
end
end
end
% 统计相连增长率值特征
for j=1:42
for i=1:36
if (b(j,i)==2&&b(j,i+1)==2)
c(j,1)=c(j,1)+1;
elseif(b(j,i)==2&&b(j,i+1)==1)
c(j,2)=c(j,2)+1;
elseif(b(j,i)==2&&b(j,i+1)==0)
c(j,3)=c(j,3)+1;
elseif(b(j,i)==2&&b(j,i+1)==-1)
c(j,4)=c(j,4)+1;
elseif(b(j,i)==2&&b(j,i+1)==-2)
c(j,5)=c(j,5)+1;
elseif(b(j,i)==1&&b(j,i+1)==2)
c(j,6)=c(j,6)+1;
elseif(b(j,i)==1&&b(j,i+1)==1)
c(j,7)=c(j,7)+1;
elseif(b(j,i)==1&&b(j,i+1)==0)
c(j,8)=c(j,8)+1;
elseif(b(j,i)==1&&b(j,i+1)==-1)
c(j,9)=c(j,9)+1;
elseif(b(j,i)==1&&b(j,i+1)==-2)
c(j,10)=c(j,10)+1;
elseif(b(j,i)==0&&b(j,i+1)==2)
c(j,11)=c(j,11)+1;
elseif(b(j,i)==0&&b(j,i+1)==1)
c(j,12)=c(j,12)+1;
elseif(b(j,i)==0&&b(j,i+1)==0)
c(j,13)=c(j,13)+1;
elseif(b(j,i)==0&&b(j,i+1)==-1)
c(j,14)=c(j,14)+1;
elseif(b(j,i)==0&&b(j,i+1)==-2)
c(j,15)=c(j,15)+1;
elseif(b(j,i)==-1&&b(j,i+1)==2)
c(j,16)=c(j,16)+1;
elseif(b(j,i)==-1&&b(j,i+1)==1)
c(j,17)=c(j,17)+1;
elseif(b(j,i)==-1&&b(j,i+1)==0)
c(j,18)=c(j,18)+1;
elseif(b(j,i)==-1&&b(j,i+1)==-1)
c(j,19)=c(j,19)+1;
elseif(b(j,i)==-1&&b(j,i+1)==-2)
c(j,20)=c(j,20)+1;
elseif(b(j,i)==-2&&b(j,i+1)==2)
c(j,21)=c(j,21)+1;
elseif(b(j,i)==-2&&b(j,i+1)==1)
c(j,22)=c(j,22)+1;
elseif(b(j,i)==-2&&b(j,i+1)==0)
c(j,23)=c(j,23)+1;
elseif(b(j,i)==-2&&b(j,i+1)==-1)
c(j,24)=c(j,24)+1;
elseif(b(j,i)==-2&&b(j,i+1)==-2)
c(j,25)=c(j,25)+1;
end
end
end
d=zeros(42,5);
% 累加求和
for i=1:42
for j=1:25
if(j<6)
d(i,1)=d(i,1)+c(i,j);
elseif(j>5&&j<11)
d(i,2)=d(i,2)+c(i,j);
elseif(j>10&&j<16)
d(i,3)=d(i,3)+c(i,j);
elseif(j>15&&j<21)
d(i,4)=d(i,4)+c(i,j);
else
d(i,5)=d(i,5)+c(i,j);
end
end
end
运行程序得到C和D矩阵数值，整理得表3-4。
从而获得一步转移概率矩阵E：
%一步转移概率矩阵
f=b(:,37);
e=zeros(5,5);
for i=1:42
for j=1:25
if(j<6)
if(d(i,1)==0)
e(5,j)=0;
else
e(1,j)=c(i,j)/d(i,1);
end
elseif(j>5&&j<11)
if(d(i,2)==0)
e(5,j-5)=0;
else
e(2,j-5)=c(i,j)/d(i,2);
end
elseif(j>10&&j<16)
if(d(i,3)==0)
e(5,j-10)=0;
else
e(3,j-10)=c(i,j)/d(i,3);
end
elseif(j>15&&j<21)
if(d(i,4)==0)
e(5,j-15)=0;
else
e(4,j-15)=c(i,j)/d(i,4);
end
else
if(d(i,5)==0)
e(5,j-20)=0;
else
e(5,j-20)=c(i,j)/d(i,5);
end
end
end
g=zeros(i,5); %预测4、5月的增长率
if(f(i,1)==2)
h=[1 0 0 0 0 ]*e
for k=1:6
h=h*e
end
elseif(f(i,1)==1)
g(i,:)=[0 1 0 0 0 ];
h=g(i,:)*e
for k=1:6
h=h*e
end
elseif(f(i,1)==0)
g(i,:)=[0 0 1 0 0 ];
h=g(i,:)*e
for k=1:6
h=h*e
end
elseif(f(i,1)==-1)
g(i,:)=[0 0 0 1 0 ];
h=g(i,:)*e
for k=1:6
h=h*e
end
elseif(f(i,1)==-2)
g(i,:)=[0 0 0 0 1 ];
h=g(i,:)*e
for k=1:6
h=h*e
end
end
end
如果目前预测对象处于状态Bi(i=1,2,3,4,5)，这时Eij就描述了目前状态Bi在未来将转向状态Bj(j=1,2,3,4,5)的可能性。按最大概率原则，选择（Ei1，Ei2，Ei3，Ei4，Ei5）中最大者对应的状态即为预测结果。
由于2011年3月5日—2011年3月15日的菜籽油价格的增长率状态为E3，即相对稳定。而经由一次转移到达5 种状态的概率分别为：E31=0.083 3，E32=0.208 3，E33=0.625 0，E34=0.083 3，E35=0。
而且，E31，E32，E34，E35与E33相比，均相差较大。
所以预测所得结果显示：2011年3月15日—2011年3月25日的价格会继续保持稳定，且增长幅度为0。
将预测结果与实际结果相比较可知，2011年3月15日—2011年3月25日的菜籽油价格增长状态为相对稳定，且增长幅度也为0，因此可以说此预测结果是准确的。
同时，在马尔可夫运算过程中，不同时期的状态概率用状态向量表示。公式为Ei( n )=Ei( n−1) E，按此公式，也可以计算出2011年3月15日—2011年3月25日的价格变动趋势，其中E为上一步转移概率矩阵。
对每个状态向量，均取其中最大的概率值，则以上结果表明：在未来两个月内，菜籽油的价格趋于稳定状态的概率较呈现其他状态的概率大得多。因此，可以说价格在2011年4、5月相对稳定。
对于其他41种食品，将依然按照上述马尔可夫链的计算方法计算，依次得出相关数据和预测数据。整理输出的结果如表3-5所示。
注：快速上升在程序中用“2”表示；缓慢上升在程序中用“1”表示；相对不变在程序中用“0”表示；缓慢下降在程序中用“-1”表示；快速下降在程序中用“-2”表示。
结果分析：
基于Markov食品物价趋势预测，使用概率矩阵完成预测2011年4、5月城市居民食品零售价格走势，误差很小。由于食品的价格增长受诸多不确定因素的影响，不能笼统地进行拟合计算分析。使用概率矩阵是根据下一次价格所出现的状态的最大可能概率来进行预测，可靠性高。
注：4—5月份预测数据结果如下所示。受篇幅所限，中间部分数据已删除，该数据仅作为参考。
1菜籽油
0.0833 0.2083 0.6250 0.0833 0
0.0799 0.1877 0.6803 0.0521 0
0.0833 0.1952 0.6648 0.0567 0
0.0832 0.1942 0.6673 0.0554 0
0.0833 0.1945 0.6666 0.0556 0
0.0833 0.1944 0.6667 0.0555 0
0.0833 0.1944 0.6667 0.0556 0
2大豆油
0 0.2500 0.6500 0.1000 0
0.0250 0.2625 0.6075 0.1050 0
0.0263 0.2819 0.5891 0.1028 0
0.0282 0.2863 0.5855 0.1000 0
0.0286 0.2891 0.5837 0.0986 0
0.0289 0.2902 0.5831 0.0978 0
0.0290 0.2908 0.5828 0.0974 0
……
42鲜牛奶
0.0303 0.0606 0.9091 0 0
0.0275 0.0551 0.9174 0 0
0.0278 0.0556 0.9166 0 0
0.0278 0.0556 0.9167 0 0
0.0278 0.0556 0.9167 0 0
0.0278 0.0556 0.9167 0 0
0.0278 0.0556 0.9167 0 0
注：每类食品输出数据的第一行为预测的2011-3-25的增长状态，后六行分别为2011-4-5、2011-4-15、2011-4-25、2011-5-5、2011-5-15、2011-5-25 的增长状态。其中每行的最大概率值表示该期最可能出现的增长状态。
本章针对食品价格波动进行分析，采用马尔可夫链模型对食品零售价格波动情况进行预测。问题一中，根据附录建立数学模型，将所涉及食品适当分类，由此分析每类食品的特点；然后分析了两种方法，分别是根据聚类分析法分类分析和根据食品的属性分类，画出了价格-时间图和增长率-时间图，并将其增长率划分为快速上升、缓慢上升、相对不变、缓慢下降、快速下降5类，对两种方法的图形进行了观察分析；最后给出了一个关于食品价格增长率特征的表，即为每类食品的特点。问题二中，使用马尔可夫链模型进行了数据分析，得出了42种食品2011年4月5日—2011年5月25日价格增幅的数据和42种食品在此期间的价格。因此，借助于强大的矩阵处理，用MATLAB求解工程问题显得更加贴合。
