《经济数学基础》课程教学资源（PPT讲稿）实验2 缉私艇追击走私船

实验2 辑私枫誼击是私好 2,】实龄目的 通过本次实验，目的在于提高学生根据实际问题建立微分方程模型的能 力，学会求微分方程解析解与数值解的方法，并会做简单的计算机仿真。 2,2实肱网题 海上边防缉私艇发现距C公里处有一走私船正以匀速沿直线行驶，缉私艇 立即以最大速度追赶，在雷达的引导下，缉私艇的方向始终指向走私船。问 缉私艇何时追赶上走私船？并求出缉私艇追赶的路线

实验2 缉私艇追击走私船 2.1 实验目的 通过本次实验，目的在于提高学生根据实际问题建立微分方程模型的能 力，学会求微分方程解析解与数值解的方法，并会做简单的计算机仿真。 2.2 实验问题 海上边防缉私艇发现距c公里处有一走私船正以匀速a沿直线行驶,缉私艇 立即以最大速度b追赶,在雷达的引导下，缉私艇的方向始终指向走私船。问 缉私艇何时追赶上走私船?并求出缉私艇追赶的路线。 x y o c

2,3鸡机型 走私船初始位置在点(0,0)， 行驶方向为y轴正方向， 缉私艇的初始位置在点(CG0), 缉私艇行驶的历程为5。 在时刻t: 走私船的位置到达点R(0,a) 缉私艇到达点D(x,) y-at ds tga b dx x-0 dt y(c)=0,y'(c)=0 r=a/b 12y dt dt dt ds -a 2 dx dx d ds dx b dx

2.3 建立模型 o c x y 走私船初始位置在点(0,0) ， 行驶方向为 y轴正方向 ， 缉私艇的初始位置在点 ( c,0) ， 缉私艇行驶的历程为 s 。 在时刻 t ： 缉私艇到达点 R ( 0 ,at ) D ( x , y ) 走私船的位置到达点 − 0 − = = xy at tg dx dy  b dt ds = dx dt a dx d y x = − 2 2 2 1 1   = = − + dx dy dx b ds ds dt dx dt  =  =   = + ( ) 0, ( ) 0 1 2 2 2 y c y c dx dy r dx d y x r = a / b

2.4棋型求解 (1)求解析解 dy dp dx" r+ d (c)=0,y'(C)=0 dp dx r=a/b p(c)=0 p+i+p 图 y(c)=0

2.4 模型求解      =  =       = + ( ) 0, ( ) 0 1 2 2 2 y c y c dx dy r dx d y x r = a /b (1) 求解析解 令： p dx dy = ， dx dp dx d y = 2 2 ，      = = + ( ) 0 1 2 p c x dx r p dp r c x p p       + + = 2 1 r x c p p       − + = − 2 1      =                −      = ( ) 0 2 1 y c x c c x dx dy r r

图 +周食 (c)=0 cr bc 当x=0时， a(1-r2 (b2-a2 c=3千米，a=0.4千米/秒， 分别取b=0.6,0.8,1.2千米/秒时， 3.5 缉私艇追赶路线的图形。 3 2 追赶时间分别为： 1.5 =9,5,2.8125(分钟) 0.5 Matlab 0 0 0.5 11.522.533.5

     =                −      = ( ) 0 2 1 y c x c c x dx dy r r = 1 b a 1) r , 2 1 1 1 1 1 1 1 2 r cr c x c r x r c y r r − +               −  −      + = + − 当 x = 0 时， 2 1 r cr y − = ， (1 ) ( ) 2 2 2 b a bc a r cr a y t − = − = = 0 0.5 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 3.5 4 c=3千米，a=0.4千米/秒， 分别取b=0.6,0.8,1.2千米/秒时， 缉私艇追赶路线的图形。 追赶时间分别为： t=9，5，2.8125（分钟） Matlab

2) >1 目9 (c)=0 2- 当X→0时，V)十0，缉私艇不可能追赶上走私船。 当X→0时，y→+©0，缉私艇不可能追赶上走私船

2) = 1 b a r      =                −      = ( ) 0 2 1 y c x c c x dx dy r r 1 1 1 1 1 2 2 1 1 − −               −  +      + = + − r cr x c c r x r c y r r 当 x →0 时， y → + ，缉私艇不可能追赶上走私船。 3） r = 1 ， ，         − − = c x c c x c y ln 2 2 1 2 2 当 x →0 时， y → + ，缉私艇不可能追赶上走私船

(2)用MATLAB软件求解析解 MATLAB软件5.3以上版本提供的解常微分方程解析解的指令是 dsolve,完整的调用格式是 dsolve ('eqn1','eqn2',.) 其中‘eqn1,egn2',.,是输入宗量，包括三部分：微分方程、初始 条件、指定变量，若不指定变量，则默认小写字母t为独立变量。 微分方程的书写格式规定：当y是因变量时，用“Dy”表示的n阶导数。 例求微分方程三X十》的通解。 dsolve ('Dy=x+x*y 图 Ans=-1+exp(1/2*x^2)*C1 Matlab y(c)=0 dso1ve('Dy=1/2*(x/e)^r-(c/x)r)','y(e)=0','x') ans=1/2*ep(-r*(1og(c)-1og()))*c^r*(1/c)^r/(r+1)*x+ 1/2*ep(r*(1og(c)-1og(x))/(-1+r)*x- 1/2*c*(-(1/c)^r*c^x+C^r*(1/C)^r*x+r+1)/(r^2-1)

（2）用MATLAB软件求解析解 MATLAB软件5.3以上版本提供的解常微分方程解析解的指令是 dsolve，完整的调用格式是: dsolve('eqn1','eqn2', .) 其中‘eqn1’ , ‘eqn2’ , .是输入宗量，包括三部分：微分方程、初始 条件、指定变量，若不指定变量，则默认小写字母t为独立变量。 微分方程的书写格式规定：当y是因变量时，用“Dny”表示y的n阶导数。 例 求微分方程 y  = x + xy 的通解。 dsolve('Dy=x+x*y','x') Ans=-1+exp(1/2*x^2)*C1      =                −      = ( ) 0 2 1 y c x c c x dx dy r r dsolve('Dy=1/2*((x/c)^r-(c/x)^r)','y(c)=0','x') ans=1/2*exp(-r*(log(c)-log(x)))*c^r*(1/c)^r/(r+1)*x+ 1/2*exp(r*(log(c)-log(x)))/(-1+r)*x- 1/2*c*(-(1/c)^r*c^r+c^r*(1/c)^r*r+r+1)/(r^2-1) Matlab

(3)用MATLAB软件求数值解 c=3,a=0.4,b=0.8,r=a/b=0.5 程序ZX.m y(c)=0 function y=zx(t,y) y=0.5*((t/3)^0.5-(3/t)^0.5) 执行下面的命令：ode23('zx',3,0,0005,0) matlab 若想看图中“0”点的坐标可执行下面的命令： [t,y]=ode23('zx',3,0.0005,0) plot(t,y) 此时缉私艇的位置坐标是(0.00050000000000,1.96013657712118) 执行下面的命令：ode45('',3,0.0005,0) 若想看图中“0”点的坐标可执行下面的命令： [t,y]=ode45('zx',3,0.0005,0) plot(t,y) 此时缉私艇的位置坐标是(0.0005,1.9675)

（3）用MATLAB软件求数值解      =                −      = ( ) 0 2 1 y c x c c x dx dy r r c=3，a=0.4，b=0.8， r = a/b = 0.5 程序zx.m function y=zx(t,y) y=0.5*((t/3)^0.5-(3/t)^0.5) 执行下面的命令：ode23('zx',3,0.0005,0) matlab 若想看图中“o”点的坐标可执行下面的命令： [t,y]=ode23('zx',3,0.0005,0) plot(t,y) 此时缉私艇的位置坐标是（0.00050000000000，1.96013657712118） 执行下面的命令: ode45('zx',3,0.0005,0) 若想看图中“o”点的坐标可执行下面的命令： [t,y]=ode45('zx',3,0.0005,0) plot(t,y) 此时缉私艇的位置坐标是（0.0005，1.9675 ）

(4)用MATLAB软件防真法 当建立动态系统的微分方程模型很困难时， 我们可以用计算机仿真法对系统进行分析研究。 所调计算机仿真就是利用计算机对实际动态系 统的结构和行为进行编程、模拟和计算，以此 来预测系统的行为效果。 走私船初始位在点(0,0)，方向为y轴正方向，缉私艇的初始位在点(C0), t三t:走私船的位置(0，t)，缉私艇的位置：(X,》】 追赶方向可用方向余弦表示为： 0-x coSQ= (0-x)广+(a-y) atk-y sm0三 (0-x)+(at-y) Xk+1-X=△xk≈b△t cosa, 时间步长为△，则在时刻1十△1时： yk+1-yk=Ayk≈b△tsin a

(4) 用MATLAB软件防真法 当建立动态系统的微分方程模型很困难时， 我们可以用计算机仿真法对系统进行分析研究。 所谓计算机仿真就是利用计算机对实际动态系 统的结构和行为进行编程、模拟和计算，以此 来预测系统的行为效果。 o c x y 方向为y轴正方向， ： 走私船的位置： 缉私艇的位置： 走私船初始位在点(0,0)， 缉私艇的初始位在点(c,0)， (0, ) t = t k at k ， ( , ) k k x y 追赶方向可用方向余弦表示为： 2 2 (0 ) ( ) 0 cos k k k k k x at y x − + − −  = 2 2 (0 ) ( ) sin k k k k k k x at y at y − + − −  = 时间步长为 t ，则在时刻 t t k +  时： + − =    cos , k 1 k k k x x x b t  k k k k y +1 − y = y  btsin 

仿真算法 第一步：设置时间步长△t,速度a,b及初始位置X0=C,%=0 第二步：计算动点缉私艇D在时刻1+1二1k+△f时的坐标（比+y) x+1=Xk+b△ ,Vk+1=y+b△ at-y Vx元+(at&-ye) x+(atk-y) compare 计算走私船R在时刻i+l=i+△f时的坐标代， 81=0，+1=aU+△) 第三步：计算缉私艇与走私船这两个动点之间的距离： =V(xx+1-元x4)2+(y+1-立+1)2 根据事先给定的距离，判断缉私艇是否已经追上了走私船，从而判断 退出循环还是让时间产生一个步长，返回到第二步继续进入下一次循环； 第四步：当从上述循环退出后，由点列(X》)和)可分别绘 制成两条曲线即为缉私艇和走私船走过的轨迹曲线

仿真算法： 第一步：设置时间步长 t , 速度a, b及初始位置 x = c 0 ， 0 y0 = 第二步：计算动点缉私艇D在时刻 t t t k+1 = k +  时的坐标 ( , ) k+1 k+1 x y 2 2 1 ( ) k k k k k k x at y x x x b t + − − + = +  2 2 1 ( ) k k k k k k x at y at y y y b t + − − ， + = +  compare 计算走私船R在时刻 t t t k+1 = k +  时的坐标 ) ~ , ~( k+1 k+1 x y 0 ~ xk+1 = ( ) ~ 1 y a t t ， k+ = k +  第三步：计算缉私艇与走私船这两个动点之间的距离： 2 1 1 2 1 1 ) ~ ) ( ~ ( k = k+ − k+ + k+ − k+ d x x y y 根据事先给定的距离，判断缉私艇是否已经追上了走私船，从而判断 退出循环还是让时间产生一个步长，返回到第二步继续进入下一次循环； 第四步：当从上述循环退出后，由点列 ( , ) k+1 k+1 x y 和 ) ~ , ~( k+1 k+1 x y 可分别绘 制成两条曲线即为缉私艇和走私船走过的轨迹曲线

取C=3千米，=0.4千米/分钟，b=0.8千米/分钟，r三☑/b=0.5, 程序zhuixiantu.m c=3,a=0.4/60;b=0.8/60 jstxb=[];jstyb=[];zscxb=[];zscyb= d=0.01;dt=2;=0; jstx-c:jsty-0;zscx-0;zscy=0; while (sqrt((jstx-zscx)2+(jsty-zscy)2)>d) t=t+dt; jstx-=jstx-b*dt*jstx/sqrt(jstx 2+(a*t-jsty)2); jstxb=[jstxb,jstx]; jsty=jsty+b*dt*(a*t-jsty)/sqrt(jstx 2+(a*t-jsty)2); jstyb=[jstyb,jsty]; zscy=a*t; zscyb=[zscyb,zscy]; end zscxb=zeros(length(zscyb)); plot (jstxb,jstyb,zscxb,zscyb,'*') matlab 历时：296秒，艇： (0.00004135527497,1.97150212967435) 船： (0,1.97333333333333)

取c=3千米，a=0.4千米/分钟，b=0.8千米/分钟， r = a/b = 0.5 ， 程序zhuixiantu.m c=3; a=0.4/60; b=0.8/60; jstxb=[];jstyb=[];zscxb=[];zscyb=[]; d=0.01;dt=2;t=0; jstx=c;jsty=0;zscx=0;zscy=0; while (sqrt((jstx-zscx)^2+(jsty-zscy)^2)>d) t=t+dt; jstx=jstx-b*dt*jstx/sqrt(jstx^2+(a*t-jsty)^2); jstxb=[jstxb,jstx]; jsty=jsty+b*dt*(a*t-jsty)/sqrt(jstx^2+(a*t-jsty)^2); jstyb=[jstyb,jsty]; zscy=a*t; zscyb=[zscyb,zscy]; end zscxb=zeros(length(zscyb)); plot(jstxb,jstyb,zscxb,zscyb,'*') matlab 历时：296秒，艇：（0.00004135527497，1.97150212967435） 船：（0，1.97333333333333）