数模算法与应用：预测模型（1）美日硫磺岛战役模型

前言

一、微分方程模型

二、美日硫磺岛战役模型

2.1、题目

2.2、过程解析

2.3、MATLAB建模程序

2.4、结果

前言

预测学是一门研究预测理论、方法及应用的新兴科学。综观预测的思维方式,其基本理论主要有惯性原理、类推原理和相关原理。预测的核心问题是预测的技术方法,或者说是预测的数学模型。随着经济预测、电力预测、资源预测等各种预测的兴起,预测对各种领域的重要性开始显现,预测模型也随着迅速发展。

预测的方法种类繁多,从经典的单耗法、弹性系数法、统计分析法,到目前的灰色预测法、专家系统法和模糊数学法,甚至刚刚兴起的神经网络法、优选组合法和小波分析法,据有关资料统计,预测方法多达 200 余种。因此在使用这些方法建立预测模型时,往往难以正确地判断该用哪种方法,从而不能准确地建立模型,达到要求的效果。不过预测的方法虽然很多,但各种方法多有各自的研究特点、优缺点和适用范围。

一、微分方程模型

当我们描述实际对象的某些特性随时间（或空间）而演变的过程、分析它的变化规律、预测它的未来性态、研究它的控制手段时，通常要建立对象的动态微分方程模型。微分方程大多是物理或几何方面的典型问题，假设条件已经给出，只需用数学符号将已知规律表示出来，即可列出方程，求解的结果就是问题的答案，答案是唯一的，但是有些问题是非物理领域的实际问题，要分析具体情况或进行类比才能给出假设条件。

做出不同的假设，就得到不同的方程。比较典型的有传染病的预测模型、经济增长预测模型、兰彻斯特（ Lanchester ）战争预测模型、药物在体内的分布与排除预测模型、人口的预测模型、烟雾的扩散与消失预测模型等。其基本规律随着时间的增长趋势呈指数形式，根据变量的个数建立微分方程模型。微分方程模型的建立基于相关原理的因果预测法。

该方法的优点是短、中、长期的预测都适合，既能反映内部规律，反映事物的内在关系，也能分析两个因素的相关关系，精度相应的比较高，另外对模型的改进也比较容易理解和实现。该方法的缺点是虽然反映的是内部规律，但是由于方程的建立是以局部规律的独立性假定为基础，故做中长期预测时，偏差有点大。

二、美日硫磺岛战役模型

J . H . Engel 用二次大战末期美日硫黄岛战役中的美军战地记录，对兰彻斯特作战模型进行了验证，发现模型结果与实际数据吻合得很好。

2.1、题目

硫黄岛位于东京以南660英里的海面上，是日军的重要空军基地。美军在1945年2月开始进攻，激烈的战斗持续了一个月，双方伤亡惨重，日方守军21500人全部阵亡或被俘，美方投入兵力73000人，伤亡20265人，战争进行到28天时美军宣布占领该岛，实际战斗到36天才停止。美军的战地记录有按天统计的战斗减员和增援情况。日军没有后援，战地记录则全部遗失。

用 A (t）和 J（t）表示美军和日军第 t 天的人数，忽略双方的非战斗减员，则

(15.1) $\left\{\begin{matrix}\frac{\mathrm{dA(t)} }{\mathrm{d} t}=-aJ(t)+u(t) & & \\\frac{\mathrm{dJ(t)} }{\mathrm{d} t}=-bA(t) & & \\ A(0)=0,J(0)=21500. & & \end{matrix}\right.$

美军战地记录给出增援 u (t）为

$u(t)=\left\{\begin{matrix}54000,0<=t<1; & & & \\ 6000,2<=t<3; & & & \\ 13000,5<=t<6; & & & \\ 0,t>=6 & & & \end{matrix}\right.$

并可由每天伤亡人数算出 A (t), t =1,2,…,36。下面要利用这些实际数据代入式（15.1)，算出 A (t）的理论值，并与实际值比较。

2.2、过程解析

利用给出的数据，对参数 a , b 进行估计。对式（15,1）两边积分，并用求和来近似代替积分，有

为估计 b，在式（15.3）中取 t =36，因为 J (36)=0，且由 A ( t ）的实际数据;

得 $\sum_{t=1}^{36}A(t)=2037000$ ，于是从式(15.3）估计出 b =0.0106。再把这个值代入式（15.3)即可算出 J ( t ), t =1,2,…,36。

然后从式（15.2）估计 a 。令 t =36，得

式中：分子是美军的总伤亡人数，为20265人，分母可由已经算出的 J(t）得到，为372500人。

由式（15.4）有 a =0.0544。把这个值代入式（15.2),得

由式（15.5）就能够算出美军人数 A (t）的理论值，与实际数据吻合得很好。可以根据式（15.3）估计日军的人数。当然也可以求微分方程组（15.1）的数值解，估计日军的人数。

2.3、MATLAB建模程序

下面画出美军人数、日军人数的按时间变化曲线和微分方程组的轨线。

求微分方程组数值解以及画图的MATLAB程序如下：

dxy=@(t,x)[-0.0544*x(2)+54000*(t>=0 & t<1)+...6000*(t>=2 & t<3)+13000*(t>=5 & t<6)-0.0106*x(1)];  %定义微分方程组右端项
[t,xy]=ode45(dxy,[0:36],[0,21500])
subplot(211), plot(t,xy(:,1),'r*',t,xy(:,2),'gD')
xlabel('时间t'),  ylabel('人数'), legend('美军','日军')
subplot(212),  plot(xy(:,1),xy(:,2))  %画微分方程组的轨线
xlabel('美军人数x'),  ylabel('日军人数y')

程序讲解：

dxy=@(t,x)[-0.0544*x(2)+54000*(t>=0 & t<1)+...6000*(t>=2 & t<3)+13000*(t>=5 & t<6)-0.0106*x(1)];  %定义微分方程组右端项

表示 dxy 为函数句柄，@是定义句柄的运算符；@(t,x)是匿名函数：t和x都是参数。

[t,xy]=ode45(dxy,[0:36],[0,21500])

ode是Matlab专门用于解微分方程的功能函数。

ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差。是一种自适应步长（变步长）的常微分方程数值解法，其整体截断误差为(Δx)^5。

解决的是Nonstiff(非刚性)常微分方程。

ode45是解决数值解问题的首选方法，若长时间没结果，应该就是刚性的，可换用ode15s试试。

2.4、结果

t =0       1       2       3       4       5       6       7       8       9       10       11       12       13       14       15       16       17       18       19       20       21       22       23       24       25       26       27       28       29       30       31       32       33       34       35       36       xy =0          21500       50540         233425/11    245486/5       269008/13    481157/9       141063/7     637775/12      254562/13    260306/5        19025       310297/5       184179/10    182863/3       177659/10    300027/5       102749/6     59091          49481/3     232843/4        15872       57364         289930/19    169651/3        58623/4     167308/3        28121/2     275102/5        40420/3     217213/4       244985/19    268087/5        12322       105925/2        82300/7     471043/9       134389/12    620927/12      202302/19    460616/9        10102       151934/3       105186/11    701947/14      153480/17    148987/3       135989/16    196857/4        55827/7     97589/2        59647/8     48403          27763/4     144118/3        51437/8     143110/3        82911/14    616134/13       48764/9     235569/5       127851/26    328019/7        57451/13    792760/17       66705/17    510762/11       17153/5     138779/3        58787/20    599470/13       12249/5     505924/11       80428/41