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

最简单的使用格式为
ode23('fname',[xs,xe],sv)
其中用单引号引起的是存储微分方程的函数文件名， xs表示自变量的初始值，xe表示自变量的终止值，sv 表示迭代初始向量值。该命令执行后自动画出微分方 程数值解曲线。
[xb,yb]=ode23('fname',[xs,xe],sv)
解：（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^2
(2) 取步长0.1，用欧拉公式求数值解
clc;clf; szy_eu=[]; y=(exp(pi)+2/pi)*pi^3; for x=pi:0.1:2*pi y=y+0.1*(3/x*y+x^3*(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,f,'b*--') hold on plot(t,szy_eu,'r-','linewidth',2)
海上边防缉私艇发现距c公里处有一走私船正以匀 速a沿直线行驶,缉私艇立即以最大速度b追赶,在雷 达的引导下，缉私艇的方向始终指向走私船。问缉 私艇何时追赶上走私船?并求出缉私艇追赶的路线 方程。 y
o
c
x
解
(1) 建立模型
走私船初始位在点(0,0)，方向为y轴正方向， 缉私艇的初始位在点(c,0)，在时刻t，走私船的 位置到达点R(0,at)，缉私艇的位置到达D(x,y)
y n1 yn hf ( xn , yn )
步骤2 对预测值再用梯形公式校正一次得yn+1，称 为校正值
yn1
h yn [ f ( x n , yn ) f ( x n1 , y n1 )] 2
整个过程（称预测-校正系统）可以表示成
yn1
h yn [ f ( x n , yn ) f ( x n h, yn hf ( x n , yn ))] 2
得到启发：只要对曲线在区间[xn,xn+1]上的平均 斜率提供一种算法，就可以得到一种计算yn+1的公 式。如果设法在区间内多预测几个点的斜率值，然 后将这些点处斜率值的加权平均值作为曲线在区间 上的平均值，由此构造出由yn计算yn+1的精度更高 的计算公式，这就是龙格-库塔方法的基本思想。
编程计算微分方程数值解
A2
A1 A0
x0 x1 x2
A3
A4
y( x n 1 ) yn 1 yn hf ( x n , yn )
这就是著名的欧拉(Euler)公式。
x3
x4
事实上，上述公式是用函数y(x)在xn处的向前差
商代替函数在这点处的导数，从而将微分方程离 散化。同样，如果用函数y(x)在xn+1处的向后差 商代替导数值，就得后退的欧拉公式
2. 用MATLAB软件求微分方程数值解
MATLAB软件求常微分方程数值解的指令有 ode23，ode23t和ode45，分别表示二三阶龙 格—库塔法，二三阶龙格—库塔梯形法和四五阶 龙格—库塔法。还有其他命令ode113，ode15i， ode15s，ode23s等，读者可通过help命令查询。
例3
使用格式: fc.m function f=fc(x,y) f=2*y+x+2; 在MATLAB命令窗中执行命令 ode23('fc',[0,1],1)
例4
求解微分方程
d 2x dx 2 1000 (1 x ) x 0 2 dt dt x(0) 2, x (0) 0
设y=y(x)是它的解，则在解y=y(x)的曲线上每一点
处的切线的斜率等于f(x,y)在该点处的值。
A2
所以，过初始点A0(x0,y0)，做切线
A1 A0
x0 x1 x2
y y0 f ( x0 , y0 )( x x0 )
x=x1处的函数值y(x1)，即
与直线x=x1交点的纵坐标y1近似代替y=y(x)在
实验内容
1.用MATLAB软件求微分方程解析解、数值解
理论.
2.编程求解微分方程数值解.
3.微分方程模型实验:缉私艇追赶走私船.
1. 微分方程的解
一阶常微分方程的初值问题:
用matlab软件求微分方程的解析解


dsolve(‘equation’,‘condition’) %求方程equation在 初始条件condition下的解，自变量默认为t dsolve('equation') %求方程equation的通解 %一阶导数用Dy表示，二阶导数用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') %解微分方程组
返回结点处的横坐标xb和纵坐标yb。
注意：
（1）命令ode求解的是形如y’=f(t,y)的微分方 程，我们称它为一阶导数可解出的微分方程。 而对于一阶导数解不出的形如 f(t,y,y’)=0 的 微分方程，可以用命令ode15i求解，有兴趣的 读者可以用help命令获得。
（2）求解微分问题时必须为其指定初值，即上 面的sv。
数学实验
Experiments in Mathematics Laboratory Mathematics
阮小娥博士 赵小艳
办公地址：理科楼214
实验十二 缉私艇追赶走私船模型实验
实验目的
1、学会用MATLAB软件求解 微分方程的初值问题.
2、了解微分方程数值解思 想，掌握两种简单的微分 方程数值解方法。 3、学会根据实际问题建立 简单微分方程数学模型。 4、了解计算机数据仿真、 数据模拟的基本方法。
y( x1 ) y1 y0 f ( x0 , y0 )( x1 x0 )
然后再过点A1(x1,y1)，以f(x1,y1)为斜率做切 线，与直线x=x2交点的纵坐标y2近似代替y=y(x)在 x=x2处的函数值y(x2)，即
y( x2 ) y2 y1 f ( x1 , y1 )( x2 x1 )
（3）ode只能直接求解一阶微分方程，高阶微分方 程必须等价地转化成一阶微分方程组，才能用ode 命令求解。比如
y’’’=f(t,y,y’,y’’)， 设y1=y, y2= y’, y3= y’’ ，
那么如上的三阶微分方程就可以表示为如下的三 个一阶方程组了。 y1'=y2 y2'=y3 y3'=f(t,y1,y2,y3)
建立m文件 erjie.m function f=erjie(t,y) f=[y(2); 1000*(1-y(1)^2)*y(2)-y(1)]; 在MATLAB命令窗中执行命令
[t,y]=ode15s('erjie',[0,3000],[2,0]) plot(t,y(:,1),'-')
4. 模型实验: 缉私艇追赶走私船
y( x n 1 ) yn 1 yn hf ( x n 1 , yn1 )
注意：这两个公式虽然具有相同的精度，但也是有区别的， 前者是计算yn+1的显示公式，后者是隐式公式。
将上两式算术平均即得梯形公式： h y( x n 1 ) yn1 yn [ f ( x n , yn ) f ( x n1 , yn 1 )] 2 梯形公式消去了欧拉公式和后退的欧拉公式误差 的主要部分，因此比它们具有更高的精度。 这是一个隐式公式，可采用迭代法求解。
5
例1 (1)
dsolve('D2y=exp(x)','x') dsolve('D2y=exp(x)') dsolve('D2y=exp(x)','y(0)=1,Dy(0)=2','x')
微分方程数值解的计算
欧拉方法、欧拉两步法、改进的欧拉方法
1) 欧拉方法：是一种求解给定初值的常微分方
程的一阶数值解方法。 对于方程
通常先用欧拉公式提供一个yn附近的初始迭代值， 然后再用梯形公式做迭代。
y
y
(0) n 1
yn hf ( xn , yn )
h (k ) yn [ f ( x n , yn ) f ( x n 1 , yn 1 )] 2
( k 1 ) n 1
2）欧拉两步法 用函数y(x)在xn处的中心差商代替微分方程中的 导数，
这就是改进的欧拉公式。
当然也可以将欧拉两步法与梯形方法结合。
再解释：
y( xn1 ) y( xn ) 称 为函数y( x )在区间 [ xn , xn1 ]上的平均斜率。 h
欧拉公式用函数y(x)在区间[xn,xn+1]上的平均斜 率代替了曲线在点xn的斜率。 欧拉两步法用区间[xn,xn+1]和区间[xn-1,xn]上的 平均斜率的算术平均值代替了曲线在xn的斜率。