12-缉私艇追走私船模型实验

1、Mathematics Laboratory 阮小娥博士 Experiments in Mathematics 赵小艳 数学实验数学实验 办公地址：理科楼办公地址：理科楼218 2 2、了解微分方程数值解思、了解微分方程数值解思 想，掌握两种简单的微分想，掌握两种简单的微分 方程数值解方法。方程数值解方法。 3、学会根据实际问题建立、学会根据实际问题建立 简单微分方程数学模型。简单微分方程数学模型。 实验目的实验目的 1 1、学会用、学会用MATLABMATLAB软件求解软件求解 微分方程的初值问题微分方程的初值问题. . 4、了解计算机数据仿真、了解计算机数据仿真、 数据模拟的基本方法。数据

2、模拟的基本方法。 实验十二实验十二 缉私艇追赶走私船模型实验缉私艇追赶走私船模型实验 w微分方程数值解微分方程数值解( (理论理论). ). w用用MATLABMATLAB软件求微分方程解析解、数值解软件求微分方程解析解、数值解. . w微分方程模型实验微分方程模型实验: :缉私艇追赶走私船缉私艇追赶走私船. . w问题的进一步扩展问题的进一步扩展. . 实验内容实验内容 1. 1. 微分方程数值解微分方程数值解 一阶常微分方程的初值问题一阶常微分方程的初值问题: : 数值解方法数值解方法: : 欧拉方法、欧拉两步法、改进的欧拉方法欧拉方法、欧拉两步法、改进的欧拉方法 w欧拉方法：欧拉方法：是

3、一种求解给定初值的常微分方是一种求解给定初值的常微分方 程的一阶数值解方法。程的一阶数值解方法。 对于方程对于方程 设设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)，即，即 )(,()( 0100011 xxyxfyyxy 然后再过点然后再过

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

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

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

7、用欧拉公式提供一个通常先用欧拉公式提供一个ynyn附近的初始迭代值，然后再用梯形附近的初始迭代值，然后再用梯形 公式做迭代。公式做迭代。 2 2）欧拉两步法欧拉两步法 用函数用函数y(x)y(x)在在xnxn处的中心差商代替微分方程中的处的中心差商代替微分方程中的 导数，导数， )(,( 2 )()( 11 nn nn xyxf h xyxy 这就是这就是欧拉两步法公式。欧拉两步法公式。 )(,(2)()( 11nnnn xyxhfxyxy 它比欧拉公式、后退的欧拉公式具有更高的精度。它比欧拉公式、后退的欧拉公式具有更高的精度。 3 3）改进的欧拉公式改进的欧拉公式 欧拉公式虽然简单，但是误差

8、较大，梯形公式虽欧拉公式虽然简单，但是误差较大，梯形公式虽 然然 提高了精度，但因其是隐式的，所以计算麻烦，提高了精度，但因其是隐式的，所以计算麻烦， 工工 作量较大。通常将这两种方法结合，得到作量较大。通常将这两种方法结合，得到改进的改进的 欧欧 拉公式，拉公式，具体如下：具体如下： ),( 1nnnn yxhfyy 步骤步骤1 1 先由欧拉公式得到先由欧拉公式得到y(xn+1)y(xn+1)的一个初始近似的一个初始近似 值，称为值，称为预测值预测值 步骤步骤2 2 对预测值再用梯形公式校正一次得对预测值再用梯形公式校正一次得yn+1yn+1，称，称 为为校正值校正值 ),(),( 2 11

9、1 nnnnnn yxfyxf h yy 整个过程（称整个过程（称预测预测- -校正系统校正系统）可以表示成）可以表示成 ),(,(),( 2 1nnnnnnnn yxhfyhxfyxf h yy 这就是这就是改进的欧拉公式。改进的欧拉公式。 当然也可以将欧拉两步法与梯形方法结合。当然也可以将欧拉两步法与梯形方法结合。 再解释：再解释： 则用函数则用函数y(x)y(x)在区间在区间xn,xn+1xn,xn+1上的平均斜率代替上的平均斜率代替 了曲线在点了曲线在点xnxn的斜率。的斜率。 上上的的平平均均斜斜率率。在在区区间间为为函函数数称称,)( )()( 1 1 nn nn xxxy h x

10、yxy 得到启发：得到启发：只要对曲线在区间只要对曲线在区间xn,xn+1xn,xn+1上的平均上的平均 斜率提供一种算法，就可以得到一种计算斜率提供一种算法，就可以得到一种计算yn+1yn+1的公的公 式。如果设法在区间内多预测几个点的斜率值，然式。如果设法在区间内多预测几个点的斜率值，然 后将这些点处后将这些点处斜率值的加权平均值斜率值的加权平均值作为曲线在区间作为曲线在区间 上的平均值，由此构造出由上的平均值，由此构造出由ynyn计算计算yn+1yn+1的精度更高的精度更高 的计算公式，这就是的计算公式，这就是龙格龙格- -库塔方法库塔方法的基本思想。的基本思想。 2. 2. 用用MAT

11、LABMATLAB软件求微分方程的解析解软件求微分方程的解析解 例例1 dsolve(D2y=exp(x),x) ans=exp(x)+C1*x+C2 dsolve(D2y=exp(x) ans =1/2*exp(x)*t2+C1*t+C2 dsolve(D2y=exp(x),y(0)=1,Dy(0)=2,x) ans =exp(x)+x 3. 3. 编程计算微分方程数值解编程计算微分方程数值解 例例2 2 （1 1）直接在命令窗执行命令）直接在命令窗执行命令 dsolve(Dy=3/x*y+x3*(exp(x)+cos(x)-2*x, y(pi)=(exp(pi)+2/pi)*pi3,x)

12、ans= x3*exp(x)+x3*sin(x)+2*x2 szy_eu=; y=(exp(pi)+2/pi)*pi3; for x=pi:0.1:2*pi y=y+0.1*(3/x*y+x3*(exp(x)+cos(x)-2*x); szy_eu=szy_eu,y； end szy_eu t=pi:0.1:2*pi; f=t.3.*exp(t)+t.3.*sin(t)+2*t.2; plot(t,szy_eu,r-,linewidth,2) hold on plot(t,f,b*-) (2) (2) 取步长取步长0.10.1，用欧拉公式求数值解，用欧拉公式求数值解 -zhao1201.m-z

13、hao1201.m 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); 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,

14、r-,linewidth,2) f=x.3.*exp(x)+x.3.*sin(x)+2*x.2; hold on plot(t,f,b*-) (3) (3) 取步长取步长0.10.1，用改进的欧拉公式求数值解，用改进的欧拉公式求数值解 -zhao1202zhao1202 4. 4. 用用MATLABMATLAB软件求微分方程数值解软件求微分方程数值解 MATLABMATLAB软件求常微分方程数值解的指令有软件求常微分方程数值解的指令有 ode23ode23，de23tde23t和和ode45ode45，分别表示二三阶龙格，分别表示二三阶龙格 库塔法，二三阶龙格库塔法，二三阶龙格库塔梯形法和四五

15、阶龙库塔梯形法和四五阶龙 格格库塔法。最简单的使用格式为库塔法。最简单的使用格式为 ode23(fname,xs,xe,sv) 其中用单引号引起的是存储微分方程的函数文件名，其中用单引号引起的是存储微分方程的函数文件名， xsxs表示自变量的初始值，表示自变量的初始值，xexe表示自变量的终止值，表示自变量的终止值，svsv 表示迭代初始向量值。该命令执行后自动画出微分方表示迭代初始向量值。该命令执行后自动画出微分方 程数值解曲线。程数值解曲线。 例例3 3 使用格式使用格式: fc.m: fc.m function f=fc(x,y)function f=fc(x,y) f=2f=2* *y

16、+x+2y+x+2 在在MATLABMATLAB命令窗中执行命令命令窗中执行命令 ode23(fc,0,1,1)ode23(fc,0,1,1) 5. 5. 模型实验模型实验: : 缉私艇追赶走私船缉私艇追赶走私船 海上边防缉私艇发现距海上边防缉私艇发现距c公里处有一走私船正以匀公里处有一走私船正以匀 速速a沿直线行驶沿直线行驶, ,缉私艇立即以最大速度缉私艇立即以最大速度b追赶追赶, ,在在 雷雷 达的引导下，缉私艇的方向始终指向走私船。问达的引导下，缉私艇的方向始终指向走私船。问 缉缉 私艇何时追赶上走私船私艇何时追赶上走私船? ?并求出缉私艇追赶的路线并求出缉私艇追赶的路线 方程。方程。

17、x y co (1) (1) 建立模型建立模型 走私船初始位在点走私船初始位在点(0,0)，方向为，方向为y轴正方向，轴正方向， 缉私艇的初始位在点缉私艇的初始位在点(c,0)，在时刻，在时刻t，走私船的，走私船的 位置到达点位置到达点R(0,at)，缉私艇的位置到达，缉私艇的位置到达D(x,y) ), 0(atR ),(yxD 可得缉私艇追击走私船过程的数学模型：可得缉私艇追击走私船过程的数学模型： (2) (2) 模型求解模型求解 求解析解求解析解 令：令： 1)1 a b r -zhao1203.m 表示追赶时间。表示追赶时间。 表示走私船跑过的距离。表示走私船跑过的距离。 取取c=3c

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

19、.5,用二三阶龙格用二三阶龙格- -库塔算法计算数值解库塔算法计算数值解: : 6 6 动态仿真动态仿真 当建立动态系统的微分方程模型很困难时，当建立动态系统的微分方程模型很困难时， 可以用计算机仿真法对系统进行分析研究。所可以用计算机仿真法对系统进行分析研究。所 谓计算机仿真就是利用计算机对实际动态系谓计算机仿真就是利用计算机对实际动态系 统的结构和行为进行编程、模拟和计算，以此统的结构和行为进行编程、模拟和计算，以此 来预测系统的行为效果。来预测系统的行为效果。 走私船初始位在点走私船初始位在点(0,0)(0,0)，方向为，方向为y y轴正方向，轴正方向， 缉私艇的初始位在点缉私艇的初始位

20、在点( (c c,0),0)， 走私船的位置到达点走私船的位置到达点), 0( k atR ，在在时时刻刻 k t 缉私艇的位置到达点缉私艇的位置到达点 ),( kk yxD 追赶方向可用方向余弦表示为追赶方向可用方向余弦表示为： ), 0( k atR k ),( kk yxD 22 )()0( 0 cos kkk k k yatx x 22 )()0( sin kkk kk k yatx yat , 1 时时tttt kk 时间步长时间步长 缉私艇的位置：缉私艇的位置： ).,( 11 kk yx则则 ,cos 1kkkk tbxxx kkkk tbyyy sin 1 ), 0( 1 k

21、at ), 0( k atR k ),( kk yxD ).,( 11 kk yx 仿真算法：仿真算法： 第二步：第二步： 计算动点缉私艇计算动点缉私艇D D在时刻在时刻 ttt kk 1 时的坐标时的坐标 22 1 )( kkk k kk yatx x tbxx 22 1 )( kkk k kk yatx yat tbyy 计算走私船计算走私船R在时刻在时刻 ttt kk 1 时的坐标时的坐标 ) , ( 11 kk yx 0 1 k x )( 1 ttay kk 第一步：设置时间步长第一步：设置时间步长 , 速度速度a, b及初始位置及初始位置 t 0, 0, 0 kycx 仿真算法步骤：

22、仿真算法步骤： 第三步第三步：计算缉私艇与走私船这两个动点之间的距离：计算缉私艇与走私船这两个动点之间的距离： 2 11 2 11 ) () ( kkkkk yyxxd 根据事先给定的距离，判断缉私艇是否已经追根据事先给定的距离，判断缉私艇是否已经追 上了走私船，从而判断退出循环还是让时间产生一上了走私船，从而判断退出循环还是让时间产生一 个步长，返回到第二步继续进入下一次循环；个步长，返回到第二步继续进入下一次循环； 第四步第四步：当从上述循环退出后，由点列当从上述循环退出后，由点列 ),( 11 kk yx 和和 ) , ( 11 kk yx 可分别绘制成两条曲线可分别绘制成两条曲线 即为缉私艇和走私船走过的轨迹曲线。即为缉私艇和走私船走过的轨迹曲线。 取取c=3c=3千米，千米，a a=0.4=0.4千米千米/ /分钟，分钟，b b=0.8=0.8千米千米/ /分钟，分钟，r=a/b=0.5r=a/b=0.5 显示船与艇行进路线程序显示船与艇行进路线程序 clc;clear; 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,j