微分方程的稳定性

⇒ 乙方胜
m=0
mc
m<0
y 0 > d = rx s rx s x 线性律 x 0 c ry s ry s y 模型
m < 0 ⇒甲方胜
0
−m d
x(t)
m = 0 ⇒ 平局
混合战争模型 甲方为游击部队，乙方为正规部队
⎧ x = −cxy
⎪ ⎨
y
=
− bx
⎪⎩ x(0)
=
x, 0
y(0)
=
y 0
然后继续下一步yi+2 的计算。
此即改进的欧拉法。
3、使用泰勒公式
以此方法为基础，有龙格-库塔法、线性多步法等方法。
4、数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）
时（k为正整数，h为步长），称它是一个k阶公式。 k越大，则数值公式的精度越高。
•欧拉法是一阶公式，改进的欧拉法是二阶公式。 •龙格-库塔法有二阶公式和四阶公式。 •线性多步法有四阶阿达姆斯外插公式和内插公式。
y (t )
n > 0,乙方胜
n = 0,平局 n < 0,甲方胜
cy 2 − 2bx = n
n
=
cy
2 0
−
2bx 0
n>0 乙方胜
⎜⎜⎝⎛
y0 x0
⎟⎟⎠⎞ 2
>
2b cx 0
⎛⎜⎜⎝
y 0
x 0
⎞⎟⎟⎠ 2
>
2r p s x xx rs x y ry 0
设 x0=100, rx/ry=1/2, px=0.1, sx=1(km2), sry=1(m2)
• 每方战斗减员率取决于双方的兵力和战斗力 • 每方非战斗减员率与本方兵力成正比 • 甲乙双方的增援率为u(t), v(t)
⎧ x(t) = f ( x, y) − αx + u(t), α > 0
模型
⎨ ⎩
y (t)
=
g
(x,
y)
−
βy
+
v (t ),
β >0
f, g 取决于战争类型
正规战争模型 双方均以正规部队作战
∫ y(xi+1 ) − y(xi ) =
xi +1 xi
f
(t,
y(t ))dt
≈
xi+1 − 2
xi
[f
(xi , y(xi )) +
f
(xi+1 ,
y(xi+1 ))]
故有公式：
⎪⎧ ⎨
y
i
+1
=
yi
+
h 2
[
f
(
xi
,
yi
)
+
f (xi+1 , yi+1 )]
⎪⎩ y0 = y(x0 )
实际应用时，与欧拉公式结合使用：
• 甲方战斗减员率只取决于乙方的兵力和战斗力
f(x, y)=−ay, a ~ 乙方每个士兵的杀伤率
a=ry py, ry ~射击率， py ~命中率
⎧ x = − ay − αx + u (t )
⎨ ⎩
y
=
−bx
−
βy
+
v (t )
g = −bx, b = rx px
⎧ x = − ay
• 忽略非战斗减员
y) ，其数值解是指由初始点x y0
0
开始
的若干离散的x值处，即对x0 < x1 < x2 < " < xn，求出准确值y(x1 ),
y(
x2
),",
y(
x
n
)
的相应近似值y 1
,
y
2
,",
y
。
n
（二）建立数值解法的一些途径
设 x i+1 − xi = h, i = 0,1,2,"n −1, 可用以下离散化方法求解微分方程： ⎧y'= f(x, y) ⎩⎨y(x0 ) = y0
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0
2
4
6
8
10
12
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]);
plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
3、结果如图
图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.
⎪ ⎨
y
=
− bx
• 假设没有增援
⎪⎩ x ( 0 ) = x 0 , y ( 0 ) = y 0
正规战争模型
⎧ x = − ay
⎪ ⎨
y
=
−bx
⎪⎩ x(0) = x0 , y (0) = y0
为判断战争的结局，不求x(t), y(t) 而在相平面上讨论 x 与 y 的关系
dy = bx dx ay
dx2 例 1 求 du = 1+ u 2 的通解.
dt 解 输入命令：dsolve('Du=1+u^2','t')
结果：u = tg(t-c)
例 2 求微分方程的特解.
⎪⎧d 2 ⎨ dx
y
2
+
4
dy dx
+
29
y
=
0
⎪⎩ y(0) = 0, y' (0) = 15
解 : 输入命令:
y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
4.4 稳定性问题
在研究许多实际问题时，人们最为关心的 也许并非系统与时间有关的变化状态，而是系 统最终的发展趋势。例如，在研究某频危种群 时，虽然我们也想了解它当前或今后的数量， 但我们更为关心的却是它最终是否会绝灭，用 什么办法可以拯救这一种群，使之免于绝种等 等问题。要解决这类问题，需要用到微分方程 或微分方程组的稳定性理论。在下两节，我们 将研究几个与稳定性有关的问题。
4.1 正规战与游击战
第一次世界大战Lanchester提出预测战 役结局的模型.
• 战争分类：正规战争，游击战争，混合战争 • 只考虑双方兵力多少和战斗力强弱 • 兵力因战斗及非战斗减员而减少，因增援而增加 • 战斗力与射击次数及命中率有关
一般模型 x(t) ~甲方兵力，y(t) ~乙方兵力
模型 假设
稳定性模型
• 对象仍是动态过程，而建模目的是研究时 间充分长以后过程的变化趋势 ——平衡状 态是否稳定。
• 不求解微分方程，而是用微分方程稳定性 理论研究平衡状态的稳定性。
一般的微分方程或微分方程组可以写成： dx = f (t, x) dt
定义 称微分方程或微分方程组
dx = f (x) dt
（3.28）
为自治系统或动力系统。
若方程或方程组f (x)=0有解xo,x = xo显然满足(3.28)。 称点xo为微分方程或微分方程组(3.28)的平衡点或奇点。
（三）用Matlab软件求常微分方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 函数 量值 值
ode45 ode23 ode113 ode15s ode23s
由待解 方程写 成的m-
ts=[t0，tf]， t0、tf为自
函数的 初值
变量的初
文件名 值和终值
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法
y=simple(y)
z=simple(z)
结果为：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t
2
2013/2/21
4.3 微分方程的数值解
ay 2 − bx 2 = k k = ay02 − bx02
y (t )
k >0
k > 0 ⇒ x = 0时y > 0 乙方胜
k =0
ka
k <0
⎜⎜⎝⎛
y0 x0
⎟⎟⎠⎞ 2
>
b a
=
rx ry
px py
平方律 模型
k < 0 ⇒甲方胜
0
−k b
x(t)
k = 0 ⇒平局
1
2013/2/21
游击战争模型 双方都用游击部队作战
⎪⎧
y (0) i+1
=
yi
+ hf (xi , yi )
⎨ ⎪⎩
y (k +1) i+1
=
yi
+
h 2
[
f
(
xi百度文库
,
yi
)
+
f
(
xi+1
,
y (k ) i +1
)]
k = 0,1,2,"
对于已给的精确度 ε，当满足
y (k +1) i +1
−
y(k) i +1
< ε 时，取 yi+1
=
y ， （k +1） i +1
1、用差商代替导数 若步长h较小，则有
y'(x) ≈ y(x + h) − y(x) h
故有公式：
⎧ ⎨ ⎩
y y
i+1 = yi + 0 = y(x0 )
hf
(
xi
,
yi
)
i = 0,1,2,", n -1
此即欧拉法。
2、使用数值积分
对方程y’=f(x,y), 两边由xi到xi+1积分，并利用梯形公式，有：
dy(2)=1000*(1-y(1)^2)*y(2)-y(1);
2、取t0=0，tf=3000，输入命令： [T,Y]=ode15s('vdp1000',[0 3000],[2 0]);
plot(T,Y(:,1),'-')
3、结果如图
例 5 解微分方程组.
⎧
y1' = y2 y3
⎪⎪ ⎨ ⎪
y2 ' = − y1 y3 y3' = −0.51y1 y2
⎪⎩ y1(0) = 0, y2 (0) = 1, y3(0) = 1
解 1、建立m-文件rigid.m如下：
function dy=rigid(t,y)
dy=zeros(3,1);
dy(1)=y(2)*y(3);
dy(2)=-y(1)*y(3);
dy(3)=-0.51*y(1)*y(2);
1
0.8
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
注意:
1、在解n个未知函数的方程组时，x0和x均为n维向量， m-文件中的待解方程组应以x的分量形式写成.
（一）常微分方程数值解的定义
在生产和科研中所处理的微分方程往往很复杂且大多 得不出一般解。而在实际上对初值问题，一般是要求得 到解在若干个点上满足规定精确度的近似值，或者得到 一个满足精确度要求的便于计算的表达式。
因此，研究常微分方程的数值解法是十分必要的。
对常微分方程
：⎨⎧ ⎩
y' = y(x
f(x, 0) =
⎧ x = − cxy
⎪ ⎨
y
=
− dxy
⎪⎩ x (0) = x0 , y (0) = y 0
游击战争模型
⎧ x = − cxy
⎪ ⎨
y
=
− dxy
⎪⎩ x (0)
=
x, 0
y(0)
=
y 0
y (t )
m>0
dy = d dx c
cy − dx = m m = cy 0 − dx 0
m > 0 ⇒ x = 0时y > 0
结果为 : y =3e-2xsin（5x）
例 3 求微分方程组的通解.
⎧ ⎪ ⎪⎪ ⎨ ⎪ ⎪ ⎪⎩
dx
dt dy
dt dz
dt
= = =
2x 4x 4x
− − −
3y 5y 4y
+ + +
3z 3z 2z
解 输入命令 ：
[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', 't')； x=simple(x) % 将x化简
2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组.
3
2013/2/21
例4
⎪⎧ ⎨
d2x dt 2
−1000(1
−
x2
)
dx dt
−
x
=
0
⎪⎩ x(0) = 2; x'(0) = 0
解: 令 y1=x，y2=y1’
则微分方程变为一阶微分方程组：
⎧ ⎪ ⎨
y2
'
=
y1' = y2 1000(1− y12
数学建模讲义
主讲人:穆学文
西安电子科技大学数学系 Email:mxw1334@126.com
2013/2/21
第四讲 微分方程模型
-------多种群增长模型
4.1 正规战与游击战 4.2 微分方程解析解 4.3 微分方程数值解 4.4 微分方程稳定性 4.5 捕鱼业的持续收获 4.6 军备竞赛 4.7 种群的相互竞争 4.8 种群的相互依存 4.9 种群的弱肉强食
( y0 / x0 )2 > 100
0
x(t) 乙方必须10倍于甲方的兵力
4.2 微分方程的解析解
求微分方程（组）的解析解命令: dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
记号: 在表达微分方程时，用字母 D 表示求微分，D2、D3 等 表示求高阶微分.任何 D 后所跟的字母为因变量，自变量可以指 定或由系统规则选定为确省. 例如，微分方程 d 2 y = 0 应表达为：D2y=0.
• 甲方战斗减员率还随着甲方兵力的增加而增加 f(x, y)=−cxy, c~ 乙方每个士兵的杀伤率
c = ry py ry~射击率 py ~命中率
py=sry /sx sx ~ 甲方活动面积 sry ~ 乙方射击有效面积
g(x, y) = −dxy, d = rx px = rxsrx / sy
• 忽略非战斗减员 • 假设没有增援
)
y2
−
y1
⎪⎩ y1(0) = 2, y2 (0) = 0
1、建立m-文件vdp1000.m如下：
function dy=vdp1000(t,y)
dy=zeros(2,1);
dy(1)=y(2);
2 1.5
1 0.5
0 -0.5
-1 -1.5
-2 -2.5
0
500
1000
1500
2000
2500
3000