数学实验：实验十二 缉私艇追赶走私船模型实验

MathematicsLaboratory阮小娥博士ExperimentsinMathematics数学实验2、了解微分方程数值解思想，掌握两种简单的微分方程数值解方法。3、学会根据实际问题建立简单微分方程数学模型。实验目的1、学会用MATLAB软件求解微分方程的初值问题.4、了解计算机数据仿真、数据模拟的基本方法。实验十二缉私艇追赶走私船模型实验用MATLAB软件求微分方程解析解、数值解理论.编程求解微分方程数值解.微分方程模型实验:缉私艇追赶走私船.实验内容1.微分方程的解一阶常微分方程的初值问题:

用matlab软件求微分方程的解析解dsolve(‘equation’,‘condition’)

%求方程equation在初始条件condition下的解，自变量默认为tdsolve('equation')%求方程equation的通解%一阶导数用Dy表示，二阶导数用D2y表示，

A=dsolve('Dy=5')B=dsolve(‘Dy=x’,‘x’)%求方程的通解，指定自变量为xC=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')

微分方程数值解的计算

欧拉方法、欧拉两步法、改进的欧拉方法欧拉方法：是一种求解给定初值的常微分方程的一阶数值解方法。对于方程设y=y(x)是它的解，则在解y=y(x)的曲线上每一点处的切线的斜率等于f(x,y)在该点处的值。所以，过初始点A0(x0,y0)，做切线与直线x=x1交点的纵坐标y1近似代替y=y(x)在x=x1处的函数值y(x1)，即然后再过点A1(x1,y1)，以f(x1,y1)为斜率做切线，与直线x=x2交点的纵坐标y2近似代替y=y(x)在x=x2处的函数值y(x2)，即A0x0x1x2A1A2这样依次类推得到上述微分方程的解在x=xn+1处函数值的近似计算公式，待求的曲线为蓝色，它的折线近似为红色A0x0x1x2x3A1A2A3A4x4特别地，如果节点之间是等分的，则有这就是著名的欧拉(Euler)公式。事实上，上述公式是用函数y(x)在xn处的向前差商代替函数在这点处的导数，从而将微分方程离散化。同样，如果用函数y(x)在xn+1处的向后差商代替导数值，就得后退的欧拉公式注意：这两个公式虽然具有相同的精度，但也是有区别的，前者是计算yn+1的显示公式，后者是隐式公式。将上两式算术平均即得梯形公式：梯形公式消去了欧拉公式和后退的欧拉公式误差的主要部分，因此比它们具有更高的精度。这是一个隐式公式，可采用迭代法求解。通常先用欧拉公式提供一个yn附近的初始迭代值，然后再用梯形公式做迭代。2）欧拉两步法用函数y(x)在xn处的中心差商代替微分方程中的导数，这就是欧拉两步法公式。它比欧拉公式、后退的欧拉公式具有更高的精度。3）改进的欧拉公式欧拉公式虽然简单，但是误差较大，并且随着迭代次数的增加，误差越来越大。梯形公式虽然提高了精度，但因其是隐式的，所以计算麻烦，工作量较大。通常将这两种方法结合，得到改进的欧拉公式，具体如下：步骤1先由欧拉公式得到y(xn+1)的一个初始近似值，称为预测值步骤2对预测值再用梯形公式校正一次得yn+1，称为校正值整个过程（称预测-校正系统）可以表示成这就是改进的欧拉公式。当然也可以将欧拉两步法与梯形方法结合。再解释：则用函数y(x)在区间[xn,xn+1]上的平均斜率代替了曲线在点xn的斜率。得到启发：只要对曲线在区间[xn,xn+1]上的平均斜率提供一种算法，就可以得到一种计算yn+1的公式。如果设法在区间内多预测几个点的斜率值，然后将这些点处斜率值的加权平均值作为曲线在区间上的平均值，由此构造出由yn计算yn+1的精度更高的计算公式，这就是龙格-库塔方法的基本思想。编程计算微分方程数值解例2

已知初值问题（1）用MATLAB软件求解析解，算出解析解在结点x=pi:0.1:2*pi处的纵坐标；（2）取步长0.1，用欧拉公式求数值解；（3）取步长0.1，用改进的欧拉公式求数值解；（4）比较两种数值解和解析解，并画出数值解和解析解曲线。解：（1）直接在命令窗执行命令dsolve('Dy=3/x*y+x^3*(exp(x)+cos(x))-2*x','y(pi)=(exp(pi)+2/pi)*pi^3','x')ans=x^3*exp(x)+x^3*sin(x)+2*x^2clc;clf;szy_eu=[];y=(exp(pi)+2/pi)*pi^3;forx=pi:0.1:2*piy=y+0.1*(3/x*y+x^3*(exp(x)+cos(x))-2*x);szy_eu=[szy_eu,y];endszy_eut=pi:0.1:2*pi;f=t.^3.*exp(t)+t.^3.*sin(t)+2*t.^2;plot(t,szy_eu,'r-','linewidth',2)holdonplot(t,f,'b*--')(2)取步长0.1，用欧拉公式求数值解zhao1201.mx=pi:0.1:2*pi;h=0.1;y=(exp(pi)+2/pi)*pi^3;szy_imeu=[y];fori=1:length(x)-1jzj=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;endszy_imeuplot(x,szy_imeu,'r-','linewidth',2)f=x.^3.*exp(x)+x.^3.*sin(x)+2*x.^2;holdonplot(t,f,'b*--')(3)取步长0.1，用改进的欧拉公式求数值解

zhao12022.用MATLAB软件求微分方程数值解MATLAB软件求常微分方程数值解的指令有ode23，ode23t和ode45，分别表示二三阶龙格—库塔法，二三阶龙格—库塔梯形法和四五阶龙格—库塔法。还有其他命令ode113，ode15i，ode15s，ode23s等，读者可通过help命令查询。一般情况下，ode45,ode23,ode113是解决非刚性（non-stiff）问题的命令，ode15s和ode23s是求解刚性（stiff）问题的命令。[xb,yb]=ode23('fname',[xs,xe],sv)返回结点处的横坐标xb和纵坐标yb。ode23('fname',[xs,xe],sv)其中用单引号引起的是存储微分方程的函数文件名，xs表示自变量的初始值，xe表示自变量的终止值，sv表示迭代初始向量值。该命令执行后自动画出微分方程数值解曲线。最简单的使用格式为注意：（1）命令ode求解的是形如y’=f(t,y)的微分方程，我们称它为一阶导数可解出的微分方程。而对于一阶导数解不出的形如

f(t,y,y‘)=0

的微分方程，可以用命令ode15i求解，有兴趣的读者可以用help命令获得。（2）求解微分问题时必须为其指定初值，即上面的sv。

（3）ode只能直接求解一阶微分方程，高阶微分方程必须等价地转化成一阶微分方程组，才能用ode命令求解。比如y‘’‘=f(t,y,y’,y‘’)，

将y,y‘,y’‘设为y1,y2,y3，那么如上的三阶微分方程就可以表示为如下的三个一阶方程组了。

y1'=y2

y2'=y3

y3'=f(t,y1,y2,y3)例3使用格式:fc.mfunctionf=fc(x,y)f=2*y+x+2;在MATLAB命令窗中执行命令ode23('fc',[0,1],1)例4求解微分方程建立m文件fff.mfunctiondy=fff(t,y)dy=zeros(2,1);dy(1)=y(2);dy(2)=1000*(1-y(1)^2)*y(2)-y(1);在MATLAB命令窗中执行命令[t,y]=ode15s('fff',[0,3000],[2,0])plot(t,y(:,1),'-')4.模型实验:缉私艇追赶走私船海上边防缉私艇发现距c公里处有一走私船正以匀速a沿直线行驶,缉私艇立即以最大速度b追赶,在雷达的引导下，缉私艇的方向始终指向走私船。问缉私艇何时追赶上走私船?并求出缉私艇追赶的路线方程。xyco解

(1)建立模型走私船初始位在点(0,0)，方向为y轴正方向，缉私艇的初始位在点(c,0)，在时刻t，走私船的位置到达点R(0,at)，缉私艇的位置到达D(x,y)可得缉私艇追击走私船过程的数学模型：(2)模型求解求解析解令：--zhao1203.m表示追赶时间。表示走私船跑过的距离。取c=3千米，a=0.4千米/分钟，分别取b=0.6,0.9,1.2千米/分钟时，缉私艇追赶路线的图形为定义m文件函数：functiony=zx(t,y)y=0.5*((t/3)^0.5-(3/t)^0.5);ode23('zx',3,0.0005,0)

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

追赶方向可用方向余弦表示为：时间步长缉私艇的位置：则仿真算法：

第二步：计算动点缉私艇D在时刻

时的坐标

计算走私船R在时刻

时的坐标第一步：设置时间步长,速度a,b及初始位置仿真算法步骤：第三步：计算缉私艇与走私船这两个动点之间的距离：

根据事先给定的距离，判断缉私艇是否已经追上了走私船，从而判断退出循环还是让时间产生一个步长，返回到第二步继续进入下一次循环；

第四步：当从上述循环退出后，由点列

和可分别绘制成两条曲线即为缉私艇和走私船走过的轨迹曲线。

取c=3千米，a=0.4千米/分钟，b=0.8千米/分钟，r=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);holdon

plot(zscx,zscy,'b*','markersize',15);pause(0.01)t=t+dt;

jstx=jstx-b*dt*jstx/sqrt(jstx^2+(a*t-jsty)^2);

jsty=jsty+b*dt*(a*t-jsty)/sqrt(jstx^2+(a*t-jsty)^2);

zscy=a*t;endfprintf('追上所需时间t=%.0f\n',t);fprintf('缉私艇的位置=(%.5f,%.5f)\n',jstx,jsty);fprintf('走私