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 = meshgrid(-2:.2:2,-1:.2:2);Z = X.*exp(-X.2
3、- 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=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)+la
5、mda*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,y)单位的百万分之一）键入x,y = meshgrid(-1:.05:1,-1:.05:1);z =exp(-x.2-2*y.2)/104);C,h = contour(x,y,z);h
6、old 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;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./sq
7、rt(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）等高线和梯度线函数的梯度线在实际中有广泛的应用，例如在温度场总热量的流动也是沿着梯度线方向的. 习题1.一块长方形的金属板，四个顶点的坐标是(1,1)，(5,1)，(1,3)，(5,3)在坐标原点处有一个火焰，它使金属板受热假定板上任意一点处的温度与该点到原点的距离成反比在(3,2)处有一个青蛙，问这只青蛙
8、应沿什么方向爬行才能最快到达较凉快的地点？(应沿由热变冷变化最骤烈的方向（即梯度方向）爬行) 第三讲 函数的极值、最值及有关的最优化问题水轮机最优化问题最优化概念反映了人类实践活动中十分普遍的现象，即要在尽可能节省人力、物力和时间前提下，争取获得在可能范围内的最佳效果，因此，最优化问题成为现代数学的一个重要课题，涉及统筹、线性规划一排序不等式等内容。3.1多元函数的偏导数1调用格式一：diff(多元函数,自变量,n)其中，n 为所求偏导数的阶数例 1 已知，求和.解 打开文件编辑窗口，在其中输入下面命令集：pzpx=diff(x2*cos(2*y),x)p2zpypx=diff(pzpx,y)
9、p2zpy2=diff(x2*cos(2*y),y,2)取名为exa9 保存，再在命令窗口中输入命令exa9，程序运行结果如下：pzpx =2*x*cos(2*y)p2zpypx =-4*x*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
10、+5*z)*xuxyz =30*cos(x2-y3+5*z)*y2*xuz3 =-125*cos(x2-y3+5*z)3.2 隐函数的导数在 Matlab 中没有直接求隐函数导数的命令，但可调用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一元函数的极（或
11、最）值例1 求 在中的最小值与最大值主程序为: f=2*exp(-x).*sin(x); fplot(f,0,8); %作图语句 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 函数若要求多元函数
12、在某一区域的极（或最）大值，可转化为求在该区域内的极（或最）小值1非线性无约束情形求极（或最）小值点或极（或最）小值的调用格式是：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(
13、2)3+10*x(1)*x(2)-x(2)2,0,0);fmax=-Fmin;Xmin,fminXmax,fmax取名为exa10保存，再在命令窗口中输入命令exa10，程序运行结果如下：Xmin =1.0016 0.8335fmin =-3.3241Xmax =-1.0000 1.0000fmax =5.00002非线性有约束情形非线性有约束优化问题的数学模型如下： 式中，和是向量，和是矩阵，和为函数，返回量和可以是非线性函数求极（或最）小值点或极（或最）小值的调用格式如下：x,fval=fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon)nonlcon参数计算
14、非线性不等式约束c(x)0,i=1,2,3 建立函数文件fun1打开文件编辑窗口，在其中输入下面命令集：function F=fun1(x) %函数文件必须是function开头F=-x(1)*x(2)*x(3);单击“保存”按钮，自动取名为fun1，再击保存 建立非线性约束函数文件yceqfunction c,ceq1=yceq1(x)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
15、,0;ub=;xmax,fmin=fmincon(fun1,x0,A,b,Aeq,beq,lb,ub,yceq1); %函数要加单引号Vmax=-fmin;xmax,Vmax取名为exa11保存，再在命令窗口中输入命令exa11，程序运行结果如下：xmax =1.00001.00001.0000Vmax =1.0000例6求三元函数满足约束条件的最小值。 建立函数文件fun3打开文件编辑窗口，在其中输入下面命令集：function F=fun3(x) %函数文件必须是function开头F=-x(1)*x(2)*x(3);单击“保存”按钮，自动取名为fun3，再击保存 建立非线性约束函数文件y
16、ceq3把约束条件改写为因为两个约束均为线性，所以符合矩阵不等式的形式，其中，.(2) 编制主程序：打开文件编辑窗口，在其中输入下面命令集：x0=10;10;10; %在此点处开始寻找A=-1 -2 -2; 1 2 2;b=0; 72; xmax,fmin=fmincon(fun3,x0,A,b); %函数要加单引号xmax,fmin取名为zuizhi3保存，再在命令窗口中输入命令zuizhi3，程序运行结果如下： xmin = 24.0000 12.0000 12.0000fmin = -3.4560e+003例 7水轮机最优化问题众所周知，三峡水电站是中国最大的水电站。已知，水是通过输水管
17、从水坝送到发电所，而输水管中水流速度的变化依赖于外界条件。若发电所有三个不同的水电涡轮，每个涡轮都有一个已知的、特定的功率输出函数来提供发电所需的功率。而发电功率又是输送到涡轮的水流速度的函数。将进入发电所的水分成体积不同的三部分，分别到达每个涡轮，因此，我们的目标就是，若给定任意一个总的水流速度，如何分流才能使总的输出功率最大。根据试验证据和伯努利方程，每个涡轮的输出功率可用二次方程式模型描述，其中，表示通过第个涡轮机的水流流速，单位为立方英尺/秒；表示第个涡轮机的输出功率，单位为千瓦；表示进入发电所的总的水流速度，单位为立方英尺/秒。1、如果三个涡轮机一起工作，需要确定通过每个涡轮机的水流
18、速度，使得总输出功率最大。当来流速度一定时，可以用拉格朗日乘子解出总的输出功率在满足约束条件下的最大值，这个最大值一定是来流速度的函数。2、当取何值时问题1中的结果才是有效的？3、当进入发电所的水流速度为2500ft3/s（1立方英尺(ft3)=0.0283立方米(m3)）,该如何分流才能使总输出功率最大？（1）用1的结果计算；（2）直接用Matlab中的非线性有约束优化问题的fmincon函数求解4、现在我们是让三个涡轮一起工作，那有没有可能在某种情况下一个涡轮工作能产生更大的电量呢？做出三个功率输出函数的图像，并讨论当来流速度为1000 ft3/s时应如何分配水流，是让三个涡轮一起工作还是
19、只用一个？（如果只用其中一个，那么该用哪一个）当来流速度为600 ft3/s时又该如何分配呢？5、对于某个来流速度而言，或许选择两个涡轮工作最好。若来流速度为1500 ft3/s，那么应该选择哪两个涡轮一起工作最好？用拉格朗日乘子计算，如何给两个涡轮分配来流才能使输出功率最大，这样分配是否比同时启动三个涡轮更有效？6、如果来流速度为3400 ft3/s，那么该如何分配？解答：1、当来流速度一定时，若同时启动三个涡轮机，那么，根据拉格朗日乘子法计算总输出功率满足约束条件的极大值。设拉格朗日函数为.由此得到微分方程组求解命令如：eq1=sym(0.1277-8.16*10(-5)*Q1)*(170
20、-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(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 = .2966598500340831
21、6291751874573960*QT+20.949784139968189047943649170643 Q3 = 54.232939483450730894493675679770+.36232674392183594637582367643717*QT取四位有效数字得： (3.1)每个分水流速度都是总水流速度的函数。显然，总输出功率也是的函数。2、根据(3.1)式，及（）的取值范围，可以解得， 因此，当同时启动三个涡轮机时，总水流速度必须满足，才能应用(3.1)式计算。3、（1）当总流速度 ft3/s时，在问题2中算出的取值范围之内，应用(3.1)式算出当 ft3/s，ft3/s， ft
22、3/s时总输出功率取得最大值，且最大值为28412千瓦。（2）计算机求解首先建立函数 建立函数文件fun2打开文件编辑窗口，在其中输入下面命令集：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(
23、3)2)*(170-1.6*10(-6)*(x(1)+x(2)+x(3)2);F=-(y1+y2+y3); 建立非线性约束函数文件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);%函数要加单
24、引号Wmax=-fmin;Qmax,Wmax运行得Qmax = 777.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.
25、(-6).*q.*q);y3=(-27.02+0.1380.*x3-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）所示，“:”线表示第一个涡轮机随水流速度的输出功率，“-”线表示第二个涡轮机随水流速度的输出功率，而“”线则表示第三个涡轮机的输出功率。观察上图可知，当来流速为1000ft3/s时，“”线在最上方
26、，如果只用一个涡轮机工作，应选择第三个涡轮机才能使输出功率最大。而当来流速度为600 ft3/s时，“:”线在最上方，即第一个涡轮机输出功率最大。当来流速度为1000 ft3/s时，若只启动第三个涡轮机，根据(3.1)式，及的取值范围解得 ft3/s， ft3/s， ft3/s时算得总的输出功率最大，最大值为千瓦。而如果只用第三个涡轮机工作，则总的输出功率为千瓦。所以应只选择第三个涡轮机单独工作最好。当来流速度为600 ft3/s时，根据的取值范围，不能同时启动三个涡轮机，若只启动一个，则选择第一个涡轮机单独工作输出功率最大，最大值为7292千瓦。5、观察图（3.1）可以看出，当来流速度超过5
27、00 ft3/s后，“-”一直处于最下方，也就是说，在同样的水流速度下，第二个涡轮机的输出功率最小。当来流速度为1500 ft3/s时，若只用两个涡轮机工作，应该选择第一和第三个涡轮机一起工作最好。此时用拉格朗日乘子法计算总输出功率满足约束条件的极大值。设拉格朗日函数由方程此时解得 (3.2)根据和的取值范围，从(3.2)式解得适合此式的的取值范围为。当=1500 ft3/s时，用Matlab求解的解得 ft3/s， ft3/s，总的输出功率最大，最大值为18208 ft3/s。而此时若同时启动三个涡轮机一起工作，根据（3.1）式计算出 ft3/s， ft3/s， ft3/s，总输出功率为16
28、539千瓦。显然，此时只启动两个涡轮机工作更有效。根据(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-1.6.*10.(-6).*q1.*q1);z=y1+y3
29、;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.(-6).*q1.*q1);y3=(-27.02+0.1380.*x3-3.84.*10.(-5).*x3.*x3)
30、.*(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(输出功率);运行得图（3. 2）。虚线表示同时启动第一和第三个涡轮机时总输出功率随总水流速度的变化情况，而实线则表示三个涡轮机同时启动时总输出功率随总水流速度的变化情况。观察图（3. 2），两条曲
31、线大约在处有个交点，说明，当时，同时启动两个涡轮机比同时启动三个能获得更大的总输出功率，而当后，三个涡轮同时启动较好。图（3. 2）总输出功率关于总来流速度的函数图像6、请同学们自己解答。习题1. 求下列函数的极小点 1) ；2) ；3) . 第1），2）题的初始点可任意选取，第3）题的初始点取为.2. 梯子长度问题一楼房的后面是一个很大的花园. 在花园中紧靠着楼房有一个温室,温室伸入花园2m,高3m,温室正上方是楼房的窗台. 清洁工打扫窗台周围,他得用梯子越过温室,一头放在花园中,一头靠在楼房的墙上. 因为温室是不能承受梯子压力的,所以梯子太短是不行的.现清洁工只有一架7m长的梯子,你认为它
32、能达到要求吗? 能满足要求的梯子的最小长度为多少?3. 陈酒出售的最佳时机问题某酒厂有批新酿的好酒,如果现在就出售,可得总收入=50万元(人民币),如果窖藏起来待来日(第年)按陈酒价格出售,第n年末可得总收入(万元),而银行利率为=0.05,试分析这批好酒窖藏多少年后出售可使总收入的现值最大. (假设现有资金万元,将其存入银行,到第年时增值为万元,则称X为的现值.)并填下表：第一种方案：将酒现在出售，所获50万元本金存入银行；第二种方案：将酒窖藏起来，待第年出售。（1） 计算15年内采用两种方案，50万元增值的数目并填入表1，2中；（2） 计算15年内陈酒出售后总收入的现值填入表3中。表1 第
33、一种方案第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年第15年4. 请解答例7中的第6个问题。第四讲 微分方程模型 导弹系统改进、交通管理色灯及自动喷水灭火装置设计实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程. 另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.
34、这就要求我们必须研究微分方程组的解法，既要研究微分方程组的解析解法（精确解），更要研究微分方程组的数值解法（近似解）. 对微分方程组的解析解法，Matlab有专门的函数可以用，本实验将作一定的介绍.然后建模实验研究导弹系统的改进等几个实际问题的微分方程建模及计算机求解.4.1 预备知识微分方程、通解、特解的概念补充一般说来，微分方程就是联系自变量、未知函数以及未知函数的某些导数之间的关系式.如果其中的未知函数只与一个自变量有关，则称为常微分方程；如果未知函数是两个或两个以上自变量的函数，并且在方程中出现偏导数，则称为偏微分方程. 本书所介绍的都是常微分方程，有时就简称微分方程或方程. 微分方程
35、的解中可以包含任意常数，其中任意常数的个数可以多到与方程的阶数相等，也可以不含任意常数. 我们把n阶常微分方程的含有n个独立的任意常数的解 ，称为该方程的通解，如果方程的解不包含任意常数，则称它为特解.4.2相关函数命令及简介1.dsolve( ) 求微分方程的解析解 Matlab语言的符号运算工具箱提供了一个常系数线性微分方程求解的实用函数dsolve( ),该函数允许用字符串的形式描述微分方程及初值、边值条件，最终得出解析解。调用格式指明自变量，其中 可以描述微分方程，又可以描述初始条件或边界条件。D4y表示，D2y(2)=3表示。2.T,Y=solver(odefun,tspan,y0)
36、,求微分方程的数值解.说明:(1)其中solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb之一.(2)odefun是显式常微分方程：（3）在积分区间上，从到，用初始条件求解。（4）要获得问题在其他指定时间点上的解，则令（要求是单调的）。3. inline（）：建立一个内联函数.格式：inline（expr, var1, var2,）,注意括号里的表达式要加引号.4.3计算实验：（自学）直接可以用Matlab求微分方程精确解的例子。 例1假设输入信号为 ,求下面微分方程的通解 （1） syms t;u=exp(-5*t)*cos(2*t
37、+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)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
38、)+92*exp(-5*t)*sin(2*t+1)+10,y(0)=3,Dy(0)=2,D2y(0)=0,D3y(0)=0)例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
39、 -20 25);comet3(x(:,1),x(:,2),x(:,3);function 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(
40、:,3);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分钟作出反应并
41、摧毁导弹.现在要求改进电子导弹系统使能自动计算出敌舰是否在有效打击范围之内.设我舰发射导弹时位置在坐标原点，敌舰在轴正向km处，其行驶速度为km/h,方向与轴夹角为，导弹水平飞行线速度为km/h. 问题的关键是求出导弹击中敌舰的时间.特别地，在导弹系统中设km/h, km/h,h. 现在求的有效范围.求解：设时刻导弹位置为.则 （4.1）易知时刻敌舰位置为，为了保持对准目标，导弹轨迹切线方向应为 （4.2）由（4.1）和（4.2）得下列微分方程， (4.3) （4.4）初始条件对给定的进行计算.当满足， （4.5）则认为已击中目标.如果,则敌舰在打击范围内,可以发射. 有两个极端情形容易算出. 若即敌舰正好背向行驶,即轴正向.那么导弹直线飞行,击中时间 得. 若即迎面驶来,类似地.一般地,应有%M函数missilefun.mfunction dy=missilefun(t,y
