预测战争模型

案例 预测战争模型
1.1问题描述 (1)
1.2分析与建模 (1)
1.3模型求解 (2)
1.4模拟求解 (5)
1.4.1运行情况1（a=0.4,b=0.10，delta_t = 0.3） (5)
1.4.2运行情况2，a=0.15; b=0.1;delta_t=0.05; (7)
1.4.3运行情况3，a=0.15; b= 0.1;delta_t=0.001; (8)
1.4.4求解程序 (10)
1.1 问题描述
在第一次世界大战期间，nchester提出了预测战争结局的数学模型。

根据战争的不同特性，他给出了三种作战模型。

在建立模型时，简化了许多因素，模型变得简单，但仍有一定的实际意义。

现考虑X、Y两方孤立交战的部队，双方均无增援部队的情况。

希望为这场战斗建立一个数学模型，应用这个模型达到如下目的：
1.预测哪一方将获胜？
2.估计获胜的一方最后剩下多少士兵？
3.计算失败的一方开始时必须投入多少士兵才能赢得这场战斗？
4.战斗的持续时间。

1.2 分析与建模
假定X部队t时刻存活的士兵数为x(t)，而Y部队在t时刻存活的士兵数为y(t)，将x(t)与y(t)都看作连续变量。

并假定双方所有士兵不是战死就是活着参加战斗，即不考虑俘虏和伤员。

关于双方的作战伤亡情况，一种合理的假设是，在Δt时间内，X部队被杀死的士兵数Δx将取决于Δt的长短，以及在Δt起始时刻与其交战的Y部队的士兵数。

假定是一种正比关系，即
Δx=－ayΔt
其中a是一个常数，，代表了Y部队的战斗力，称为“杀伤率”，更明确地说，a 是Y部队的一个士兵在单位时间内杀死X部队的名士兵数。

类似地，对于Y部队有
Δy－axΔt
令Δt→0，得到两个微分方程
,ay dt dx
−= ( a>0) ,bx dt
dy
−= ( b>0) 从而得到联立微分方程组如下：
⎪⎪⎪⎩
⎪⎪⎪⎨⎧==>−=>−=.
)0(,)0()0(,,
)0(,00y y x x b bx dt
dy
a ay dt dx
(6.3.4) 1.3 模型求解
对联立微分方程组(6.3.4)中的任一方程进行积分，直接求出方程的解是很困难的也无必要。

根据作战的实际背景，可以分析出以下几点：方程组里的变量满足x ≥0，y ≥0，有唯一平衡点（0，0）； x (t)，y (t)都是递减函数，且随着x ，y 的减小，其衰减速度也在降低。

在我们的模型中，若有一方部队的士兵数为零，就标志着战斗的结束。

将两个方程相除消去时间变量t ，得
ay
bx dx dy = 可分离变量 bxdx aydy =
对两边积分得到
2
2222c
bx ay += 或者
c bx ay +=2
2
代入初始条件，有
c bx ay +=2
020
（6.3.5） )()(20220
2x x b y y a −=−在相平面（xy 平面）上，轨线是双曲线的一部分，如图6.3.1所示。

为预测何方部队获胜，将剩下多少士兵，先考虑一种特殊情况。

战斗开始时双方投入兵力满足 ，解曲线方程（6.3.5）化为
2
02
0bx ay = 或 22bx ay =x a
b y =
方程的轨线是一条过原点，斜率为a b /的直线，称为等战斗力直线。

这种情
况下，战斗将无情地从（x 0，y 0）点进行到（0，0）点，表明这是一场势均力敌的，导致相互毁灭的战争。

若，可从相位图观察到点（x 2020bx ay >0，y 0）位于等战斗力直线的左上方，可以判断Y 部队将获胜，在轨线方程中令y=0可以验证这个结论，若令x=0 , 得
a bx ay y /)(2
020−=，
即为根据模型预测出的Y 部队获胜时的幸存士兵数。

如果战败的一方希望转败为胜，那么开始他们应投入多少兵力呢？ 现假设双方初始兵力分别为x 0=10000，y 0=5000，，假设在一个小时（单位时间）内每个Y 部队的士兵杀死0.15个X 部队士兵，而X 部队的每个士兵在一个小时内杀死0.1个Y 部队士兵，则方程组（6.3.4）为
⎪⎪⎪⎩
⎪⎪⎪⎨⎧====−=−=,5000)0(,10000)0(,
1.0,15.000y y x x x dt
dy
y dt dx
（6.3.6） 由于，，有 ，模型预测X 部队将获胜，Y 部队若要获胜最初投入兵力y 62201075.3)5000(15.0×=×=ay 622
1010)10000(1.0×=×=bx 02
020>−ay bx 0必须满
足2020x a b y >，即应满足：822020103
11000015.01.0×=×=>x a b y 或 y 0>8165。

怎样估算战斗的持续时间不太容易，现在先不去求解方程组，用分析方法做
出估计。

有
10001.000−=−==x dt dy
t ， 意味着战斗开始时Y 部队的士兵以每小时1000人的速度被歼灭，如果一直持续这种速度，仔细思考实际情不会如此，因为X 部队的士兵数也在减少，故战斗至少持续5000/1000=5（小时）。

战斗结束时X 部队余下的士兵数为
79061075.31010/)(662020≈×−×=−b ay bx 名。

此时，Y 军队士兵被歼灭的速度为 1.79079011.0−=×−≈结束
dt dy
这是Y 部队士兵被歼的最慢速度，若保持不变，有
y =－790.1t + 5000,
令y=0, 解得t=5000/790.1≈6.32小时，应为Y 部队被歼灭的最长时间。

分析结果表明，战斗会持续5～6.32小时，取中间值约为5.7小时。

通过求解微分方程组可得到确切的答案。

将方程组（6.3.6）的第一个方程
两边微分，得
⎪⎪⎩⎪⎪⎨⎧−=−=.1.0,15.02
2x dt
dy
dt dy
dt x d 把其中的第二个方程代入第一方程，有
0015.02=−x dt
x
d 二次积分得到方程的解为
t t Be Ae x 015.0015.0−+=，t ≥0， （6.3.7）
其中A 和B 是积分常数，代入初始条件x 0=10000，y 0=5000，并令t=0,有
A+B=10000 （6.3.8）
对解函数（6.3.7）两边求导
][015.0015.0015.0t t Be Ae dt
dx
−−=，t ≥0 因 dt dx
y 15.01−= 0
015.01=−=t dt dx y
)(15
.0015
.050000B A y −−
== （6.3.9） 联立求解（6.3.8）和（6.39），解得A ≈1938.14，B ≈8061.86。

在任意时刻Y 军队的士兵数为
dt dx
t y 15.01)(−=
]86.806114.1938[15
.0015
.0015.0015.0t t e e −−−=
，t ≥0， 令y=0 ，算得16.4015
.02≈=a
b e
t ，解得t ≈5.82小时，证明我们前面估算的战斗持续时间相当准确。

前面已得出Y 部队要赢得这场战斗，开始时必须再增加3165名士兵。

假定战斗开始后的某个时刻到达增援部队（如空降伞兵），设在任一时刻使Y 部队战斗力发生改变从而赢得胜利所需的增援人数为N ，则 y x y x a b
N −=−≥
015
.01
.0 可以算出以下结果 战斗持续时间t （小
时）
0 1 2 3 4 5
增援士兵数N(个)
3165 3577 4043 4570 5166 5839
1.4 模拟求解
微分方程
)(t f dt
dy
=
则可得
dt t f dy )(=
则利用微分一阶近似计算（即泰勒一阶展式）：
)()()(t y t t y t t f dy −Δ+≈Δ=
t t f t y t t y Δ+≈Δ+)()()(
特别的，若令，
t t t i i Δ=−+1则
t t f t y t y i i i Δ+≈+)()()(1
要讨论的问题：
1. 这个模型怎么理解？
2. 模型怎么求解？（解析解还是数值解）
3. 结果的图形表示：采用哪些图形来表现结果？
1.4.1 运行情况1（a=0.4,b=0.10，delta_t = 0.3）
观察战斗持续了多长时间？
给定一组数据后图形（随时间t 的变化）：
给定一组数据后图形（X，Y军队数量对比图）：
运行结果：
结束时间：t= 10.80,x= -110,y= 1228 1.4.2运行情况2，a=0.15; b=0.1;delta_t=0.05;
程序结果输出：
结束时间：t= 5.85,x= 7882,y= -12 1.4.3运行情况3，a=0.15; b= 0.1;delta_t=0.001;
程序结果输出：
结束时间：t= 5.82,x= 7905,y= -0
思考：
与解析解结果对比，看有多大差别？
1.4.4求解程序
function sim_zhandou
%战斗模型（数学建模问题）
%
%Programmed by: Yong Zhang
%2005.04.06
%pp.131-135
%init
%a=0.4; b=0.1;delta_t = 0.3;
a=0.15; b=0.1; delta_t=0.001;
x0 = 10000;%初始化双方战斗人数（开战前）
y0 = 5000;%常数
x = x0;%x,y存储动态变换人数
y = y0;
k = 0;
%模拟时间轴上离散取值
%（假设战斗最长持续100小时，可以调整）
vec_t = delta_t:delta_t:100;
for cur_t=vec_t,
k = k + 1;%计算第k-1个离散时间点上的人数
%简单的差分法求解微分方程的数值解
x=x-a*y*delta_t;
y=y-b*x*delta_t;
vx(k) = x;%Keep these values
vy(k) = y;
if vx(end)<0 | vy(end)<0,%中止条件（只要一方士兵人数为0或小于0）
disp(sprintf('结束时间：
t=%10.2f,x=%6.0f,y=%6.0f',cur_t,x,y))
break
end
end
%vec_t,vx,vy
t = [0,delta_t:delta_t:cur_t];
close all
%作出军队数量变换与时间t的关系（分别画出曲线）
plot(t,[x0,vx],'r-*',t,[y0,vy],'k-o')
title(sprintf('The curves of two functions (*: Troop X;O: Troop Y,a=%5.2f,b=%5.2f)',a,b))
xlabel('t (unit: Hour)')
hold on
text(t(end)+delta_t*2,vx(end),'X')
text(t(end)+delta_t*2,vy(end),'Y')
hold off
%作战预测模型轨线
figure
plot(vx,vy,'-*')
xlabel('\it{x} (Troop X)')
ylabel('\it{y} (Troop Y)')。