Matlab在数学建模中的应用

1、第七章 Matlab在数学建模中的应用这一部分介绍数学建模中经常使用的软件：Matlab软件。Matlab是Matrix Laboratory的缩写，目前它已经成为国际上最流行的科学与工程计算的软件工具之一。现在Matlab 语言的功能越来越强大，集数值计算、符号计算和图形可视化于一体。尤其体现在运用简单而直接的符号代数方法来表示关系式，并执行运算，从而使用户可以将大部分精力集中在运算逻辑的推理上，而不必在繁杂的运算上耗费太多的精力。同时，Matlab提供了大量的函数解决数学模型问题。§7.1 Matlab的基础操作7.1.1 Matlab 的进入与界面 一、安装与启动：双击Matl

2、ab 的安装包，同其它应用软件类似，按照安装向导或提示进行安装即可。成功安装后，在Windows 桌面上就会出现Matlab图标，双击图标，就进入Matlab界面。二、界面（四个窗口） （1）命令窗口（Command Window）：在命令窗口中可以直接输入命令行，以实现计算或绘图功能。（2）工作空间 （Workspace）：显示当前Matlab的内存中使用的变量的信息，包括变量名、变量数组大小、变量字节大小和变量类型。 （3）命令历史（Command History）：显示所有执行过的命令。利用该窗口，一方面可以查看曾执行过的命令；另一方面，可以重复利用原来输入的命令行，这只需在命令历史窗口

3、中直接双击某个命令行，就可执行该命令行。 （4）当前目录(Current Folder)：该窗口显示当前工作目录下所有文件的文件名和文件类型，可以在窗口上方的小窗口中修改工作目录。图1 Matlab界面三、运行方式（两种） （1）命令行方式：通过直接在命令窗口输入命令行来实现计算或作图功能。不足：在处理比较复杂的问题和大量数据时相当困难。 （2）M文件方式：先在一个以 m为扩展名的M文件中输入一系列数据和命令，然后让Matlab执行这些命令。两种格式：i) 脚本M文件：命令的简单叠加； ii) 函数 M文件。建立 M 文件的方法：File New M-File，编缉程序，再以m为扩展名存储。

4、7.1.2 变量的Matlab表示变量由Matlab语句直接输入，不需要提前声明，输入格式为：“变量=表达式”例如，需要得到一个新变量t，其值为1.34，命令及结果如下：>> t=1.34t = 1.34注意：Matlab的命令行用“>>”表示开始输入，输入后回车即得结果。如果在命令行以英文的分号“;”结束，Matlab同样会进行相应的运算，但不显示运算结果。如：>> t=1.34;事实上大部分数值计算问题可归结为矩阵的计算。通常矩阵与数组的意义相同，都是指含有m行n列数字的矩形结构。在Matlab语言中表示一个矩阵式很容易，例如，对于矩阵，采用下面的命令：

5、>> A=1,2,3;4 5 6;7 8 9 A= 123 456 789注意到矩阵的内容由方括号括起来表示，而方括号内的分号表示矩阵的换行，逗号或空格表示同一行矩阵元素间的分隔。给出了上面的命令，就可以在MATLAB的工作空间中建立一个矩阵变量A。再比如，行向量和列向量可以类似下面的方法直接输入：>> A1=1 2 3;>> A2=4;5;6;7;我们可以用下标来访问某个或某些矩阵元素。例如，A矩阵的第2行第3列元素可以利用下面的命令得到：>> A(2,3)ans = 6注意这里的ans是缺省的变量名，用来显示运行的结果。我们还可以通过下面的命

6、令提取矩阵A的子矩阵：>> A(1,3,2,3)ans = 2 3 8 9可以看出，由上面的语句可以提取出矩阵A的第1,3行和第2、3列，而后构成了一个子矩阵。>> A(:,2,3)ans = 2 3 5 6 8 9在上面的语句中用的了冒号“:”可以提取出矩阵A的所有行和第2、3列，而后构成了一个子矩阵。MATLAB语言定义了独特的冒号表达式来给行向量赋值，其基本使用格式为：as1:s2:s3；其中，sl为起始值，s2为间距，s3为终止值。如果s2的值为负值，则要求s1的值大于s3的值，否则结果为一个空向量。如果省略了s2的值，则步长为默认值1。例如：>>

7、a=1:0.1:2;>> b=1:8;上面的语句分别将向量1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0赋值给变量a，将向量1 2 3 4 5 6 7 8赋值给变量b。还有一些函数命令可以得到矩阵。例如 MATLAB用函数eye(m,n)来表示m行n列的单位矩阵，用函数eye(n,n)来表示n阶单位方阵，还有可以用zeros(m,n)表示m行n列的零矩阵。>> I=eye(4,4);%得到四阶单位矩阵>> J=eye(4,3)J = 1 0 0 0 1 0 0 0 1 0 0 0>> P=zeros(3,4

8、); %得到三行四列的零矩阵注：上述命令中的百分号“%”是注释命令，它后面的文字为注释，不参与运算。7.1.3 变量的相关运算MATLAB提供了许多数学函数，如 abs绝对值函数、sin为正弦函数、sqrt为开平方、exp为以e为底的指数函数、log为自然对数函数等。MATLAB所支持的部分常用函数如下表： 函 数名 称函 数名 称sin(x)正弦函数asin(x)反正弦函数cos(x)余弦函数acos(x)反余弦函数tan(x)正切函数atan(x)反正切函数abs(x)绝对值max(x)最大值min(x)最小值sum(x)元素的总和sqrt(x)开平方exp(x)以e为底的指数l

9、og(x)自然对数 log10 (x)以10为底的对数sign(x)符号函数fix(x)取整表1 常用基本函数表例如，计算，可以采用如下命令：>> y=2*sin(0.3*pi)/(1+sqrt(5)y = 0.5000注意，乘法命令的“*”不可缺少，函数命令用到的小括号不可缺少。在MATLAB中这些函数作用在矩阵或向量上，函数作用在矩阵或向量的每个元素上，得到一个同型矩阵。例如：>> A=pi pi/2 0;B=cos(A)B= -1.000 0.0000 1.0000对于矩阵变量，还有矩阵间的运算。如，矩阵的加法与减法是指对两个大小相等的矩阵，其相应的位置上的数做加

10、法或减法运算，命令如下：>> B=1,2,1;3,3,3;A=1,2,3;4,5,6;>> C=A+Bans = 2 4 4 7 8 9 矩阵乘积，只有当的列数等于的行数或者、两者中有一个是常数时， 才有意义，运算同线性代数中的定义相同。>>B=1,3,4;5,7,8;9,11,12;A=1,2,3;4,5,6;>>A*B;在数学公式中一般把一个矩阵的转置记作，Matlab所用命令为A。>> A=1,2,3;4,5,6;>> Aans = 1 4 2 5 3 6当矩阵为方阵时，矩阵的乘方运算在数学上表述成，若为正整数，则乘

11、方表达式的结果可以将矩阵自乘次得出。例如：>> A=1,2,3;4,2,5;6,2,3;>> A3;Matlab中定义了一种特殊的运算，即所谓的点运算。两个矩阵之间的点运算是它们对应元素的直接运算。例如表示和矩阵的相应元素之间直接进行乘法运算，然后将结果赋给矩阵，即。点乘积运算要求和矩阵的维数相同。>> A=3,-2,3; B=1,2,3; >>D=A.*BD = 3 -4 9>>E=A.2E = 9 4 9还有一些针对向量的函数，可用下列函数找出其元素个数、最大值、最小值等。例如：length(x)求向量x的元素个数、max(x)求

12、向量x的元素的最大值、min(x)求向量x的元素的最小值。7.1.4 Matlab程序设计用Matlab语言编写一些可以调用的程序，称为M文件。M文件有两类：函数M文件和脚本M文件。函数M文件可以输入参数，也可返回输出参数。而脚本M文件没有输入参数与输出参数，可以看作一些命令集成到一起。下面重点介绍函数M文件。一、函数M文件的建立从MATLAB命令窗口的File菜单中选择New菜单项，再选择M-file命令。将需要运行的函数语句编辑到一个文件中，保存函数文件名。在MATLAB命令窗口输入参数和函数名，就可以就会执行此函数。函数文件由function语句引导，其格式为：function 

13、;输出参数=函数名(输入参数)函数体注：一般保存的文件名与函数名相同。 例如， 建立函数，并求时的值。M文件内容如下：function z=f1(x,y)z=2*sin(x)/(1+y);保存文件名为f1.m，在命令窗口调用命令：>>z=f1(3,2)z = 0.0941二、逻辑语句在M文件中，很多时候结构稍复杂，需要用到逻辑语句，常见结构如下。（一） 选择结构 1. 形式一   if  条件n          

14、          语句组  n                 end如果条件成立，则执行语句组。例如，设，求。需要建立M文件f2.m定义函数，再在Matlab命令窗口输入f2(2)、f2(-1)即可。M文件内容如下：function f=f2(x) if x>1 f = x2+1 end if x<=1

15、 f = 2*x end 2. 形式二   if  条件1 n      语句组1 n                 else n      语句组2 nend如果条件1成立，则执行语句组1；否则，则执行语句组2。3. 形式三&

16、#160;  if  条件1 n                    语句组1 n                  elseif  条件2 n&#

17、160;                    语句组2 n                    n        &

18、#160;         elseif  条件m n                    语句组m n             

19、0;   else n                    语句组m+1 n                 end（二） 循环结构 n 实现循环结

20、构的语句：for语句和while语句。  n 1  For语句：  n 格式：   for 循环变量=表达式1:表达式2:表达式3 n                    循环体语句  n     

21、           end n 注：其中表达式1的值为循环变量的初值，表达式2的值为步长，表达式3的值为循环变量的终值。步长为1时，表达式2可以省略。例如，对，分别求的值。编写M文件f3.m如下： for n=1:10 x(n) = sin(n*pi/10); endx在命令窗口调用如下：>>f3可得x=0.3090 0.5878 0.8090 0.9511 1.0000 0.9511 0.8090 0.5878 0.3090 0.0000。 

22、;n 2 While语句  n 格式为：  n          while (条件) n              循环体语句  n         

23、0;end 例如，设银行年利率为11.25%，将10000元钱存入银行，问多长时间会连本带利翻一番？编写M文件f4.m如下：money = 10000 years = 0 ;while money < 20000 years = years +1; money = money * (1+11.25/100); end years money在命令窗口调用如下：>>f4可得years =7，money =2.1091e+004。注意，对所有循环结构的循环体，内部可以包括一个循环结构，称为循环的嵌套。7.1.5 Matlab的图形绘制功能一、二维图形 原理：Matlab

24、作图是通过描点、连线来实现的，故在画一个曲线图形之前，必须先取得该图形上的一系列的点的坐标（即横坐标和纵坐标），然后将该点集的坐标传给Matlab函数画图。命令：plot(X,Y,S)。其中，X、Y是向量，分别表示点集的横坐标和纵坐标。而S表示线型，可以省略，默认为简单连线。常用“r,b,y”分别表示红色、蓝色和黄色，“-”表示连线，“:”短虚线。特殊的，plot(X ,Y1,S1, X,Y2, S2, ,X ,Yn,Sn)是将多条线画在一起例如，可画出一条正弦曲线。>>x=0:0.001:10; % 0到10的1000个点的x座标>>y=sin(x);

25、60;        % 对应的1000个点的y座标 >>plot(x,y);       % 绘图 图2 正弦曲线图>>x=0:0.001:10; y=sin(x); Y=sin(10*x); >>plot(x,y,:'r',x,Y,'b')    % 同时画两个函数，颜色分别是红色虚线和蓝色连线。图3 两个

26、函数的图像绘制一些离散点，比如用“*”来表示，可以用plot(x,y,'*')。除了plot命名，还有fplot可以用来自动的画一个已定义的函数分布图，而无须产生绘图所需要的一组数据，其语法为fplot(fun,a,b)。二、三维图形（一）空间曲线命令：plot3(X,Y,Z,S)其中，“X,Y,Z”是n维向量，分别表示曲线上点集的横坐标、纵坐标、函数值，而S仍表示线型。例如，在区间上画出参数曲线。输入命令: >>t=0:pi/50:10*pi; >>plot3(sin(t),cos(t),t)>>rotate3d %旋转 结果：图4 空间曲

27、线图（二）曲面图命令：surf(x,y,z)其中，“x,y,z”为数据矩阵，分别表示数据点的横坐标、纵坐标、函数值。曲面图一般先得到在x-y平面的矩形网格点上的z轴坐标值。 >>x= -2:0.05:2; % 在x轴上取点>>y= -2:0.05:2; % 在y轴上取点  >>xx,yy=meshgrid(x, y); %形成 xx和yy矩阵 >>zz=xx.*exp(-xx.2-yy.2); % 计算函数值 

28、0;>>surf(xx, yy, zz); % 画出立体曲面图 图5 空间曲面图§7.2 应用实例Matlab提供了许多现成函数解决数学模型问题，这些函数都可以通过Matlab中的帮助命令查看使用说明。下面通过一些实例说明。7.2.1 Matlab在规划问题中的应用 问题一（线性规划）：假设某厂计划生产甲、乙两种产品，现库存主要材料有A类3600公斤，B类2000公斤，C类3000公斤。每件甲产品需用材料A类9公斤，B类4公斤，C类3公斤。每件乙产品，需用材料A类4公斤，B类5公斤，C类10公斤。甲单位产品的利润

29、70元，乙单位产品的利润120元。问如何安排生产，才能使该厂所获的利润最大。 解：设、分别为生产甲、乙产品的件数。为该厂所获总润，建立数学模型： MATLAB程序如下： >>f=-70 -120;%建立目标函数矩阵f   >>A=9 4 4 5;3 10 %建立条件约束对应系数矩阵A   >>b=3600;2000;3000;%建立条件约束对应向量b   >>lb=0&

30、#160;0; ub=; %给出下界，上界为空 >>x,fval,exitflag=linprog(f,A,b,lb,ub)%调用线性规划程序linprog   >>maxf=-fval%求出最大值运行结果：    x =200.0000 240.0000    fval =       -4.2800e+004  

31、0; exitflag =           1      maxf =      4.2800e+004问题二（二次规划）：MATLAB程序如下： >>f=1;-3%建立目标函数一次项系数矩阵>>H=6 -2;-2 8%建立目标函数二次项系数矩阵>>A=2 1

32、;-1 4%建立条件约束矩阵  >>b=2;3 %建立向量  >>X,fval,exitflag=quadprog(H,f,A,b) %使用quadprog二次规划函数来求值    可得结果： ，。问题三（多目标规划）：某钢铁厂准备用5000万用于A、B两个项目的技术改造投资。设、分别表示分配给项目A、B的投资。据专家预估计，投资项目A、B的年收益分别为70%和66%。同时，投资后总的风险损失将随着总投资和单项投资的增加而增加，已知总的风险损失为，问应

33、如何分配资金才能使期望的收益最大，同时使风险损失为最小 。解：数学模型如下通过线性加权构造目标函数，并转化为最小值问题：利用Matlab软件求解，首先编辑目标函数M文件f5.m，M文件如下：function  f=f5(x) %建立目标函数 f=-0.5*(70*x(1)+66*x(2)+0.5*(0.02*x(1)2+0.01*x(2)2+0.04*(x(1)+x(2)2); 调用单目标规划函数fmincon，求最小值问题的函数，命令行如下：>>x0=1000,1000 %构造初始值  &#

34、160; >>A=1 1; b=5000;%建立条件矩阵 >>lb=zeros(2,1); %构造下界 >>x,fval, exitflag=fmincon(f5,x0, A,b,lb,) %使用函数fmincon求解即可7.2.2 Matlab在拟合问题中的应用问题一（多项式拟合）：通过数据拟合Malthus人口指数增长模型中的参数。从17901980年间美国每隔10年的人口记录如下：年份1790180018101820183018401850人口(百万)3.95.3

35、7.29.612.917.123.2年份1860187018801890190019101920人口(百万)31.438.650.262.976.092.0106.5年份193019401950196019701980人口(百万)123.2131.7150.7179.3204.0226.5解： Malthus 模型的基本假设是：人口的增长率为常数记为 ，记时刻的人口为，且初始时刻的人口为，可以得到其解为，左右取对数。Matlab计算如下：>>t=0:19; >>x=3.9 5.3 7.2 9.6

36、60;12.9 17.1 23.2 31.4 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5; >>plot(t,x,'+')  %画散点图 >>y=log(x); p=polyfit(t,y,1)最终计算结果，最终可得。问题二（多元拟合）：电影院调查分别投放电视广告费用和

37、报纸广告费用对每周电影院收入的影响，得到下面的数据，建立收入和费用间关系。每周收入9690959295959494电视广告费用1.521.52.53.32.34.22.5报纸广告费用5242.533.52.53解：编写下面的程序，调用lsqlin()函数来解决该问题。>>x1=1.5,2.0,1.5,2.5,3.3,2.3,4.2,2.5'>>x2=5.0,2.0,4.0,2.5,3.0,3.5,2.5,3.0'>>y=96,90,95,92,95,95,94,94'>>M=x1,x2,ones(size(x1);

38、0;%构造调用lsqlin 函数所需的M 矩阵 >>X=lsqlin(M,y)%调用lsqlin()函数求解二元线性拟合问题 得到结果X = 1.2985    2.3372   83.2116，所以拟合结果为。问题三（非线性拟合）：反应动力学中的Hougen-Waston 模型是非线性模型的一个经典例子，其模型如下：，其中为反应速率。三个决定因素分别为(氢气)，(n-戊烷)，(异戊烷) ，以下是一组实验数据，试建立回归模型，求出未知

39、参数。氢气 470   285   470   470    470    100    100    470   100    100    100    285    

40、;285  -戊烷 300    80    300    80     80      190    80     190   300    300    80 

41、60;   300    190     戊烷 10     10    120   120     10      10      65     65 

42、60;   54    120    120    10     120 反应速率 8.55   3.79  4.82  0.02   2.75   14.39   2.54  4.35   13 &

43、#160;  8.5    0.05   11.32  3.13解：输入 >>data=470 285 470 470 470 100 100 470 100 100 100 285 285; 300 80 300 80 80 190 80 190 300&

44、#160;300 80 300 190; 10 10 120 120 10 10 65 65 54 120 120 10 120'%反应物数据 >>rate=8.55 3.79 4.82 0.02 2.75 14.39 2.54 4.35 13 8.5 0.05 11.32 3.13

45、9; %反应速率数据 >>f=(b,x)( b(1)*x(:,2)-x(:,3)/b(5)./(1+b(2)*x(:,1)+b(3)*x(:,2)+b(4)*x(:,3); %拟合函数 >>b0=1;0.1;0.2;0.1;2; %参数求解的初始值 >>betafit=nlinfit(data,rate,f,b0) %调用nlinfit 函数求解得到结果为：betafit =    1.2526  &

46、#160; 0.0628    0.0400    0.1124    1.1914，所以原模型为：。7.2.3 Matlab在数据统计中的应用问题：一道工序用自动化车床连续加工某种零件，由于刀具损坏等原因，该工序会出现故障，工序出现故障是完全随机的，假定在生产任一零件时出现故障的机会均相同。工作人员通过检查零件来确定工序是否出现故障。现积累有 100次刀具故障记录，故障出现时该刀具完成的零件数如下表：100 次刀具故障记录（已完成的零件数）459 362

47、 624 542 509 584 433 748 815 505612 452 434 982 640 742 565 706 593 680926 653 164 487 734 608 428 1153 593 844527 552 513 781 474 388 824&#

48、160;538 862 659775 859 755 49 697 515 628 954 771 609402 960 885 610 292 837 473 677 358 638699 634 555 570 84 416 606 1062 484 120447 654 564

49、60;339 280 246 687 539 790 581621 724 531 512 577 496 468 499 544 645764 558 378 765 666 763 217 715 310 851试确定刀具的平均寿命，同时判断该刀具出现故障时完成的零件数属于何种分布。解：由问题提供的 100 次刀具故障记录， 我们通过作直方图来近似判断刀

50、具寿命所服从的概率分布，首先在 MATLAB 命令框中录入数据：>> x1=459 362 624 542 509 584 433 748 815 505; >> x2=612 452 434 982 640 742 565 706 593 680; >> x3=926 653 164 487 

51、;734 608 428 1153 593 844; >> x4=527 552 513 781 474 388 824 538 862 659; >> x5=775 859 755 49 697 515 628 954 771 609; >> x6=402 96

52、0 885 610 292 837 473 677 358 638; >> x7=699 634 555 570 84 416 606 1062 484 120; >> x8=447 654 564 339 280 246 687 539 790 581; >

53、> x9=621 724 531 512 577 496 468 499 544 645; >> x10=764 558 378 765 666 763 217 715 310 851; >> x=x1 x2 x3 x4 x5 x6 x7 x8 x9 

54、;x10; 接着用以上输入数据作刀具寿命的频数直方图： >> hist(x,10) 图6 刀具寿命的频数直方图对于近似推断刀具总体寿命的概率分布形式，我们用MATLAB非参数检验命令来加以验证其总体分布类型，以提供初步结论成立的更加可靠的依据： >> normplot(x) 图7 正态分布检验图  可以看到数据基本分布在一条直线上，由normplot命令的定义可初步确定刀具寿命为正态分布。在基本确定所给刀具寿命数据的分布后，就可以估计该分布的某些参数： >>&

55、#160;muhat,sigmahat,muci,sigmaci = normfit(x) 得到muhat = 594，sigmahat = 204.130，muci =553.496，634.504，sigmaci =179.228，237.133。估计出该刀具寿命的均值为594，方差204，均值0.95的置信区间为553.496，634.504，方差的0.95的置信区间为179.228，237.133。结果分析：由上述过程我们可以推断刀具寿命服从正态分布，在总体分布形式已知而方差未知的情形下，我们来检

56、验前面所估计的参数是否可信： >> h,sig,ci = ttest(x,594) h = 0，sig = 1，ci =553.4962，634.5038； 检验结果：布尔变量h = 0，表示不拒绝接受假设，说明提出的假设寿命均值为594是合理的；95%的置信区间为553.4962，634.5038，它完全包括594，且精度很高；sig = 1，远远超过0.5，不能拒绝假设。所以，可以确定刀具平均寿命为594。7.2.4 Matlab

57、在微分方程中的应用 问题一：一只猎犬发现其正东方100米处有一只野兔，野兔以速度向其正北方100米处的洞穴逃跑，猎犬向野兔追去，速度是。 问：1.猎犬的运动轨迹方程； 2.追上兔子的时间、地点。解：模型1：设猎犬的运动轨迹为。图8 猎犬运动大致轨迹由题意得化为微分方程组 当时停止。求解此微分方程组，首先建立方程组的M文件“f6.m”如下：function z=f6(t,y)v=1;z(1)=2*v*(100-y(1)/sqrt(v*t-y(2)2+(100-y(1)2);z(2)=2*v*(v*t-y(2)/sqrt(v*t-y(2)2+(100-y(1)2);z=z(1) z(2)'

58、命令窗口： >>tintvl=0 70; y0=0 0;>>t,y=ode23('f6',tintvl,y0);>>plot(y(:,1),y(:,2),100*ones(size(t),t)得到图像如下：图9 猎犬运动轨迹图一模型2: 由题意得，消去解得故得二阶常微分方程。当时停止, 其中。用M文件f7.m来定义此方程组，内容如下：function z=f7(x,y) z(1)=y(2); z(2)=sqrt(1+y(2).2)./(2*(100-x); z=z(1) z(2)'Matlab命令窗口如下：>>y0=0 0

59、'>>x,y=ode23('f7',0 100,y0);>>plot(x,y(:,1),'o') 得到图像如下：图10 猎犬运动轨迹图二问题二：在城市道路的十字路口，都会设置红绿交通灯。为了让那些正行驶在交叉路口或离交叉路口太近而又无法停下的车辆通过路口，红绿灯转换中间还要亮起一段时间的黄灯。 对于一名驶近交叉路口的驾驶员来说，万万不可处于这样进退两难的境地：要安全停车但又离路口太近；要想在红灯亮之前通过路口又觉得距离太远。那么，黄灯应亮多长时间才最为合理呢？已知城市道路法定速度为，交叉路口的宽度为，典型的车身长度统一定为，一般情

60、况下驾驶员的反应时间为，地面的磨擦系数为。（假设，）分析：对于驶近交叉路口的驾驶员，在他看到黄色信号后要做出决定：是停车还是通过路口。如果他以法定速度（或低于法定速度）行驶，当决定停车时，他必须有足够的停车距离。当驾驶员决定通过路口时，必须有足够的时间让他能完全通过路口。这包括做出停车决定的反应时间以及通过停车所需的最短距离的驾驶时间，能够很快看到黄灯的驾驶员可以利用刹车距离将车停下来。于是，黄灯状态所应持续的时间包括驾驶员的反应时间，他通过交叉路口的时间以及通过刹车距离所需要的时间。解：由于城市道路法定速度为，交叉路口的宽度为，典型的车身长度统一定为。考虑到车通过路口实际上指的是车的尾部必须

61、通过路口。因此，通过路口的时间为。现在我们来计算刹车距离：设为汽车的重量，为磨擦系数，由牛顿力学知，地面对汽车的磨擦力为，其方向与汽车运动的方向相反。汽车在停车过程中，由牛顿第一动力定理有，其中m为汽车质量，为汽车的加速度，是汽车所受的摩擦力。这里加速度是停车距离关于时间的二阶导数，所以行驶距离与时间的关系可由方程来确定，化简同时，我们知道，当时，距离。初速度是距离在0时刻的一阶导数，于是可以给出方程的初始条件。求解方程，Matlab命令行输入：>> x=dsolve('D2x=-ug','x(0)=0,Dx(0)=v0','t')

62、x = t*v0 - (t2*ug)/2即得到停车距离关于时间的解析式。停车时速度为 0，即，可得到汽车刹车所用的时间，从而刹车距离。设黄灯闪烁时间为，则的表达式为。取、，命令行如下：>>I=9;L=4.5;u=0.2;T=1;g=9.8;>>v0=40*1000/3600;%速度转化为米每秒>>A=v0/(2*u*g)+(I+L)/v0+TA = 5.0495也就是说，车速为，黄灯时间为5.05秒。7.2.5 Matlab在模拟问题中的应用问题：在正方形ABCD的四个顶点各有一个人。设在初始时刻时，四人同时出发匀速以沿顺时针走向下一个人。如果他们始终对准下

63、一个人为目标行进，最终结果会如何？并且作出各自的运动轨迹。解：该问题可以通过计算机模拟来实现，需要将时间离散化。设时间间隔为，则时刻表示时间。设第个人时刻的位置坐标为：。 对前面3个人表达式为：其中 ，。对第4个人表达式为：其中 ，。Matlab实现程序f8.m如下：%模拟运动n=240;x=zeros(4,n);y=zeros(4,n);dt=0.05; %时间间隔v=10; %速度x(1,1)=100; y(1,1)=0; %第1个人初始坐标x(2,1)=0; y(2,1)=0; %第2个人初始坐标x(3,1)=0; y(3,1)=100; %第3个人初始坐标x(4,1)=100; y(4,1)=100; %第4个人初始坐标for i=2:n for j=1:3 d=sqrt(x(j+1,i-1)-x(j,i-1)2+(y(j+1,i-1)-y(j,i-1)2); %第j+1个人和第j个人距离 cosx=(x(j+1,i-1)-x(j,i-1)/d; %求cos值 sinx=(y(j+1,i-1)-y(j,i-1)/d; %求sin值 x(j,i)=x(j,i-1)+v*dt*cosx; %求新x坐标