基于MATLAB的差分方程战斗模型

MATLAB实现战争模型的数值解法

摘要

战争模型是一种非常典型的连续系统仿真模型，战争系统中变量的关系具有很高的学习研究的价值。本文将通过MATLAB软件来实现简单的战斗模型模拟，通过差分方程的形式实现数值解的求解。除此之外，本文还将探讨在有增援条件下的战争模型的变化，并将进一步分析增援总人数和单位时间增援人数的关系

模型简介

假定 X 部队t时刻存活的士兵数为x(t)，而 Y 部队在t时刻存活的士兵数为 y(t)，将
x(t)与 y(t)都看作连续变量。并假定双方所有士兵不是战死就是活着参加战斗，即不考虑
俘虏和伤员。
关于双方的作战伤亡情况，一种合理的假设是，在Δt 时间内，X 部队被杀死的士兵数
Δx 将取决于Δt 的长短，以及在Δt起始时刻与其交战的 Y部队的士兵数。假定是一种正比
关系，即
其中a 是一个常数，代表了Y 部队的战斗力，称为“杀伤率”，更明确地说，a 是Y 部队的
一个士兵在单位时间内杀死X 部队的名士兵数。类似地，对于Y 部队有
根据微分方程的形式，我们可以得到如下的差分方程

代码实现

现在根据之前得到的差分方程来进行离散化的数值求解：

变量	含义	变量值
t	离散的时刻	为一个矩阵，初值为0
dt	时间步长	0.001
a	Y军队的战斗力	0.15
b	X军队的战斗力	0.1
x	X军队的剩余人数	矩阵，初始值为10000
y	Y军队的剩余人数	矩阵，初始值为5000
totallr	已经增援的数量	初始值为0
p	每隔固定的时间增援的人数	变量
q	增援的最大人数	变量
time_interval	增援时间间隔	设为100倍dt，0.1

当Y方没有增援时，p=0

clear all
t=0;
a=0.15;b=0.1;
dt = 0.001 ;
x(1)=10000;y(1)=5000;t(1)=0;
r=2.5;%假定y获得增援的时刻
totalr=0;
p=0;q=6000;
time_interval=0.1
for i=1: 1e5
x(i+1)=x(i)-a*y(i)*dt; %用导数定义理解近似解法：函数增量/自变量增量apprx 导数
y(i+1)=y(i)-b*x(i)*dt;
t(i+1)=t(i)+dt;
   if t(i+1)>r&&rem(t(i+1),time_interval)<0.000001&&totalr<q
    y(i+1)=y(i+1)+p;
    totalr=totalr+p;   
   end
if x(i+1)<1 ||y(i+1)<1,%战斗结束条件:假设直到一方剩余士兵数量<1
T=t(i+1);
break
end
end
T%战争结束的时刻
lastx=x(i+1)%结束战争时X的剩余人数
lasty=y(i+1)%结束战争时Y的剩余人数
t_=1:length(x);
plot(t_,x,'r',t_,y,'g')

Y会输。当p=700,q=6000时：
Y赢了。
现在来通过全局搜索法探究,当Y要赢的时候p 和q要满足的大致关系

clear all
t=0;
a=0.15;b=0.1;
dt = 0.001 ;
x(1)=10000;y(1)=5000;t(1)=0;
r=2;%假定y获得增援的时刻
totalr=0;
avalable=[];
num=0;
time_interval=0.1;
for p=0:1:800;
for q=0:25:15000
x(1)=10000;y(1)=5000;t(1)=0;
totalr=0;
check=0;
for i=1: 1e5
x(i+1)=x(i)-a*y(i)*dt; %用导数定义理解近似解法：函数增量/自变量增量apprx 导数
y(i+1)=y(i)-b*x(i)*dt;
t(i+1)=t(i)+dt;
   if t(i+1)>r&&rem(t(i+1),time_interval)<0.000001&&totalr<q%注意这个地方rem不能让他判断是否等于0，因为小数除法存在计算机误差，可以整除时余数不是绝对的0，而是0.0..(若干个0)..0xxx.
    y(i+1)=y(i+1)+p;
    totalr=totalr+p;   
   end
if x(i+1)<1%X输了的时候直接结束这次循环
    check=1; 
    break
end

if x(i+1)<1 ||y(i+1)<1,%战斗结束条件:假设直到一方剩余士兵数量<1
T=t(i+1);
break
end

end
if check==1
    avalable=[avalable;p,q];
    num=num+1
    break
end  
end
end
T
lastx=x(i+1)
lasty=y(i+1)
t_=1:length(x);
%plot(t_,x,'r',t_,y,'g')
plot(avalable(:,1),avalable(:,2))
xlabel('每次增援')
ylabel('增援总人数')

结果如上图所示。这里出现锯齿状的线条推测是因为差分方程的精度不够，导致可能的周期性波动。经查阅资料得知，理论上的曲线应为一大致趋势相同但平滑的曲线，所以这次实验有待改进。