缉私艇追走私船模型实验课件.pptx

1、2 2、了解微分方程数值解思、了解微分方程数值解思想，掌握两种简单的微分想，掌握两种简单的微分方程数值解方法。方程数值解方法。3、学会根据实际问题建立、学会根据实际问题建立简单微分方程数学模型。简单微分方程数学模型。实验目的实验目的1 1、学会用、学会用MATLABMATLAB软件求解软件求解微分方程的初值问题微分方程的初值问题. .4、了解计算机数据仿真、了解计算机数据仿真、数据模拟的基本方法。数据模拟的基本方法。 实验十二实验十二 缉私艇追赶走私船模型实验缉私艇追赶走私船模型实验1.1.用用MATLABMATLAB软件求微分方程解析解、数值解软件求微分方程解析解、数值解理论理论. . 2.

2、2.编程求解微分方程数值解编程求解微分方程数值解. . 3.3.微分方程模型实验微分方程模型实验: :缉私艇追赶走私船缉私艇追赶走私船. . 实验内容实验内容1. 1. 微分方程的解微分方程的解一阶常微分方程的初值问题一阶常微分方程的初值问题: : 4 用用matlab软件求软件求微分方程的解析解微分方程的解析解ldsolve(equation,condition) %求方程求方程equation在在初始条件初始条件condition下的解，自变量默认为下的解，自变量默认为tldsolve(equation) %求方程求方程equation的通解的通解l%一阶导数用一阶导数用Dy表示，二阶导数

3、用表示，二阶导数用D2y表示，表示， A=dsolve(Dy=5) B=dsolve(Dy=x,x) %求方程的通解，指定自变求方程的通解，指定自变 量为量为x C=dsolve(D2y=1+Dy) D=dsolve(D2y=1+Dy,y(0)=1,Dy(0)=0) x,y=dsolve(Dx=y+x, Dy=2*x,x(0)=0, y(0)=1) %解微分方程组解微分方程组例例1 (1)dsolve(D2y=exp(x),x)dsolve(D2y=exp(x)dsolve(D2y=exp(x),y(0)=1,Dy(0)=2,x) 微分方程数值解的计算微分方程数值解的计算 欧拉方法、欧拉两步法

4、、改进的欧拉方法欧拉方法、欧拉两步法、改进的欧拉方法1)1)欧拉方法：欧拉方法：是一种求解给定初值的常微分方是一种求解给定初值的常微分方程的一阶数值解方法。程的一阶数值解方法。对于方程对于方程设设y=y(x)是它的解，则在解是它的解，则在解y=y(x)的曲线上每一点的曲线上每一点处的切线的斜率等于处的切线的斜率等于f(x,y)在该点处的值。在该点处的值。所以，过初始点所以，过初始点A0(x0,y0)，做切线做切线)(,(0000 xxyxfyy 与直线与直线x=x1x=x1交点的纵坐标交点的纵坐标y1y1近似代替近似代替y=y(x)在在x=x1x=x1处的函数值处的函数值y(x1)，即，即)(

5、,()(0100011xxyxfyyxy 然后再过点然后再过点A1(x1,y1)，以以f(x1,y1)为斜率做切为斜率做切线，与直线线，与直线x=x2交点的纵坐标交点的纵坐标y2y2近似代替近似代替y=y(x)在在x=x2处的函数值处的函数值y(x2)，即，即A0 x0 x1x2A1A2)(,()(1211122xxyxfyyxy 这样依次类推得到上述微分方程的解在这样依次类推得到上述微分方程的解在x=xn+1x=xn+1处处函数值的近似计算公式，函数值的近似计算公式，)(,()(111nnnnnnnxxyxfyyxy 待求的曲线为蓝色，它的折线近似为红色待求的曲线为蓝色，它的折线近似为红色A

6、0 x0 x1x2x3A1A2A3A4x4特别地，如果节点之间是等分的，则特别地，如果节点之间是等分的，则),()(11nnnnnyxhfyyxy 这就是著名的这就是著名的欧拉欧拉(Euler)(Euler)公式。公式。事实上，上述公式是用函数事实上，上述公式是用函数y(x)y(x)在在xnxn处的向前差处的向前差商代替函数在这点处的导数，从而将微分方程离商代替函数在这点处的导数，从而将微分方程离散化。同样，如果用函数散化。同样，如果用函数y(x)y(x)在在xn+1xn+1处的向后差处的向后差商代替导数值，就得商代替导数值，就得后退的欧拉公式后退的欧拉公式),()(1111 nnnnnyxh

7、fyyxy 注意：这两个公式虽然具有相同的精度，但也是有区别的，注意：这两个公式虽然具有相同的精度，但也是有区别的，前者是计算前者是计算yn+1yn+1的显示公式，后者是隐式公式。的显示公式，后者是隐式公式。A0 x0 x1x2x3A1A2A3A4x4将上两式算术平均即得将上两式算术平均即得梯形公式：梯形公式：),(),(2)(1111 nnnnnnnyxfyxfhyyxy梯形公式消去了欧拉公式和后退的欧拉公式误差梯形公式消去了欧拉公式和后退的欧拉公式误差的主要部分，因此比它们具有更高的精度。的主要部分，因此比它们具有更高的精度。这是一个隐式公式，可采用迭代法求解。这是一个隐式公式，可采用迭代

8、法求解。通常先用欧拉公式提供一个通常先用欧拉公式提供一个ynyn附近的初始迭代值，附近的初始迭代值，然后再用梯形公式做迭代。然后再用梯形公式做迭代。),()0(1nnnnyxhfyy ),(),(2)(11)1(1knnnnnknyxfyxfhyy 2 2）欧拉两步法欧拉两步法用函数用函数y(x)y(x)在在xnxn处的中心差商代替微分方程中的处的中心差商代替微分方程中的导数，导数，)(,(2)()(11nnnnxyxfhxyxy 这就是这就是欧拉两步法公式。欧拉两步法公式。)(,(2)()(11nnnnxyxhfxyxy 它比欧拉公式、后退的欧拉公式具有更高的精度。它比欧拉公式、后退的欧拉公

9、式具有更高的精度。3 3）改进的欧拉公式改进的欧拉公式 欧拉公式虽然简单，但是误差较大，并且欧拉公式虽然简单，但是误差较大，并且随着迭代次数的增加，误差越来越大。梯形公随着迭代次数的增加，误差越来越大。梯形公式虽然提高了精度，但因其是隐式的，所以计式虽然提高了精度，但因其是隐式的，所以计算麻烦，工作量较大。通常将这两种方法结合，算麻烦，工作量较大。通常将这两种方法结合，得到得到改进的欧拉公式，改进的欧拉公式，具体如下：具体如下：),(1nnnnyxhfyy 步骤步骤1 1 先由欧拉公式得到先由欧拉公式得到y(xn+1)y(xn+1)的一个初始近似的一个初始近似值，称为值，称为预测值预测值步骤步

10、骤2 2 对预测值再用梯形公式校正一次得对预测值再用梯形公式校正一次得yn+1yn+1，称，称为为校正值校正值),(),(2111 nnnnnnyxfyxfhyy整个过程（称整个过程（称预测预测- -校正系统校正系统）可以表示成）可以表示成),(,(),(21nnnnnnnnyxhfyhxfyxfhyy 这就是这就是改进的欧拉公式。改进的欧拉公式。当然也可以将欧拉两步法与梯形方法结合。当然也可以将欧拉两步法与梯形方法结合。再解释：再解释：欧拉公式用函数欧拉公式用函数y(x)y(x)在区间在区间xn,xn+1xn,xn+1上的平均斜上的平均斜率代替了曲线在点率代替了曲线在点xnxn的斜率。的斜率

11、。上的平均斜率。上的平均斜率。在区间在区间为函数为函数称称,)()()(11 nnnnxxxyhxyxy得到启发：得到启发：只要对曲线在区间只要对曲线在区间xn,xn+1xn,xn+1上的平均上的平均斜率提供一种算法，就可以得到一种计算斜率提供一种算法，就可以得到一种计算yn+1yn+1的公的公式。如果设法在区间内多预测几个点的斜率值，然式。如果设法在区间内多预测几个点的斜率值，然后将这些点处后将这些点处斜率值的加权平均值斜率值的加权平均值作为曲线在区间作为曲线在区间上的平均值，由此构造出由上的平均值，由此构造出由ynyn计算计算yn+1yn+1的精度更高的精度更高的计算公式，这就是的计算公式

12、，这就是龙格龙格- -库塔方法库塔方法的基本思想。的基本思想。欧拉两步法用区间欧拉两步法用区间xn,xn+1xn,xn+1和区间和区间xn-1,xnxn-1,xn上的上的平均斜率的平均斜率的算术平均值算术平均值代替了曲线在代替了曲线在xnxn的斜率。的斜率。编程计算微分方程数值解编程计算微分方程数值解例例2 2 已知初值问题已知初值问题333(cos )2,22|()xxyyxexxxxye （1）用）用MATLAB软件求解析解，算出解析解在结软件求解析解，算出解析解在结点点x=pi:0.1:2*pi处的纵坐标；处的纵坐标；（2）取步长）取步长0.1，用欧拉公式求数值解；，用欧拉公式求数值解；

13、（3）取步长）取步长0.1，用改进的欧拉公式求数值解；，用改进的欧拉公式求数值解；（4）比较两种数值解和解析解，并画出数值解和）比较两种数值解和解析解，并画出数值解和解析解曲线。解析解曲线。解：（解：（1 1）直接在命令窗执行命令）直接在命令窗执行命令dsolve(Dy=3/x*y+x3*(exp(x)+cos(x)-2*x, y(pi)=(exp(pi)+2/pi)*pi3,x)ans= x3*exp(x)+x3*sin(x)+2*x2clc;clf;szy_eu=;y=(exp(pi)+2/pi)*pi3;for x=pi:0.1:2*piy=y+0.1*(3/x*y+x3*(exp(x)

14、+cos(x)-2*x);szy_eu=szy_eu,y;endszy_eu t=pi:0.1:2*pi;f=t.3.*exp(t)+t.3.*sin(t)+2*t.2;plot(t,f,b*-) hold onplot(t,szy_eu,r-,linewidth,2)(2) (2) 取步长取步长0.10.1，用欧拉公式求数值解，用欧拉公式求数值解x=pi:0.1:2*pi;h=0.1; y=(exp(pi)+2/pi)*pi3;szy_imeu=y; for i=1:length(x)-1 jzj=y+h*(3/x(i)*y+x(i)3*(exp(x(i)+cos(x(i)-2*x(i);

15、yy=y+h/2*(3/x(i)*y+x(i)3*(exp(x(i)+cos(x(i)- 2*x(i)+3/x(i+1)*jzj+x(i+1)3*(exp(x(i+1)+cos(x(i+1)-2*x(i+1); szy_imeu=szy_imeu,yy; y=yy; end szy_imeu plot(x,szy_imeu,r-,linewidth,2) f=x.3.*exp(x)+x.3.*sin(x)+2*x.2; hold on plot(x,f,b*-)(3) (3) 取步长取步长0.10.1，用改进的欧拉公式求数值解，用改进的欧拉公式求数值解 2. 2. 用用MATLABMATLAB

16、软件求微分方程数值解软件求微分方程数值解 MATLAB MATLAB软件求常微分方程数值解的指令有软件求常微分方程数值解的指令有 ode23 ode23，ode23tode23t和和ode45ode45，分别表示二三阶龙，分别表示二三阶龙格格库塔法，二三阶龙格库塔法，二三阶龙格库塔梯形法和四五库塔梯形法和四五阶龙格阶龙格库塔法。还有其他命令库塔法。还有其他命令ode113ode113，ode15iode15i，ode15sode15s，ode23sode23s等，读者可通过等，读者可通过helphelp命令查询。命令查询。 xb,yb=ode23(fname,xs,xe,sv) 返回结点处的横

17、坐标返回结点处的横坐标xbxb和纵坐标和纵坐标ybyb。ode23(fname,xs,xe,sv) 其中用单引号引起的是存储微分方程的函数文件名，其中用单引号引起的是存储微分方程的函数文件名，xsxs表示自变量的初始值，表示自变量的初始值，xexe表示自变量的终止值，表示自变量的终止值，svsv表示迭代初始向量值。该命令执行后自动画出微分方表示迭代初始向量值。该命令执行后自动画出微分方程数值解曲线。程数值解曲线。最简单的使用格式为最简单的使用格式为注意：注意：（1 1）命令）命令odeode求解的是形如求解的是形如y y=f(t,y)=f(t,y)的微分方的微分方程，我们称它为一阶导数可解出的

18、微分方程。程，我们称它为一阶导数可解出的微分方程。而对于一阶导数解不出的形如而对于一阶导数解不出的形如 f(t,y,yf(t,y,y)=0)=0 的的微分方程，可以用命令微分方程，可以用命令ode15iode15i求解，有兴趣的求解，有兴趣的读者可以用读者可以用helphelp命令获得。命令获得。（2 2）求解微分问题时必须为其指定初值，即上）求解微分问题时必须为其指定初值，即上面的面的svsv。 （3 3）odeode只能直接求解一阶微分方程，高阶微分方只能直接求解一阶微分方程，高阶微分方程必须等价地转化成一阶微分方程组，才能用程必须等价地转化成一阶微分方程组，才能用odeode命令求解。比

19、如命令求解。比如 y y=f(t,y,y=f(t,y,y,y,y) )， 设设y1=y, y2= yy1=y, y2= y, y3= y, y3= y ， 那么如上的三阶微分方程就可以表示为如下的三那么如上的三阶微分方程就可以表示为如下的三个一阶方程组了。个一阶方程组了。 y1=y2y1=y2 y2=y3y2=y3 y3=f(t,y1,y2,y3)y3=f(t,y1,y2,y3)例例3 3使用格式使用格式: fc.m: fc.mfunction f=fc(x,y)function f=fc(x,y)f=2f=2* *y+x+2;y+x+2;在在MATLABMATLAB命令窗中执行命令命令窗中执

20、行命令ode23(ode23(fcfc,0,1,1),0,1,1)例例4 4 求解微分方程求解微分方程建立建立m m文件文件 erjie.merjie.mfunction f=erjie(t,y)f=y(2); 1000*(1-y(1)2)*y(2)-y(1);在在MATLABMATLAB命令窗中执行命令命令窗中执行命令t,y=ode15s(erjie,0,3000,2,0)plot(t,y(:,1),-) 0)0(, 2)0(0)1(1000222xxxdtdxxdtxd4. 4. 模型实验模型实验: : 缉私艇追赶走私船缉私艇追赶走私船 海上边防缉私艇发现距海上边防缉私艇发现距c公里处有一

21、走私船正以匀公里处有一走私船正以匀速速a沿直线行驶沿直线行驶, ,缉私艇立即以最大速度缉私艇立即以最大速度b追赶追赶, ,在在雷雷达的引导下，缉私艇的方向始终指向走私船。问达的引导下，缉私艇的方向始终指向走私船。问缉缉私艇何时追赶上走私船私艇何时追赶上走私船? ?并求出缉私艇追赶的路线并求出缉私艇追赶的路线方程。方程。xyco解解 (1) (1) 建立模型建立模型走私船初始位在点走私船初始位在点(0,0)，方向为，方向为y轴正方向，轴正方向，缉私艇的初始位在点缉私艇的初始位在点(c,0)，在时刻，在时刻t，走私船的，走私船的位置到达点位置到达点R(0,at)，缉私艇的位置到达，缉私艇的位置到达

22、D(x,y), 0(atR ),(yxD可得缉私艇追击走私船过程的数学模型：可得缉私艇追击走私船过程的数学模型：(2) (2) 模型求解模型求解求解析解求解析解令：令：1)1 abr表示追赶时间。表示追赶时间。表示走私船跑过的距离。表示走私船跑过的距离。a=0.4;c=3;for b=0.6:0.3:1.2r=a/b;x=0:0.05:c;y=c/2*(1/(r+1)*(x/c).(1+r)-1/(1-r)*(x/c).(1-r)+c*r/(1-r2);plot(x,y,r-,linewidth,2)pausehold onendaxis(0 4 0 3.5)gtext(b=0.6)gtext

23、(b=0.9)gtext(b=1.2)取取c=3c=3千米，千米，a=0.4a=0.4千米千米/ /分钟，分别取分钟，分别取b=0.6,0.9,b=0.6,0.9,1.21.2千米千米/ /分钟时，编程绘制缉私艇追赶路线图形分钟时，编程绘制缉私艇追赶路线图形. .1)2 abr1)3 r定义定义m文件函数：文件函数：function y=zx(t,y)y=0.5*(t/3)0.5-(3/t)0.5); ode23(zx,3,0.0005,0) (3)(3)用用MATLABMATLAB软件求数值解软件求数值解设设c=3c=3千米，千米，a=0.4a=0.4千米千米/ /分钟，分钟，b=0.8b=

24、0.8千米千米/ /分钟，分钟，r=a/b=0.5,r=a/b=0.5,用二三阶龙格用二三阶龙格- -库塔算法计算数值解库塔算法计算数值解: :(4) (4) 动态仿真动态仿真 当建立动态系统的微分方程模型很困难时，当建立动态系统的微分方程模型很困难时，可以用计算机仿真法对系统进行分析研究。所可以用计算机仿真法对系统进行分析研究。所谓计算机仿真就是利用计算机对实际动态系谓计算机仿真就是利用计算机对实际动态系统的结构和行为进行编程、模拟和计算，以此统的结构和行为进行编程、模拟和计算，以此来预测系统的行为效果。来预测系统的行为效果。走私船初始位在点走私船初始位在点(0,0)(0,0)，方向为，方向

25、为y y轴正方向，轴正方向，缉私艇的初始位在点缉私艇的初始位在点( (c c,0),0)，走私船的位置到达点走私船的位置到达点), 0(katR，在在时时刻刻kt缉私艇的位置到达点缉私艇的位置到达点 ),(kkyxD追赶方向可用方向余弦表示为追赶方向可用方向余弦表示为： ), 0(katRk ),(kkyxD22)()0(0coskkkkkyatxx 22)()0(sinkkkkkkyatxyat ,1时时ttttkk 时间步长时间步长缉私艇的位置：缉私艇的位置： ).,(11 kkyx则则 ,cos1kkkktbxxx kkkktbyyy sin1 ), 0(1 kat), 0(katRk

26、),(kkyxD).,(11 kkyx仿真算法：仿真算法： 第二步：第二步： 计算动点缉私艇计算动点缉私艇D D在时刻在时刻 tttkk 1时的坐标时的坐标 221)(kkkkkkyatxxtbxx 221)(kkkkkkyatxyattbyy 计算走私船计算走私船R在时刻在时刻 tttkk 1时的坐标时的坐标 ),(11 kkyx01 kx)(1ttaykk 第一步：设置时间步长第一步：设置时间步长, 速度速度a, b及初始位置及初始位置t 0, 0,0 kycx仿真算法步骤：仿真算法步骤：第三步第三步：计算缉私艇与走私船这两个动点之间的距离：计算缉私艇与走私船这两个动点之间的距离： 211

27、211)()( kkkkkyyxxd 根据事先给定的距离，判断缉私艇是否已经追根据事先给定的距离，判断缉私艇是否已经追上了走私船，从而判断退出循环还是让时间产生一上了走私船，从而判断退出循环还是让时间产生一个步长，返回到第二步继续进入下一次循环；个步长，返回到第二步继续进入下一次循环； 第四步第四步：当从上述循环退出后，由点列当从上述循环退出后，由点列 ),(11 kkyx和和),(11 kkyx可分别绘制成两条曲线可分别绘制成两条曲线即为缉私艇和走私船走过的轨迹曲线。即为缉私艇和走私船走过的轨迹曲线。 取取c=3c=3千米，千米，a a=0.4=0.4千米千米/ /分钟，分钟，b b=0.8

28、=0.8千米千米/ /分钟，分钟，r=a/b=0.5r=a/b=0.5显示船与艇行进路线程序显示船与艇行进路线程序clc;clear;clf;c=3; a=0.4/60; b=0.8/60;d=0.01;dt=2;t=0;jstx=c;jsty=0;zscx=0;zscy=0;while (sqrt(jstx-zscx)2+(jsty-zscy)2)d) plot(jstx,jsty,ro,markersize,15); hold on plot(zscx,zscy,b*,markersize,15); pause(0.01) t=t+dt; jstx=jstx-b*dt*jstx/sqrt(

29、jstx2+(a*t-jsty)2); jsty=jsty+b*dt*(a*t-jsty)/sqrt(jstx2+(a*t-jsty)2); zscy=a*t;endfprintf(追上所需时间追上所需时间t=%.0fn,t);fprintf(缉私艇的位置缉私艇的位置=(%.5f,%.5f)n,jstx,jsty);fprintf(走私船的位置走私船的位置=(%.5f,%.5f)n,zscx,zscy);结果分析结果分析计算机仿真法计算的结果依赖于时计算机仿真法计算的结果依赖于时间迭代步长的选取和程序终止条件间迭代步长的选取和程序终止条件的设定，修改终止条件的设定和减的设定，修改终止条件的设定

30、和减小时间迭代步长可以提高计算精度，小时间迭代步长可以提高计算精度，减小误差。减小误差。 上机任务上机任务李继成（书）李继成（书） Page 175 练习练习2 第第1、2、3题题 Page 179 练习练习3 第第3、4、2题题 Page 186 练习练习4 第第4、1、3题题 学号末位除以学号末位除以3 3，余数为，余数为0 0做第一列，余数为做第一列，余数为1 1做做第二列，余数为第二列，余数为2 2做第三列。做第三列。 1 如图，有一半径为4a的大圆，里面有有一个小圆，半径为a，现在，大圆固定而小圆在大圆内相切滚动，设M起点的坐标为(4a,0),确定小圆上一点M的轨迹。M2改变人、狗的速度，曲线是否变化？改变人、狗的速度，曲线是否变化？改变人跑步的曲线，比如椭圆、摆线、星形线、心形线改变人跑步的曲线，比如椭圆、摆线、星形线、心形线等，狗的曲线会是什么样？等，狗的曲线会是什么样？