第四章 微分方程与差分方程方法

第四章 微分方程与差分方程方法
第一节 微分方程模型
我们在数学分析中所研究的函数,是反映客观现实世界运动过程中量与量之间的一种关系,但我们在构造数学模型时,遇到的大量实际问题往往不能直接写出量与量之间的关系,却能比较容易地建立这些变量和它们的导数(或微分)间的关系式,这种联系着自变量、未知函数及其导数(或微分)的关系式称为微分方程.
§4.1.1 微分方程简介
这一节,我们将介绍关于微分方程的一些基本概念. 一、微分方程的阶数
首先我们具体的来看一个微分方程的例子.
例4-1 物体冷却过程的数学模型
将某物体放置于空气中,在时刻0=t ,测量得它的温度为C u 00150=,10分钟后测量得温度为C u 01100=.我们要求决定此物体的温度u 和时间t 的关系,并计算20分钟后物体的温度.这里我们假定空气的温度保持为C u 024=α. 解:根据物理学中的牛顿冷却定律可知,热量总是从温度高的物体向温度低的物体传导;一个物体的温度变化速度与这一物体的温度与其所在介质温度的差值成
正比.设物体在时刻t 的温度为)(t u u =,则温度的变化速度可以用dt
du
来表示.我
们得到描述物体温度变化的微分方程
)(αu u k dt
du
--= (4.1.1) 其中0>k 是比例常数.
方程(4.1.1)中含有未知函数u 及它的一阶导数dt
du
,这样的方程,我们称为一
阶微分方程.微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数.
方程
)(33t f cy dt dy
b dt
y d =++ (4.1.2)
中未知函数最高阶导数的阶数是三阶,则方程(4.1.2)称为三阶微分方程. 二、常微分方程与偏微分方程
如果在微分方程中,自变量的个数只有一个,我们称这种微分方程为常微分
方程;自变量的个数为两个或两个以上的微分方程称为偏微分方程.
方程
0222222=∂∂+∂∂+∂∂z T
y T x T (4.1.3)
就是偏微分方程的例子,其中T 是未知函数,x 、y 、z 都是自变量.而方程(4.1.1)(4.1.2)都是常微分方程的例子. 三、线性与非线性微分方程
如果n 阶常微分方程
0),,,,(=n n dx
y
d dx dy y x F (4.1.4)
的左端为关于未知函数y 及其各阶导数的线性组合,则称该方程为线性微分方程,否则称为非线性方程.一般的n 阶线性微分方程具有形式
)()()()(1111x f y x a dx dy
x a dx
y d x a dx y d n n n n n n =++++--- (4.1.5)
其中)1( )(),(n i x f x a i =是关于x 的已知函数.
当()0f x =时,称(4.1.5)为n 阶齐次线性微分方程;当()0f x ≠时,称(4.1.5)为
n 阶非齐次线性微分方程.
例如,方程(4.1.2)是三阶线性微分方程.而方程
2
0dy dy t y dx dx ⎛⎫
++= ⎪⎝⎭ (4.1.6)
是一阶齐次非线性方程. 四、微分方程初边值问题
我们把含有n 个独立的任意常数12,,n c c c 的解
12(,,,)n y x c c c ϕ=
称为n 阶微分方程(4.1.4)的通解.为了确定微分方程一个特定的解,我们通常给出
这个解所必须满足的条件,这就是定解条件.常见的定解条件是初始条件.n 阶微分方程(4.1.4)的初始条件是指如下的n 个条件:
1(1)(1)
0000001
()()(),,n n n dy x d y x y x y y y dx dx
---=== 求解微分方程满足初始条件的解,称为初值问题.求解初值问题的过程,就是
通过初始条件确定通解中的常数,从而求得满足初始条件的特解的过程.
若方程所给出的定解条件,既有自变量初始时刻的值,也有自变量取终值时的值,则称该问题为边值问题.
§4.1.2微分方程的求解及Matlab 实现
线性微分方程和低阶特殊微分方程往往可以通过解析解的方法求解,但一般的非线性微分方程是没有解析解的,即使可以求得解析解,参数也很难确定,需要用数值解的方式求解.具体的求解方法很多,本节我们主要介绍如何利用matlab 来求解微分方程的解析解和数值解.
§4.1.2.1微分方程的解析解
一、基本理论
有些简单的常微分方程可用一些技巧,如分离变量法、积分因子法、常数变异法、降阶法等,化为可直接积分的方程从而求得显示解.例如,一阶常系数线性常微分方程
(0)dy
ay b a dt
=+≠ 可化为
dy
dt ay b
=+ 两边积分可得()y t 的通解为 at b
y Ce a
=-.
一般的常系数线性微分方程
1111()n n n n n n d y d y dy
a a a y f x dx dx dx
---++++= (4.1.7)
的解满足叠加原理,即方程(4.1.7)的通解是对应齐次方程
11110n n n n n n d y d y dy
a a a y dx dx dx
---++++= (4.1.8)
的通解与非齐次方程(4.1.7)的一个特解的和.
一阶常系数线性常微分方程总可用这一思路求得显示解.
高阶线性常系数微分方程可用特征根法求对应齐次微分方程(4.1.8)的通解.(4.1.8)对应的代数特征方程
121210n n n n n s a s a s a s a ---+++++=
若方程的特征根i s 均可以求出,且两两相异,则方程(4.1.7)的解可以表示为
1212()()n s x s x s x n y x C e C e C e x γ=++++
其中,i C 为待定系数,()x γ是满足方程(4.1.8)的一个特解,这个特解可以通过常数变异法,比较系数法得到.i s 有重根的情况也可以写出相应的解析解形式,这里我
们不做过多介绍,因为我们下面要介绍一种更简单的计算机求解法. 二、用matlab 求解线性常系数微分方程解析解 dsolve 其中,i f 既可以描述微分方程,也可以描述初边值条件.在描述微分方程时,用字母D 代表微分,后面加数字表示微分的阶数,例如D4y 表示(4)()y t ,也可以用D2y(0)=3 这类记号来表示初值条件 ''(0)3y =.一般而言任何D 后面所跟的字母为因变量,自变量默认为t ,也可以在函数调用时指明.
例4-2 求方程2223t d x dx
x e dt dt
---=的通解.
解:输入命令
>>syms t,x; %定义符号变量
x=dsolve(‘D2x -2*Dx-3*x=exp(-t)’) %求解方程
latex(x) %用 LATEX 语句显示结果 运行结果
x
=exp(-t)*C2+exp(3*t)*C1-1/4*t*exp(-t) ans
={e^{3t}}{\it C2}+{e^{-t}}{\it C1}-1/4\,t{e^{-t}} 原方程的通解为
3121
4
t t t x C e C e te --=+-
其中,i C 为任意常数.若给出初始或边界条件,则可以通过这些条件建立方程,求出i C 的值.
仍考虑上面的微分方程,假设已知 (0)3,(0)2dx
x dt
==,则可以通过下面的命令求满足该方程的特解.
>> x=dsolve(‘D2x -2*Dx-3*x=exp(-t)’,’x(0)=3’,’Dx(0)=2’) 运行结果
x
=27/16*exp(-t)+21/16*exp(3*t)-1/4*t*exp(-t)
所以，满足这一初值问题的特解是
321271
16164
t t t x e e te --=+- 对于线性方程组我们可以通过定义多维的因变量同样求解. 例4-3 试求解下面的线性微分方程
233453442dx
x y z dt dy
x y z dt dz
x y z dt ⎧=-+⎪⎪⎪=-+⎨⎪⎪=-+⎪⎩
解:输入命令
>>[x,y,z]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z')； x=collect(simple(x)) %化简x 的表达式 y=collect(simple(y)) z=collect(simple(z)) 结果为
部分特殊的非线性微分方程也可以用dsolve()函数来求解析解.这样的方程描述方式和前面介绍的线性微分方程是一致的.下面我们通过例子来演示非线性方程的解析解求解问题,同时还将演示不能求解的例子. 例4-4 试求出一阶非线性微分方程2()(1())dx
x t x t dt
=-的解析解. 解: 输入命令
>>x=dsolve('Dx=x*(1-x^2)') 运行结果
x =
1/(1+exp(-2*t)*C1)^(1/2) -1/(1+exp(-2*t)*C1)^(1/2)
即该方程的解析解为()x t =但是如果稍微改变原方程,例如将等号右侧加上1,则可以用下面语句试解方程.我们会发现该方程是没有解析解的.
223222312231t t
t t t t t
x C e C e y C e C e C e z C e C e ----=+=++=+
>>x=dsolve('Dx=x*(1-x^2)+1') 运行结果
Warning: Explicit solution could not be found; implicit solution returned. > In dsolve at 312 x =
t+Int(-1/(_a-_a^3+1),_a = .. x)+C1 = 0
从上例我们可以看出,可以求得解析解的方程是很有限的,对于一般的非线性方程只能用数值解法去求解.下面我们介绍微分方程数值解的求解方法.
§4.1.2.2微分方程的数值解
一、基本理论 (一) 数值解的定义
考虑一阶常微分方程组初值问题
00
0'(,)
()f y f t y t t t y t y =⎧<<⎨=⎩ (4.1.9)
其中 12(,,,)m T y y y y = ,12(,,,)m T f f f f = ,120000(,,,)m T
y y y y = ,这里T 表示
转置.微分方程的数值解,是指不求出解()y t 的解析表达式,而是寻求()y t 在一系列离散结点01n f t t t t <<<≤ 上的近似值(0,1,,)k y k n = .称1k k k h t t +=-为步长.通常,取为常数0()/n h t t n =-,此时结点变为等距结点,有
10k k t t h t kh -=+=+.
(二) 求解数值解的方法 1.欧拉法
基本思想:在结点处用差商近似替代导数
()()
'()(,())(,)k k k k k k k y t h y t y t f t y t f t y h
+-≈=≈
这样导出欧拉格式的计算公式
100(,)
0,1,1() k k k k y y hf t y k n y y t +=+⎧=-⎨=⎩ (4.1.10)
欧拉法可以求解各种形式的微分方程,但是它只有一阶精度,所以实际应用
效率比较差.
2.改进的欧拉法
基本思想: 使用数值积分离散化.对方程'(,)y f t y =,两边由k t 到1k t +积分,得到
1
1()()(,())k k
t k k t y t y t f t y t dt ++-=⎰
对右端的定积分用数值积分方法做离散化,可得计算公式,如用矩形公式可得欧拉公式,若用梯形公式可得改进的欧拉公式,下面我们用梯度公式来做数值积分
11111()()(,())[(,())(,())]2
k k
t k k k k k k k k t t t
y t y t f t y t dt f t y t f t y t +++++--=≈+⎰
故有如下迭代计算公式:
11100[(,)(,)]2
()k k
k k k k h y y f t y f t y y y t +++⎧
⎪⎨⎪⎩
=++= (4.1.11) 实际应用时,与欧拉公式结合使用:
(0)1
(1)()
111(,)[(,)(,)] 0,1,2,2
k k k k l l k k k k k k y y hf t y h y y f t y f t y l +++++⎧=+⎪⎨=++=⎪⎩ (4.1.12) 对于已给的精度ε,当满足(1)()11 l l k k y y ε+++-<时,取(1)
11l k k y y +++=,然后继续下一
步2i y +的计算. 3.Runge-Kutta 方法
基本思想:利用Taylor 展式进行离散化.假设(4.1.9)中的(,)f t y 充分光滑,将
1()k y t +在k t 点作Taylor 展开:
2()
1()()'()''()()2!!
p p k k k k k h h y t y t hy t y t y t p +=+++++
其中
()()
'()(,())
''()[(,())]'()[(,())]t t y p p t y t f t y t y t f t y t f f f y t f t y t ===+⋅=
对照标准形式1(,;)k k k k y y h t y h +=+Φ.若取
1()
(,;)'()''()()2!!
p p h h t y h y t y t y t p -Φ=+++
并用k y 代替()k y t ,则得到一个p 阶近似公式
10
0(,;)
0,1,1()k k k k y y h t y h k n y y t +=+Φ⎧=-⎨=⎩ (4.1.13)
显然p=1时,式(4.1.13)就是(4.1.10),即为欧拉方法.但当2p ≥时,如果直接利用公式(4.1.13)求解数值解问题,就需要计算(,)f t y 的高阶微商.这个计算量是很大的.因此,直接利用式(4.1.13)构造高阶公式是不实用的.
通常,在求解高阶精度的数值解时,我们并不是直接使用Taylor 级数,而是利用它的思想,即计算(,)f t y 在不同结点的函数值,然后作这些函数值的线性组合,构造近似公式,式中有一些可供选择的参数.将近似公式与Taylor 展开式相比较,使前面的若干项密合,从而使近似公式达到一定的精度.
若线性组合选取的函数值个数为n ,计算达到的精度阶数为p ,这样的Runge-Kutta 方法就是我们通常所说的p 阶n 级Runge-Kutta 方法.
后来,德国学者Felhberg 对传统的Runge-Kutta 方法进行了改进,使原来的定步长算法变为自动变换步长的Runge-Kutta-Felhberg 方法,保证了更高的精度和数值稳定性,它已经成为求解数值解问题中最常见的方法.具体的算法我们在这里不做介绍,下面我们主要来讨论如何通过matlab 软件来运用这一方法求解微分方程的数值解问题.
二、用matlab 求解微分方程数值解 (一) 一阶微分方程组初值问题
Matlab 下求解一阶微分方程组初值问题(4.1.9)的数值解的最常见的方法
ode45()函数,该函数实现了前面介绍的变步长四阶五级Runge-Kutta-Felhberg 算法,可以采用变步长的算法求解微分方程. [t, y ] 中t 表示自变量,是标量.y 表示函数,可以是标量也可以是向量.0y 表示初值.当 y 为 n 维向量时,表明求解的问题是有n 个未知函数的方程组,则0y 也应为n 维向量.0[,]f t t 表示微分方程的求解区间,即自变量的取值范围.如果只给出
一个值 f t 则表示初始时刻为00t =.另外,该函数还允许0f t t >,即可以认为0t 为终值时刻,f t 为初始时刻,则0y 表示状态变量的终值,而该函数可以直接求解这样一个终值问题.
odefun 表示待解的微分方程(4.1.9)中(,)f t y 所写成的m-文件的名称,m-文件的编写格式为
也可以不把函数单独用m-文件保存,而是直接用inline()函数输入,从而获得函数名,编写格式为
其中,t 是自变量,即使所求方程是自治系统,也需要写出t ,否则matlab 在变量传递过程中将出现问题.y 为未知函数,ydot 表示未知函数的导数.
如果求解的问题是微分方程组问题,即我们之前讨论过的y 为向量的情况,则在书写上述m-文件或inline()函数时,待解方程(,)f t y 应以ydot 的分量形式写成,或者仍用单个方程的编写格式,多个方程用[ ]括起来,并以;隔开,在后面的例题中我们将具体的看到.
options 表示控制选项.在微分方程求解中有时需要对求解算法及控制条件进行进一步设置,这可以通过求解过程中的options 变量来进行修改.初始options 变量可以通过odeset()函数来获取,该变量是一个结构体变量,其中有众多成员变量,表4-1列出了常用的一些成员变量.
表4-1 微分方程求解函数的控制参数表
修改option 变量有两种方式: ①用odeset()函数设置,其格式为
达到的目的是,将相对误差上限设置为较小的710-
②直接修改options 的成员变量,也可以达到这样的目的,其格式为 12,,,m p p p 表示附加参数.即在函数表达式中,除了,t y 还有其他的可变的参数,引入这样的附加参数,使得微分方程的某些参数可以选择不同的值,对不同值求解时,考虑附加参数可以避免每次修改模型文件.但是值得注意的是,这种情况下的函数定义格式有特殊要求,格式为
或
注意,这里的flag 变量不能省略,用于占位. 例4-5 解微分方程组
33
'' 030(0)1,(0)0.5x x y y x y
t x y ⎧=--⎪=-<<⎨⎪==⎩
解: 该方程是非线性微分方程,不能用我们之前讨论的求解解析解的方法来处理,只能用数值解法求解.首先,我们将待解方程写成m -文件,用文件名eg1fun.m 保存.我们将2个未知变量x ,y ,以分量的形式写成一个向量变量x . %函数eg1fun.m function f=eg1fun(t,x) f(1)=-x(1)^3-x(2); f(2)=x(1)-x(2)^3;
f=f(:); %保证f 为列向量
这个m-文件还有另一种等价的书写方式.
%函数eg1fun.m的另一种书写方式
function f=eg1fun(t,x)
f=[-x(1)^3-x(2);x(1)-x(2)^3];
这样我们就完成了对待解方程的描述,然后输入命令
>>[t,x]=ode45('eg1fun',[0 30],[1;0.5]); %求解方程
subplot(1,2,1); %分2部分绘图,编辑第一部分
plot(t,x(:,1),t,x(:,2),':'); %在第一个图中绘制x,y的函数图象title('函数曲线图'); %图象名称
xlabel('t'); %横轴说明
legend('x(t)','y(t)'); %说明图例
axis square %显示正方形的坐标系
subplot(1,2,2); %编辑第二部分图象
plot(x(:,1),x(:,2)); %在第二个图中绘制相平面图
title('相平面图'); %图象名称
xlabel('x'); %横轴说明
ylabel('y'); %纵轴说明
axis square %显示正方形的坐标系
运行结果:
图4-1例4-5函数曲线和相平面图
运行程序,得到函数图象和相平面图.根据相平面上的相轨线,可以看出解从初值开始的变化趋势,这对于研究微分方程平衡点的稳定性是很有意义的,在后面的章节中我们将具体讨论. 例4-6 解微分方程组
⎪⎪
⎩⎪
⎪⎨
⎧
===-=-==1
)0(,1)0(,0)0(51.0'''321213312321y y y y y y y y y y y y 解:此问题仍为非线性微分方程组,我们求它的数值解,不妨将求解区间取为
[0,20]t =.下面我们用inline()函数来描述方程组,这样就不用单独编写m -文件来
描述方程.
输入命令:
>>eg2fun=inline('[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)]','t','y'); tf=20;
y0=[0;1;1];
[t,y]=ode45(eg2fun,[0,tf],y0);
plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+'); %绘出解的图象 legend('y1','y2','y3'); %说明图例 运行得到函数曲线图
图4-2 例4-6函数曲线图
例4-7 编写带有附加参数的matlab 函数,来描述食饵-捕食者模型
11112222211
2'()()'()()(0)(0)1
x t x a b x x t x a b x x x =-⎧⎪
=--⎨⎪==⎩
并对11223,2, 2.5,1a b a b ====和11222,1,1,0.5a b a b ====分别求解,画出解的曲线图和相轨线图.
解:选定附加参数为1122,,,a b a b ,编写如下m -文件来描述题述方程
%函数eg3fun.m
function xdot=eg3fun(t,x,flag,a1,b1,a2,b2) xdot=[x(1)*(a1-b1*x(2));-x(2)*(a2-b2*x(1))];
注意,这里的flag 变量不能省略.这样我们就完成了对方程的描述.我们不妨取自变量区间为[0,10],对方程求数值解.
输入命令 >>tf=10;x0=[1;1];
a1=3;b1=2;a2=2.5;b2=1;
[t,x]=ode45(‘eg3fun’,[0,tf],x0,[],a1,b1,a2,b2);
subplot(1,2,1);plot(t,x(:,1),t,x(:,2),':');title('函数曲线图'); xlabel('t');legend('x1','x2');axis square;
subplot(1,2,2);plot(x(:,1),x(:,2));title('相平面图'); xlabel('x1');ylabel('x2');axis square;
运行结果
图4-3 11223,2, 2.5,1a b a b ====
当11222,1,1,0.5a b a b ====时,只需要修改上述命令中的给a1,b1,a2,b2的赋值语句,就可以得到需要的结果.
图4-4
11222,1,1,0.5a b a b ====
(二) 一阶微分方程组边值问题
对于一阶常微分方程组的边值问题
'(,)
((),())0y f t y a t b g y a y b =⎧<<⎨=⎩ (4.1.14)
这里y ,f ,g 都可以是向量.
用其中odefun 是微分方程组函数,bcfun 为边值条件函数,sinit 是由bvpinit 得到的粗略解网络.求得的边值问题的解sol 是一个结构体变量,sol.x 为求解结点,sol.y 为数计算由bvp4c 得到的解在ti 的值. 例4-8 求解边值问题
1221
11
'' ()0,()20y y y y a t b y a y b =⎧⎪
=-<<⎨⎪=+=⎩ 解:输入命令
>>sinit=bvpinit(0:4,[1;0]);
odefun=inline('[y(2);-abs(y(1))]','t','y'); bcfun=inline('[ya(1);yb(1)+2]','ya','yb'); sol=bvp4c(odefun,bcfun,sinit); t=linspace(0,4,101); y=deval(sol,t);
plot(t,y(1,:),sol.x,sol.y(1,:),'o',sinit.x,sinit.y(1,:),'x'); legend('解曲线','求解点','粗略解');
运行结果
图4-5 例4-8函数曲线图
(三) 高阶微分方程
由于matlab 只能处理一阶微分方程组问题,对于高阶微分方程的初边值问题,我们必须通过变量代换,使它转换为一阶方程组,才能用matlab 求解数值解.我们考虑高阶常微分方程的初值问题
()(1)(,,',,)n n y f t y y y -= (4.1.15) 且初始值(1)(0),'(0),(0)n y y y - 为已知.
令(1)12,',,n n x y x y x y -=== ,通过这样的变量代换,高阶微分方程(4.1.15)就转化为一阶方程组的形式
1223
12'' '(,,,,)
n n x x x x x f t x x x =⎧⎪=⎪⎨
⎪⎪=⎩ (4.1.16) 且初值条件为(1)12(0)(0),(0)'(0),,(0)(0)n n x y x y x y -=== .这样,我们就可以通过之前讨论过的一阶微分方程组的数值解法来研究这一问题.
对于高阶微分方程组,只要依次将每个未知变量的各阶导数均单独定义变量,
就可以同样的转化为一阶微分方程组. 例4-9 用数值解法求解微分方程
2''(1)'0
(0)0.2,'(0)0.7
y y y y y y ⎧+-+=⎨
=-=-⎩ 解:首先我们需要将这个2阶微分方程转化为一阶微分方程组,我们令
1x y =,2'x y =,则原方程转换为
122
2121
1
2''(1)(0)0.2,(0)0.7
x x x x x x x x =⎧⎪=---⎨⎪=-=-⎩ 对[0,20]区间求数值解,输入命令
>>tf=20;x0=[-0.2;-0.7];
eg5fun=inline('[x(2);-((x(1))^2-1)*x(2)-x(1)]','t','x'); [t,y]=ode45(eg5fun,[0,tf],x0); plot(t,y);
legend('y','y'''); 运行得到函数图象
图4-6 例4-9函数曲线图
例4-10 求解微分方程组(竖直加热板的自然对流问题)
23232
22
223202.10(0)0,(0)0,(0)0.68,(0)1,(0)0.5d f d f
df f T d d d d T dT f d d df d f dT
f T d d d η
ηηη
ηηηη⎧⎛⎫+-+=⎪ ⎪⎝⎭⎪⎪⎪
+=⎨⎪⎪=====-⎪⎪⎩
解:这是一个高阶方程组问题,我们需要通过变量代换转化为一阶微分方程组,从
而寻求数值解.我们令
2123452
,,,,df d f dT
y f y y y T y d d d ηηη
===== 原方程转化为如下一阶方程组
1
2232
3132454515
1
2345,32, 2.1(0)0,(0)0,(0)0.68,(0)1,(0)0.5
dy dy y y d d dy y y y y d dy
dy y y y d d y y y y y ηηη
η
η⎧==⎪⎪⎪=-+-⎪
⎨⎪==-⎪⎪⎪=====-⎩ 输入命令
>>y0=[0;0;0.68;1;-0.5];tf=5;
eg6fun=inline('[y(2);y(3);-3*y(1)*y(3)+2*y(2)^2-y(4);y(5);-2.1*y(1)*y(5)]','t','y'); [t,y]=ode45(eg6fun,[0,tf],y0); plot(t,y(:,1),t,y(:,4),':'); legend('f','T'); 运行得到函数图象
图4-7 例4-10函数曲线图
§4.1.3导弹追踪问题
一、导弹追击问题
一艘敌舰在某海域内沿正北方向航行时,我方战舰恰位于敌舰的正西方向1km 处.我舰向敌舰发射导弹,导弹头始终对准敌舰,敌舰速度为0v =1km/min,导弹速度为敌舰速度的5倍.问需要多长时间,在何处导弹击中敌舰?
以我舰位置为坐标原点,正北方向为y 轴方向建立坐标系,设t 时刻导弹所处的位置为((),())P x t y t ,敌舰所处的位置为0(1,)Q v t .
图4-8 导弹追击问题
由于导弹头始终对准敌舰,因此直线PQ 是导弹运行轨迹OP 在P 点的切线,
即
01v t y dy dx x
-=- 从而
0(1)
dy
v t x y dx
=-+ (4.1.17) 又因为导弹速度为敌舰速度的5倍,所以OP 弧的长度为AQ 长度的5倍. 即
00
5v t =⎰
(4.1.18) 联立(4.1.17)和(4.1.18)得到
5(1)5dy x y dx =-+⎰
两边对x 求导,得到微分方程
2
25(1)d y
x d x
=-
(4.1.19) 该问题的初始条件为
(0)0,'(0)0y y ==
首先,我们将二阶方程化为一阶方程组,令12,dy
y y y dx
==
,(4.1.19)转化为 ()1
21/2
2221215(1)(0)0,(0)0dy y dx y dy
dx
x y y ⎧=⎪⎪
⎪+⎪=
⎨-⎪⎪==⎪⎪⎩ (4.1.20) 对于这一问题,我们可以求得方程的解析解,从而确定导弹击中敌舰位置和时
间,也可以通过求数值解的方法得到同样的结论. 1) 解析解.
关于2y 的方程
()1/2222215(1)(0)0
y dy dx x y ⎧+⎪=
⎨-⎪
=⎩ 分离变量,有
()
2
1/2
225(1)
1dy dx
x y =
-+
两边积分,并代入初始条件,得
21ln(ln(1)5
y x +=--
即
15
2(1)y x -
=- (4.1.21) 又因为
15
2(1)y x ==-- (4.1.22)
(4.1.21)和(4.1.22)两式相加,得到
115
521((1)(1))2y x x -=---
下面,我们求解关于1y 的方程
1115
521
1((1)(1))2(0)0dy y x x dx
y -⎧==---⎪
⎨⎪=⎩ 直接积分,并代入初始条件得到
46
551555
(1)(1)81224
y y x x ==--+-+
当1x =时,524y =
,即当敌舰航行到点(1,524
)处时,被导弹击中.被击中的时间为0
12.5y
t s v =
=. 2) 数值解.
首先,编写出描述方程(4.1.20)的m-文件 function dy=equ20(x,y)
dy=[y(2);1/5*sqrt(1+y(2)^2)/(1-x)];
由于方程在1x =处没有定义,我们取求解区间为[0,1-1e-8],输入命令 xf=1-1e-8;
[x,y]=ode45('equ20',[0,xf],[0;0]) plot(x,y(:,1));
hold on y=0:0.01:2; plot(1,y,'*');
运行结果 y =
0 0
0.0000 0.0001 …… …… 0.2083 16.1097 0.2083
20.0688
图4-9 导弹击中曲线
说明当1x =时,0.2083y =,导弹击中敌舰的位置大致在(1,0.2083)处,击中的时间0
y
t v =
=0.2083min=12.498s.这与求解析解得到的结果是一致的. 比较上面2种求解方法,我们发现,解析解可以得到解的真实值,但是方法比较复杂,技巧性强,只能处理比较简单的方程问题;而数值解虽然得到的解与真实值之间会有一些允许范围内的误差,但是方法比较容易掌握,更具一般性. 二、导弹系统的改进
现根据情报,这种敌舰能在我舰发射导弹后T 分钟做出反应并摧毁导弹.因此,我们要求改进电子导弹系统,使其根据敌舰与我舰的距离,行使方向和速度,能自动判断出敌舰是否在有效打击范围之内.
对于更一般的情况,我们做出如下假设,设敌舰在我舰正东方向d km 处,行驶
速度为0v km/min,行驶方向与正东方向的夹角为θ,导弹的飞行速度为v km/min.问题的关键是计算出导弹击中敌舰所需要的时间*t ,并将*t 与T 比较,若*t <T ,则敌舰在打击范围内.
我们仍以我舰位置为坐标原点,以正北方向为y 轴建立坐标系,设t 时刻导弹所处的位置为((),())P x t y t ,敌舰所处位置为00(cos ,sin )Q d v t v t θθ+
.
图4-10 导弹击中问题
由于导弹头始终对准敌舰,因此直线PQ 是导弹运行轨迹OP 在P 点的切线,即
00sin cos v t y dy
dx d v t x
θθ-=
+- (4.1.23) 又因为导弹速度为敌舰速度的
0v v 倍,所以OP 弧的长度为AQ 长度的0
v
v 倍.即
v x v t v =⎰
(4.1.24) 联立方程(4.1.23)和(4.1.24),消去0v t ,再对方程两边对x 求导,得到微分方程
21/2
202
2
(1())(s i n c o s )
()(s i n c o s )c o s (())v d y d y d y v dx dx dy dy dx d x d x y dx dx
θθθθθ+-=--+-+ (4.1.25) 同样,我们令12,dy
y y y dx
==,(4.1.25)转化为一阶微分方程组
1
221/22
022222112(1())(sin cos )()(sin cos )cos (())(0)0,(0)0dy y dx v y y dy v dx d x y y d x y y y θθθθθ⎧=⎪⎪
⎪+-⎪=⎨--+-+⎪
⎪==⎪⎪⎩
(4.1.26) 从图象上来判断,我们知道,当P 和Q 两点的运动曲线相遇时,导弹击中敌舰.因此我们可以认为,若**,x y 满足方程(4.1.25)且
**()tan y x d θ=- (4.1.27) 则点(**,x y )为导弹击中敌舰的击中点.再根据Q 点的表达式
00(cos ,sin )Q d v t v t θθ+,可以计算出击中时间*
*
0sin y t v θ
=
.若*t <T ,则敌舰在打击范围内,可以发射.
从上述分析,我们知道,当0v 、v 、d 、θ、T 已知时,根据(4.1.25)和(4.1.27)可以计算出导弹击中敌舰所需要的时间*t ,且当*t <T 时,敌舰在打击范围内.这种对于每组0v 、v 、d 、θ、T 分别计算从而作出判断的方法,称为在线算法.如果在敌舰和我舰的情况已知,即0v 、v 、T 已知,计算出所有在打击范围内的d 和θ,之后要判断敌舰是否在打击范围,只需要在计算结果中查询,这种算法称为离线算法.从原理上来看,我们发现,在线算法灵活,容易调整参数和模型,但是速度比较慢;离线算法事先计算好,实时使用查询方式,不需要计算,速度快. (1)在线算法
现已知0v =90km/h,v=450km/h,d=50km,4
π
θ=
,T=0.1h.首先,我们编写带参数的m-文件描述方程(4.1.26).为了使方程始终有意义,我们通常在分母上加上一个无穷小量1e-8.
%方程(4.1.26)
function dy=equ26(x,y,v0,v,d,theta) dy(1)=y(2);
dy(2)=((v0/v)*(1+(y(2))^2)^(1/2)*(sin(theta)-y(2)*cos(theta))^2)
/((d-x)*(sin(theta)-y(2)*cos(theta))+cos(theta)*(y(2)*(d-x)+y(1))+1e-8); dy=dy(:);
输入命令
clear;close;
v0=90;v=450;d=50;theta=pi/4;T=0.1;
%初始化参数
[x,y]=ode45(@equ26,[0,60],[0;0],[],v0,v,d,theta);
%计算导弹运行轨迹方程的数值解
z=(x-d)*tan(theta);
%计算敌舰运行的轨迹
n=length(x);
for i=1:n
if abs(z(i)-y(i,1))<1e-6
xk=x(i);
yk=y(i,1);
break;
end
end
xk
yk
%求出击中点xk,yk坐标
for i=1:n
if z(i)>0
z1(i)=z(i);
else
z1(i)=0;
end
end
plot(x,y(:,1),x,z1(:),xk,yk,'*');
legend('导弹运行轨迹','敌舰运行轨迹','击中点');
tk=yk/(v0*sin(theta)+1e-8)
%计算击中时间
if tk<T
disp(['敌舰在打击范围内，击中地点在',
num2str(xk),num2str(yk),'击中时间为', num2str(tk)]);
else
disp(['敌舰不在打击范围内']);
end
得到结果
击中点xk=58.4056,yk=8.4056,击中时间tk =0.1321,所以敌舰不在打击范围内,应等接近一些再发射.
图4-11 在线算法
(2)离线算法
在介绍离线算法之前,我们先来看方程(4.1.26)的另一种表达.因为导弹的线速度
v = (4.1.28) 方程(4.1.28)与(4.1.23)联立,可以建立以t 为参数的关于x ,y 的参数方程
dx
dt dy dt
⎧==⎪⎪⎪⎪
⎨⎪==
⎪⎪⎪⎩ (4.1.29)
当()x t 满足0()cos x t d v t θ≥+,则导弹击中了敌舰.由于我们所要求的打击范围是在t T <中讨论的,所以考虑以t 为参数的方程(4.1.29)能使求解过程大大简化.
现已知0v =90km/h,v=450km/h,T=0.1h,要计算出所有在打击范围内的d 和θ.依题意,0θπ≤≤.据此,我们可以确定d 的取值范围.当0θ=时,敌舰正好背向行驶,
导弹直线运行,击中时间
0/()t d v v T =-<
求得min 0()36km d T v v =-=.当θπ=时,敌舰迎面驶来,导弹直线运行,击中时间
0/()t d v v T =+<
则max 0()54km d T v v =+=.所以3654d ≤≤.这样,我们对于所有可能的d 和θ的取值,计算击中所需时间,从而对不同的θ,得到d 的临界值.具体应用时直接查询判断.
编写如下m -文件描述方程(4.1.29) %方程(4.1.29)
function dy=equ29(t,y,v0,v,d,theta)
dydx=(v0*t*sin(theta)-y(2))/(d+v0*t*cos(theta)-y(1)+1e-8); dy(1)=v/(1+dydx^2)^(1/2); dy(2)=v/(1+dydx^(-2))^(1/2); dy=dy(:);
输入命令 clear;close;
v0=90;v=450;d=50;theta=pi/4;T=0.1; i=1;
for d=54:-1:36
for theta=0:0.1:pi
[t,y]=ode45(@equ29,[0,T],[0;0],[],v0,v,d,theta); if max(y(:,1)-d-v0*t*cos(theta))>0 range(i,:)=[d,theta]; i=i+1; break; end end end figure;
plot(range(:,1),range(:,2)); xlabel('d');ylabel('theta');
运行得到临界曲线,即在该曲线上方的d 和θ值,所对应的敌舰位置在打击范围内,曲线下方不在打击范围内.
图4-12 离线算法
因此d=50km,4
π
θ=
,不在打击范围内.
§4.1.4微分方程稳定性理论简介
在处理实际问题时,对于有些微分方程模型我们不仅要得到问题的解,有时还需要研究解的稳定性,即解对初始值的连续依赖性,如果解在一定范围内是稳定的,那么初始条件发生一些小的扰动(如实验测量误差等),对问题的解不会造成影响.另外,还有一些问题我们并不需要求解,而通过解的变化趋势的研究,并分析一些特殊解的稳定性就可以解决问题.本节仅介绍几类特殊,但常用的常微分方程稳定性分析的基本方法. 一、单个常微分方程的平衡点及稳定性
若微分方程
)(x f dt
dx
= (4.1.30) 方程右端不显含自变量t ,称为自治方程.代数方程的实根0x x =称为方程(4.1.30)的平衡点(或奇点).注意到,平衡点也是方程的解(奇解).
如果从一定范围内的初始条件出发,方程(4.1.30)的解)(t x 都满足
,)(lim 0x t x t =+∞
→
则称平衡点0x 是稳定的.
对于一些不易求解的问题,我们可以不求方程(4.1.30)的解,不用定义来判断平衡点0x 的稳定性,下面我们来介绍这种方法.
将)(x f 在0x 点作泰勒(Taylor )展开,只取一次项,得到方程(4.1.30)的近似线性方程
),)(('00x x x f dt
dx
-= (4.1.31) 0x 也是方程(4.1.31)的平衡点,方程(4.1.31)的通解为0|)('|0)(x ce t x t x f +=.关于0x 点稳定性有如下结论:
① 若0)('0<x f ,则0x 对于方程(4.1.31)和(4.1.30)都是稳定的; ② 若0)('0>x f ,则0x 对于方程(4.1.31)和(4.1.30)都是不稳定的. 二、二阶常微分方程组的平衡点及稳定性
现在讨论二阶微分方程组
(;,),(;,).dx
f t x y dt
dy g t x y dt
⎧=⎪⎪⎨⎪=⎪⎩ (4.1.32)
它的解
(),()x x t y y t ==
在以,,t x y 为坐标的欧氏空间中决定了一条曲线,如果把时间t 看作参数,仅考虑
,x y 为坐标的欧氏空间,此空间称为方程(4.1.32)的相平面(若方程组是高阶,则称
为相空间).对于右端函数不显含时间t 的自治系统 ⎪⎩⎪⎨⎧==).,(),
,(y x g dt
dy
y x f dt dx
(4.1.33)
代数方程组
⎩
⎨
⎧==.0),(,
0),(y x g y x f 的实根00,y y x x ==称为方程(4.1.33)的平衡点,记作),(000y x P .它也是方程(4.1.33)的解.
如果从一定的范围内的初始条件出发,方程(4.1.33)的解)(),(t y t x 都满足
,)(lim ,)(lim 00y t y x t x t t ==+∞
→+∞
→
则称平衡点0P 是稳定的,否则称0P 是不稳定的.
与单个方程的讨论类似,对于不易求解的微分方程组,我们可以通过研究与
其对应的近似线性方程组,从而得到平衡点和稳定性.在这里,我们省略证明过程,只给出判别平衡点0P 是否稳定的判别准则.令
,)()()()
( ,)()(000000y
P g x
P g y
P f x P f q y P g x P f p ∂∂∂∂∂∂∂∂=
⎥⎦
⎤⎢⎣⎡∂∂+∂∂-= 关于0P 点稳定性有如下结论:
① 当0>p 且0>q 时,平衡点0P 是稳定的; ② 当0<p 或0<q 时,平衡点0P 是不稳定的.。