部分的讲函数的等高线

1、第二讲 函数的等高线、梯度线及有关的作图问题鲨鱼袭击目标的前进途径等高线和梯度线有广泛的实际应用，例如在地理学中绘制地貌图，在气象学中绘制气象图等等.本实验通过鲨鱼袭击目标这一例子介绍二元函数的等高线和梯度线的绘制，最后介绍用等高线来做一元隐函数的图形及微分方程的积分曲线.2.1 等高线的绘制二元函数在空间表示的是一张曲面，这个曲面与平面的交线在面上的投影曲线称为函数的一条等高线，我们可以用Matlab作出等高线的图形.等高线的作图命令为“Contour”，最基本格式为Contour二元函数，自变量1，自变量1最小值，自变量1最大值，自变量2，自变量2最小值，自变量2最大值例1 作出在区间上的

2、等高线.解 X,Y = meshgrid(-2:.2:2,-2:.2:3);Z = X.*exp(-X.2-Y.2);C,h = contour(X,Y,Z);set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)colormap cool运行后见图（2.1）.2.2 矢量场图矢量场图（又称速度图）是指由指令quiver实现的.它主要用于描写函数在点的梯度大小和方向。其一般的调用格式为： quiver(X,Y,DZX,DZY)例2 作出函数的等高线和矢量场.解 X,Y = me

3、shgrid(-2:.2:2,-1:.2:2);Z = X.*exp(-X.2 - Y.2);请预览后下载！DX,DY = gradient(Z,.2,.2);% 求二元函数矩阵Z的梯度指令，0.2 为x、y方向上的计算步长. DX,DY是.C,h = contour(X,Y,Z);hold onquiver(X,Y,DX,DY)colormap hsvhold off运行后见图（2.2）. 图（2.1）等高线及其标注图（2.2）等高线和矢量场2.3 梯度线的描绘请预览后下载！设为平面曲线，如果上任意一点处的切线与函数在该店处的梯度位于同一直线上，则称为的梯度线。现在来讨论如何作出函数的梯度线

4、。下面我们一等步长的折线段来近似模拟函数的梯度线。设步长为，从点()出发，沿梯度方向前进得到点，即,再从出发沿梯度线向前进得到点,依次得到一列点，利用"plot"做出此点集的图形，即得梯度线的图形. 例3.作出函数的梯度线.clear allt=cputime;syms x yS= sym(x2 -y2);Sx=diff(S,'x');Sy=diff(S,'y');x0=1;y0=1;lamda=0.01;i=1;sx(1)=x0;sy(1)=y0;for i=2:400 fx=subs(Sx,x,y,sx(i-1),sy(i-1); fy=

5、subs(Sy,x,y,sx(i-1),sy(i-1); sx(i)=sx(i-1)+lamda*fx./sqrt(fx.2+fy.2); sy(i)=sy(i-1)+lamda*fy./sqrt(fx.2+fy.2);请预览后下载！endplot(sx,sy)cputime-t运行后见图（2.3）. 图（2.3）梯度线2.4 鲨鱼袭击目标的前进途径海洋生物学家发现，当鲨鱼在海水中察觉到血液的存在时，就会沿着血液浓度增加得最快的方向前进去寻找目标.根据在海水中世纪测试的结果，如果以流血目标处作为原点在海面上建立直角坐标系，则在海面上点P(x,y)处的血液浓度近似等于（x,y的单位为m,f(x,

6、y)单位的百万分之一）键入x,y = meshgrid(-1:.05:1,-1:.05:1);z =exp(-x.2-2*y.2)/104);C,h = contour(x,y,z);hold on 运行后，作为函数的等高线，得到图（2.4）.由题设条件和梯度的性质可知，鲨鱼袭击目标的前进途径即为的梯度线，下面作出的梯度线，有前面梯度线的绘制可知， 请预览后下载！ 图（2.4）的等高线syms x yS= sym(exp(-x.2-2*y.2)/104);Sx=diff(S,'x');Sy=diff(S,'y');x0=1;y0=1;lamda=0.01;i=1

7、;sx(1)=x0;sy(1)=y0;for i=2:400 fx=subs(Sx,x,y,sx(i-1),sy(i-1); fy=subs(Sy,x,y,sx(i-1),sy(i-1); sx(i)=sx(i-1)+lamda*fx./sqrt(fx.2+fy.2); sy(i)=sy(i-1)+lamda*fy./sqrt(fx.2+fy.2);endplot(sx,sy)hold on得到图（2.5）：请预览后下载！图（2.5）鲨鱼袭击目标的前进途径再运用hold on 命令把等高线和梯度线在同一坐标系中显示，得图（6）. 图（2.6）等高线和梯度线函数的梯度线在实际中有广泛的应用，例如

8、在温度场总热量的流动也是沿着梯度线方向的. 习题1.一块长方形的金属板，四个顶点的坐标是(1,1)，(5,1)，(1,3)，(5,3)在坐标原点处有一个火焰，它使金属板受热假定板上任意一点处的温度与该点到原点的距离成反比在(3,2)处有一个青蛙，问这只青蛙应沿什么方向爬行才能最快到达较凉快的地点？(应沿由热变冷变化最骤烈的方向（即梯度方向）爬行) 请预览后下载！第三讲 函数的极值、最值及有关的最优化问题水轮机最优化问题最优化概念反映了人类实践活动中十分普遍的现象，即要在尽可能节省人力、物力和时间前提下，争取获得在可能范围内的最佳效果，因此，最优化问题成为现代数学的一个重要课题，涉及统筹、线性规

9、划一排序不等式等内容。3.1多元函数的偏导数1调用格式一：diff('多元函数','自变量',n)其中，n 为所求偏导数的阶数例 1 已知，求和.解 打开文件编辑窗口，在其中输入下面命令集：pzpx=diff('x2*cos(2*y)','x')p2zpypx=diff(pzpx,'y')p2zpy2=diff('x2*cos(2*y)','y',2)取名为exa9 保存，再在命令窗口中输入命令exa9，程序运行结果如下：pzpx =2*x*cos(2*y)p2zpypx =-4*x

10、*sin(2*y)p2zpy2 =-4*x2*cos(2*y)即,.2调用格式二：syms x y z diff(f,自变量,n)例2 已知，求解 在命令行中依次输入：请预览后下载！syms x y zu=sin(x2-y3+5*z);ux=diff(u,x);uxy=diff(ux,y);uxyz=diff(uxy,z);uz3=diff(u,z,3);ux,uxyz,uz3运行结果如下：ux =2*cos(x2-y3+5*z)*xuxyz =30*cos(x2-y3+5*z)*y2*xuz3 =-125*cos(x2-y3+5*z)3.2 隐函数的导数在 Matlab 中没有直接求隐函数导

11、数的命令，但可调用Maple 中求隐函数导数的命令，调用格式如下：maple('implicitdiff(f(u,x,y,z,，)=0,u,x)')例3 求由多元方程x2 + y2 + z2 = xyz所确定的隐函数解 在命令行中输入：pzpx=maple('implicitdiff(x2+y2+z2-x*y*z=0,z,x)')运行结果是：pzpx =(2*x-y*z)/(-2*z+x*y)即xy zx yz3.3一元函数的极（或最）值例1 求 在中的最小值与最大值主程序为: f='2*exp(-x).*sin(x)' fplot(f,0,8)

12、; %作图语句请预览后下载！ xmin,ymin=fminbnd (f, 0,8) f1='-2*exp(-x).*sin(x)' xmax,zmin=fminbnd (f1, 0,8)ymax=-zmin运行结果： xmin = 3.9270 ymin = -0.0279 xmax = 0.7854 ymax = 0.6448 3.4 多元函数的极（或最）值在 Matlab 中同样有求多元函数的极（或最）小值的函数，但由于多元函数的形式比较复杂，不同情况用到不同的Matlab 函数若要求多元函数在某一区域的极（或最）大值，可转化为求在该区域内的极（或最）小值1非线性无约束情形

13、求极（或最）小值点或极（或最）小值的调用格式是：x,fval=fminsearch(f,x0) %fminsearch是不能设定约束范围的f是被最小化的目标函数名，x0 是求解的初始值向量例 求二元函数的最值点和最值解 打开文件编辑窗口，在其中输入下面命令集：%必须对自变量进行转化x=x(1)，y=x(2)Xmin,fmin=fminsearch('2*x(1)3+4*x(1)*x(2)3-10*x(1)*x(2)+x(2)2',0,0);Xmax,Fmin=fminsearch('-2*x(1)3-4*x(1)*x(2)3+10*x(1)*x(2)-x(2)2'

14、;,0,0);fmax=-Fmin;Xmin,fminXmax,fmax取名为exa10保存，再在命令窗口中输入命令exa10，程序运行结果如下：Xmin =1.0016 0.8335fmin =-3.3241Xmax =-1.0000 1.0000请预览后下载！fmax =5.00002非线性有约束情形非线性有约束优化问题的数学模型如下： 式中，和是向量，和是矩阵，和为函数，返回量和可以是非线性函数求极（或最）小值点或极（或最）小值的调用格式如下：x,fval=fmincon('fun',x0,A,b,Aeq,beq,lb,ub,nonlcon)nonlcon参数计算非线性不

15、等式约束c(x)<=0 和非线性等式约束ceq(x)=0例 5 求表面积为的体积最大的长方体体积解 设长方体的长、宽、高分别为x1、x2、x3，则f(x)=-x(1)*x(2)*x(3)，S.t x(1)*x(2)+x(2)*x(3)+x(3)*x(1)-3=0，x(i)>0,i=1,2,3 建立函数文件fun1打开文件编辑窗口，在其中输入下面命令集：function F=fun1(x) %函数文件必须是function开头F=-x(1)*x(2)*x(3);单击“保存”按钮，自动取名为fun1，再击保存 建立非线性约束函数文件yceqfunction c,ceq1=yceq1(x

16、)c=x(1)*x(2)+x(2)*x(3)+x(3)*x(1)-3;ceq1=;保存方法同上，自动取名为yceq1，再击保存 编制主程序：打开文件编辑窗口，在其中输入下面命令集：x0=3;3;3; %给长宽高一个初值A=;b=;Aeq=;beq=;请预览后下载！lb=0,0,0;ub=;xmax,fmin=fmincon('fun1',x0,A,b,Aeq,beq,lb,ub,'yceq1'); %函数要加单引号Vmax=-fmin;xmax,Vmax取名为exa11保存，再在命令窗口中输入命令exa11，程序运行结果如下：xmax =1.00001.0000

17、1.0000Vmax =1.0000例6求三元函数满足约束条件的最小值。 建立函数文件fun3打开文件编辑窗口，在其中输入下面命令集：function F=fun3(x) %函数文件必须是function开头F=-x(1)*x(2)*x(3);单击“保存”按钮，自动取名为fun3，再击保存 建立非线性约束函数文件yceq3把约束条件改写为因为两个约束均为线性，所以符合矩阵不等式的形式，其中，.(2) 编制主程序：打开文件编辑窗口，在其中输入下面命令集：x0=10;10;10; %在此点处开始寻找A=-1 -2 -2; 1 2 2;b=0; 72; xmax,fmin=fmincon('

18、fun3',x0,A,b); %函数要加单引号xmax,fmin请预览后下载！取名为zuizhi3保存，再在命令窗口中输入命令zuizhi3，程序运行结果如下： xmin = 24.0000 12.0000 12.0000fmin = -3.4560e+003例 7水轮机最优化问题众所周知，三峡水电站是中国最大的水电站。已知，水是通过输水管从水坝送到发电所，而输水管中水流速度的变化依赖于外界条件。若发电所有三个不同的水电涡轮，每个涡轮都有一个已知的、特定的功率输出函数来提供发电所需的功率。而发电功率又是输送到涡轮的水流速度的函数。将进入发电所的水分成体积不同的三部分，分别到达每个涡轮，

19、因此，我们的目标就是，若给定任意一个总的水流速度，如何分流才能使总的输出功率最大。根据试验证据和伯努利方程，每个涡轮的输出功率可用二次方程式模型描述，其中，表示通过第个涡轮机的水流流速，单位为立方英尺/秒；表示第个涡轮机的输出功率，单位为千瓦；表示进入发电所的总的水流速度，单位为立方英尺/秒。1、如果三个涡轮机一起工作，需要确定通过每个涡轮机的水流速度，使得总输出功率最大。当来流速度一定时，可以用拉格朗日乘子解出总的输出功率在满足约束条件下的最大值，这个最大值一定是来流速度请预览后下载！的函数。2、当取何值时问题1中的结果才是有效的？3、当进入发电所的水流速度为2500ft3/s（1立方英尺(

20、ft3)=0.0283立方米(m3)）,该如何分流才能使总输出功率最大？（1）用1的结果计算；（2）直接用Matlab中的非线性有约束优化问题的fmincon函数求解4、现在我们是让三个涡轮一起工作，那有没有可能在某种情况下一个涡轮工作能产生更大的电量呢？做出三个功率输出函数的图像，并讨论当来流速度为1000 ft3/s时应如何分配水流，是让三个涡轮一起工作还是只用一个？（如果只用其中一个，那么该用哪一个）当来流速度为600 ft3/s时又该如何分配呢？5、对于某个来流速度而言，或许选择两个涡轮工作最好。若来流速度为1500 ft3/s，那么应该选择哪两个涡轮一起工作最好？用拉格朗日乘子计算，

21、如何给两个涡轮分配来流才能使输出功率最大，这样分配是否比同时启动三个涡轮更有效？6、如果来流速度为3400 ft3/s，那么该如何分配？解答：1、当来流速度一定时，若同时启动三个涡轮机，那么，根据拉格朗日乘子法计算总输出功率满足约束条件的极大值。设拉格朗日函数为.由此得到微分方程组求解命令如：请预览后下载！eq1=sym('(0.1277-8.16*10(-5)*Q1)*(170-1.6*10(-6)*QT2)+Q4=0');eq2=sym('(0.1358-9.38*10(-5)*Q2)*(170-1.6*10(-6)*QT2)+Q4=0');eq3=sym(

22、'(0.1380-7.68*(10(-5)*Q3)*(170-1.6*(10(-6)*QT2)+Q4=0');eq4=sym('Q1+Q2+Q3-QT=0');Q1,Q2,Q3,Q4=solve(eq1,eq2,eq3,eq4)运行得：Q1 = .34101340604408089070665757782322*QT-75.182723623418919942437324850413 Q2 = .29665985003408316291751874573960*QT+20.949784139968189047943649170643 Q3 = 54.232939

23、483450730894493675679770+.36232674392183594637582367643717*QT取四位有效数字得： (3.1)每个分水流速度都是总水流速度的函数。显然，总输出功率也是的函数。2、根据(3.1)式，及（）的取值范围，可以解得， 请预览后下载！因此，当同时启动三个涡轮机时，总水流速度必须满足，才能应用(3.1)式计算。3、（1）当总流速度 ft3/s时，在问题2中算出的取值范围之内，应用(3.1)式算出当 ft3/s，ft3/s， ft3/s时总输出功率取得最大值，且最大值为28412千瓦。（2）计算机求解首先建立函数 建立函数文件fun2打开文件编辑窗口

24、，在其中输入下面命令集：function F=fun2(x) %函数文件必须是function开头y1=(-18.89+0.1277*x(1)-4.08*10(-5)*x(1)2)*(170-1.6*10(-6)*(x(1)+x(2)+x(3)2);y2=(-24.51+0.1358*x(2)-4.69*10(-5)*x(2)2)*(170-1.6*10(-6)*(x(1)+x(2)+x(3)2);y3=(-27.02+0.1380*x(3)-3.84*(10(-5)*x(3)2)*(170-1.6*10(-6)*(x(1)+x(2)+x(3)2);F=-(y1+y2+y3); 建立非线性约束

25、函数文件yceq2function c,ceq2=yceq2(x)c=x(1)+x(2)+x(3)-2500;ceq2=; 编制主程序：打开文件编辑窗口，在其中输入下面命令集：x0=20;30;10; %从该点处开始解决问题 A=;b=;Aeq=;beq=; %如果没有不等式存在,则令A=,b= lb=250,250,250;ub=1110,1110,1225;Qmax,fmin=fmincon('fun2',x0,A,b,Aeq,beq,lb,ub,'yceq2');%函数要加单引号Wmax=-fmin;Qmax,Wmax请预览后下载！运行得Qmax = 77

26、7.3508 762.5994 960.0498Wmax = 2.8412e+0044、做出每个涡轮机单独工作时其输出功率关于来流速度的函数图像，程序如下:x1=250:0.01:1110;x2=250:0.01:1110;x3=250:0.01:1250;q=2500;y1=(-18.89+0.1277.*x1-4.08.*10.(-5).*x1.*x1).*(170-1.6.*10.(-6).*q.*q);y2=(-24.51+0.1358.*x2-4.69.*10.(-5).*x2.*x2).*(170-1.6.*10.(-6).*q.*q);y3=(-27.02+0.1380.*x3-

27、4.08.*10.(-5).*x3.*x3).*(170-1.6.*10.(-6).*q.*q);plot(x1,y1,':') hold onplot(x2,y2,'-')plot(x3,y3,'-')xlabel('流速');ylabel('输出功率');运行得图（3.1）请预览后下载！图（3.1）三个涡轮机单独工作时其输出功率关于来流速度的函数图像如图（3.1）所示，“:”线表示第一个涡轮机随水流速度的输出功率，“-”线表示第二个涡轮机随水流速度的输出功率，而“”线则表示第三个涡轮机的输出功率。观察上图可知，

28、当来流速为1000ft3/s时，“”线在最上方，如果只用一个涡轮机工作，应选择第三个涡轮机才能使输出功率最大。而当来流速度为600 ft3/s时，“:”线在最上方，即第一个涡轮机输出功率最大。当来流速度为1000 ft3/s时，若只启动第三个涡轮机，根据(3.1)式，及的取值范围解得 ft3/s， ft3/s， ft3/s时算得总的输出功率最大，最大值为千瓦。而如果只用第三个涡轮机工作，则总的输出功率为千瓦。所以应只选择第三个涡轮机单独工作最好。当来流速度为600 ft3/s时，根据的取值范围，不能同时启动三个涡轮机，若只启动一个，则选择第一个涡轮机单独工作输出功率最大，最大值为7292千瓦。

29、5、观察图（3.1）可以看出，当来流速度超过500 ft3/s后，“-”一直处于最下方，也就是说，在同样的水流速度下，第二个涡轮机的输出功率最小。当来流速度为1500 ft3/s时，若只用两个涡轮机工作，应该选择第一和第三个涡轮机一起工作最好。此时用拉格朗日乘子法计算总输出功率请预览后下载！满足约束条件的极大值。设拉格朗日函数由方程此时解得 (3.2)根据和的取值范围，从(3.2)式解得适合此式的的取值范围为。当=1500 ft3/s时，用Matlab求解的解得 ft3/s， ft3/s，总的输出功率最大，最大值为18208 ft3/s。而此时若同时启动三个涡轮机一起工作，根据（3.1）式计算

30、出 ft3/s， ft3/s， ft3/s，总输出功率为16539千瓦。显然，此时只启动两个涡轮机工作更有效。根据(3.2)式和(3.1)式，绘出总输出功率关于总来流速度的函数图像。程序如下：l=0.01;%步长q0=953;qe=3636;q1=q0:l:qe;x1=-65.0253+0.4848.*q1;x3=65.0253+0.5152.*q1;y1=(-18.89+0.1277.*x1-4.08.*10.(-5).*x1.*x1).*(170-1.6.*10.(-6).*q1.*q1);y3=(-27.02+0.1380.*x3-3.84.*10.(-5).*x3.*x3).*(170

31、-1.6.*10.(-6).*q1.*q1);z=y1+y3;请预览后下载！plot(q1,z,':','LineWidth',2)hold onq1=q0:l:qe;x1=-75.1827+0.341.*q1;x2=-20.9498+00.2967.*q1;x3=54.233+0.3623.*q1;y1=(-18.89+0.1277.*x1-4.08.*10.(-5).*x1.*x1).*(170-1.6.*10.(-6).*q1.*q1);y2=(-24.51+0.1358.*x2-4.69.*10.(-5).*x2.*x2).*(170-1.6.*10.(

32、-6).*q1.*q1);y3=(-27.02+0.1380.*x3-3.84.*10.(-5).*x3.*x3).*(170-1.6.*10.(-6).*q1.*q1);y=y1+y2+y3;h=z-y;aa=find(h=min(abs(h);%找交点所对应的迭代次数q2=l*aa+q0 %计算交点所对应的流速y(aa) %计算交点所对应的输出功率text(q2,y(aa),'*(2.1038e+003,2.3898e+004)')plot(q1,y,'-','LineWidth',2)xlabel('流速');ylabel(

33、'输出功率');运行得图（3. 2）。虚线表示同时启动第一和第三个涡轮机时总输出功率随总水流速度的变化情况，而实线则表示三个涡轮机同时启动时总输出功率随总水流速度的变化情况。观察图（3. 2），两条曲线大约在处有个交点，说明，当时，同时启动两个涡轮机比同时启动三个能获得更大的总输出功率，而当后，三个涡轮同时启动较好。请预览后下载！图（3. 2）总输出功率关于总来流速度的函数图像6、请同学们自己解答。习题1. 求下列函数的极小点 1) ；2) ；3) . 第1），2）题的初始点可任意选取，第3）题的初始点取为.2. 梯子长度问题一楼房的后面是一个很大的花园. 在花园中紧靠着楼房有

34、一个温室,温室伸入花园2m,高3m,温室正上方是楼房的窗台. 清洁工打扫窗台周围,他得用梯子越过温室,一头放在花园中,一头靠在楼房的墙上. 因为温室是不能承受梯子压力的,所以梯子太短是不行的.现清洁工只有一架7m长的梯子,你认为它能达到要求吗? 能满足要求的梯子的最小长度为多少?请预览后下载！3. 陈酒出售的最佳时机问题某酒厂有批新酿的好酒,如果现在就出售,可得总收入=50万元(人民币),如果窖藏起来待来日(第年)按陈酒价格出售,第n年末可得总收入(万元),而银行利率为=0.05,试分析这批好酒窖藏多少年后出售可使总收入的现值最大. (假设现有资金万元,将其存入银行,到第年时增值为万元,则称X

35、为的现值.)并填下表：第一种方案：将酒现在出售，所获50万元本金存入银行；第二种方案：将酒窖藏起来，待第年出售。（1） 计算15年内采用两种方案，50万元增值的数目并填入表1，2中；（2） 计算15年内陈酒出售后总收入的现值填入表3中。表1 第一种方案第1年第2年第3年第4年第5年第6年第7年第8年第9年第10年第11年第12年第13年第14年第15年表2 第二种方案第1年第2年第3年第4年第5年第6年第7年第8年第9年第10年第11年第12年第13年第14年第15年表3 陈酒出售后的现值第1年第2年第3年第4年第5年请预览后下载！第6年第7年第8年第9年第10年第11年第12年第13年第14

36、年第15年4. 请解答例7中的第6个问题。第四讲 微分方程模型 导弹系统改进、交通管理色灯及自动喷水灭火装置设计实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程. 另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程组的解法，既要研究微分方程组的解析解法（精确解），更要研究微分方程组的数值解法（近似解）.请预览后下载！ 对微分方程组的解析解法，Matlab有专门的函数可以用，本实验将作一定的介绍.然后建模实验研究导弹系统的改进等几个实际问题的微分方程建模及计算机求解.4.1 预备知识微分方程、通解、特解的概念补充一般说来，微

37、分方程就是联系自变量、未知函数以及未知函数的某些导数之间的关系式.如果其中的未知函数只与一个自变量有关，则称为常微分方程；如果未知函数是两个或两个以上自变量的函数，并且在方程中出现偏导数，则称为偏微分方程. 本书所介绍的都是常微分方程，有时就简称微分方程或方程. 微分方程的解中可以包含任意常数，其中任意常数的个数可以多到与方程的阶数相等，也可以不含任意常数. 我们把n阶常微分方程的含有n个独立的任意常数的解 ，称为该方程的通解，如果方程的解不包含任意常数，则称它为特解.4.2相关函数命令及简介1.dsolve( ) 求微分方程的解析解 Matlab语言的符号运算工具箱提供了一个常系数线性微分方

38、程求解的实用函数dsolve( ),该函数允许用字符串的形式描述微分方程及初值、边值条件，最终得出解析解。调用格式指明自变量，其中 可以描述微分方程，又可以描述初始条件或边界条件。D4y表示，D2y(2)=3表示。2.T,Y=solver(odefun,tspan,y0),求微分方程的数值解.说明:(1)其中solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb之一.请预览后下载！(2)odefun是显式常微分方程：（3）在积分区间上，从到，用初始条件求解。（4）要获得问题在其他指定时间点上的解，则令（要求是单调的）。3. inline

39、（）：建立一个内联函数.格式：inline（expr, var1, var2,）,注意括号里的表达式要加引号.4.3计算实验：（自学）直接可以用Matlab求微分方程精确解的例子。 例1假设输入信号为 ,求下面微分方程的通解 （1） syms t;u=exp(-5*t)*cos(2*t+1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*usyms t y;y=dsolve('D4y+10*D3y+35*D2y+50*Dy+24*y=','87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10')

40、latex(y)例2 假设 方程（1）满足初始条件 求特解>>syms t;u=exp(-5*t)*cos(2*t+1)+5;uu=5*diff(u,t,2)+4*diff(u,t)+2*usyms t y;y=dsolve('D4y+10*D3y+35*D2y+50*Dy+24*y=','87*exp(-5*t)*cos(2*t+1)+92*exp(-5*t)*sin(2*t+1)+10','y(0)=3','Dy(0)=2','D2y(0)=0','D3y(0)=0')请预览后下载！

41、例3 求解下面线性微分方程组 x,y=dsolve('D2x+2*Dx=x+2*y-exp(-t)','Dy=4*x+3*y+4*exp(-t)')例4Lorenz模型其中设. 令初值为.试求解该微分方程组.求解：编程一：function lo()t_final=100;x0=0;0;1e-10;t,x=ode45(lorenzeq,0,t_final,x0);plot(t,x),figure;plot3(x(:,1),x(:,2),x(:,3);axis(10 42 -20 20 -20 25);comet3(x(:,1),x(:,2),x(:,3);func

42、tion xdot=lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3)-10*x(2)+10*x(3) -x(1)*x(2)+28*x(2)-x(3);编程二：f1=inline('-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3)','t','x');请预览后下载！t_final=100;x0=0;0;1e-10;t,x=ode45(f1,0,t_final,x0);plot(t,x),figure;plot3(x(:,1),x(:,2),x(:,3)

43、;axis(10 42 -20 20 -20 25);comet3(x(:,1),x(:,2),x(:,3);例1 求解微分方程例2求微分方程在初始条件下的特解，并画出解函数的图形。例3求微分方程组在初始条件下的特解，并画出解函数的图形。例4求解微分方程初值问题的数值解，求解范围为0,0.5.例5求解描述振荡器的经典的VerderPol微分方程4.4建模实验：例1 导弹系统的改进海军方面要求改进现有的舰对舰系统.目前的电子系统能迅速测出敌舰的种类、位置以及敌舰行驶的速度和方向，且导弹自动制导系统能保证在发射后的任意时刻都能对准目标.根据情报，这种敌舰能在我军舰发射导弹后T分钟作出反应并摧毁导弹

44、.现在要求改进电子导弹系统使能自动计算出敌舰是否在有效打击范围之内.设我舰发射导弹时位置在坐标原点，敌舰在轴正向km处，其行驶速度为km/h,方向与轴夹角为，导弹水平飞行线速度为km/h. 问题的关键是求出导弹击中敌舰的时间.特别地，在导弹系统中设km/h, 请预览后下载！ km/h,h. 现在求的有效范围.求解：设时刻导弹位置为.则 （4.1）易知时刻敌舰位置为，为了保持对准目标，导弹轨迹切线方向应为 （4.2）由（4.1）和（4.2）得下列微分方程， (4.3) （4.4）初始条件对给定的进行计算.当满足， （4.5）则认为已击中目标.如果,则敌舰在打击范围内,可以发射. 有两个极端情形容

45、易算出. 若即敌舰正好背向行驶,即轴正向.那么导弹直线飞行,击中时间 得. 若即迎面驶来,类似地.一般地,应有请预览后下载！%M函数missilefun.mfunction dy=missilefun(t,y,a,b,d,theta)dydx=(a*t*sin(theta)-y(2) +1e-8)/.(abs(d+a*t*cos(theta)-y(1)+1e-8);dy(1)=b/(1+dydx2)0.5;dy(2)=b/(1+dydx(-2)0.5;dy=dy(:);在指令窗口执行下列脚本文件得-5.7410. %M脚本missilefunªclear;close;a=90;b=4

46、50;d=50;theta=pi/2;t,y=ode45(missilefun,0,0.1,0 0,a,b,d,theta);plot(y(:,1),y(:,2);max(y(:,1)-d-a*t*cos(theta)由于在内，（5）式不成立，所以敌舰不在有效打击范围，应等近一些再发射.例2. 交通管理色灯中黄灯应亮多长时间？让我们来考虑这样一个问题：红绿灯在亮红灯之前黄灯应亮多长时间？在交通管理中，定期的亮一段时间黄灯是为了让那些正行驶在交叉路口上或距交叉路口太近以致无法停下的车辆通过路口。这样，红绿灯应保持足够长时间的黄灯，便那些无法停下的驾驶员有机会在黄灯亮的时候通过路口。对于一位驶近交

47、叉路口的驾驶员来说，万万不可处于这样的进退两难的境地：要安全停车，离路口太近，而要在红灯亮之前通过路口又显得太远了。关于这个时间的计算，有一个粗糙的“经验方法”：对法定迫近速度的每个10mi/h亮1s黄灯。我们来看理论计算是否证实这个经验方法。驶近交叉路口的驾驶员，在看到黄色信号后要作出决定：是停车还是通过路口。如果他以法定速度（或低于法定速度）行驶，当决定停车时，他必须有足够的停车距离。当决定通过路口时，他必须有足够的时间使他能够完全通过路口，这包括作出停车决定的时间（反应时间）以及通过停车所需的最短距离的驾驶时间。能够很快看到黄灯的驾驶员可以利用刹车距离将车停下。于是，黄灯状态应持续的时间

48、包括驾驶员的反应时间、他通过交叉路口的时间以及停车所需的时间（刹车距离）。有了这么多的时间，驾驶员就能够在刹车距离内安全停车。请预览后下载！如果法定速度为，交叉路口的宽度为,典型的车身长度为，那么通过路口的时间为。（注意车的尾部必须通过路口，这样，路口的实际长度就是.)现在，我们来计算刹车距离。注意到实际的刹车和停车过程是相当复杂的，驾驶员首先开加速器，然后不同程度地用劲踩住刹车踏板（也许用了气动刹车）使车子减速，直到停止。这些过程的大部分是很难用精确地模型来描述的，因而我们将回避它们。通过在刹车过程引入一个抵抗摩擦力，我们将假定刹车效果可用模型来反映。设为汽车重量，为摩擦系数。根据定义，对汽

49、车的制动力为,其方向与运动方向相反. 汽车在停车过程中，行驶的距离可通过解下面的微分方程求得，这个方程反映的是在常力作用下的直线运动： (4.6) 其中是重力加速度。赋予的正确条件是时，以及. 于是，刹车距离就是直到时汽车驶过的距离。 求解：在的条件下对（4.6）积分，得 (4.7)因此，当时，速度为零。在条件下对（4.7）积分，得请预览后下载！ (4.8)时，的值是 （4.9）上述求解的过程可以用计算机求解：（huangdeng1.m）syms t x f g v0 y;x=dsolve('D2x=','-f*g','Dx(0)=v0',

50、9;x(0)=0')x =-1/2*f*g*t2+v0*t注意，对必须用到的单位要谨慎：对我们所涉及的距离，交通工程师是用来度量的（ft为英尺的简写，1英尺=0.3048米），而速度则通常用（mi/h=0.44704m/s）来度量。重力加速度= 32.16在作任何计算之前，读者应将转换成. 我们将计算一下黄灯状态其中是驾驶员的反应时间。于是 假设 另外，我们将接受公路工程师提出的具有代表性的(参看练习1）当=30，40以及50时，黄灯时间如表1所示。表中同时也给出了经验法的值。请预览后下载！黄灯周期与法定速度的关系编写程序如下（huangdeng2.m）v0=1:0.1:50;f=0.

51、2;g=9.8/0.3048;I=30;L=15;T=1;A=(v0*0.44704/0.3048)./(2*f*g)+(I+L)./(v0*0.44704/0.3048)+T;text(30,5.46,'*')plot(v0,A)xlabel('v_0(mi/h)');ylabel('A(s)');关于的图像描绘出来，如图 4.1所示。图4.1 黄灯周期与法定速度的关系 表 1（mi/h) A 经验法30 5.46s 3s40 6.35s 4s50 7.34s 5s 我们注意到，经验法的结果一律比我们预测的黄灯状态短些。这使人想起，许多交叉路口

52、的设计很可能使车辆在红绿灯转为红灯时正处于交叉路口上。请预览后下载！ 即使给了充分的停车时间，仍有许多汽车司机企图加速想抢在红灯亮之前冲过十字路口。这当中的大多数驾驶员不知道（有些人甚至不予注意）什么时候绿灯转红。有一种“倒数”型的红绿灯可以部分地解决这个问题，在黄灯亮的最后几秒钟内，黄灯上显示出一串倒着数的数字，它们准确地向驾驶员警告红绿灯何时将变色。这种系统几年前就在休斯顿试用了，它成功地降低了事故发生率。参看练习6. 习题1 (广告效应)某公司生产一种耐用消费品,市场占有率为时开始做广告,一段时间的市场跟踪调查后,该公司发现:单位时间内购买人口百分比的相对增长率与当时还没有买的百分比成正比,且估得此比例系数为0.5. (1)建立该问题的数学模型,分别求其解析解和数值解,并作比较. (2)厂家问:要做多少时间广告,可使市场占有率达到