945-微分方程式で解く硫黄島の戦い.sagews

Hiroshi TAKEMOTO([email protected])

微分方程式で解く硫黄島の戦い

大学時代数学セミナーに微分方程式で硫黄島の戦いがシミュレーションできるという記事がありました。

そこで、Googleで「硫黄島の戦い」と「微分方程式」で検索すると、ランチェスターのモデルを使って 硫黄島の戦いを表した本：「微分方程式 その数学と応用」あることが分かりました。

また、Lanchester iwo jimaをキーワードとするとVerifying Lanchester's Combat Model Battle of Iwo Jimaが見つかり、クイック・ビューでその論文を見ることができました。 この論文では、

• xを米軍の兵士の数
• yを日本軍の兵士の数

とし、硫黄島の戦いを以下のような微分方程式で表すことができます。 $$\begin{array}{l l l} \frac{dx}{dt} &=& - a y + f(t) \\ \frac{dy}{dt} &=& - b x \\ \end{array}$$ a, bは米軍、日本軍の戦闘効率係数とし、$f(t)$は米国軍の補強関数で、以下のようになっています。 $$f(t) = \left\{ \begin{array}{l l} 54,000 & 0 \le t \lt 1 \\ 0 & 1 \le t \lt 2 \\ 6,000 & 2 \le t \lt 3 \\ 0 & 3 \le t \lt 5 \\ 13,000 & 5 \le t \lt 6 \\ 0 & t \ge 6 \\ \end{array} \right.$$

微分方程式をsageで表現する

$f(t)$を任意の関数にすることは難しいので、$\delta(t)$として、一瞬で補強が完了するような式を sageを使って定義したのが以下の式です。

# 微分方程式で解く硫黄島の戦い（ランチェスターのモデル）
var('t a b c')
assume(c>0)
x = function('x', t)
y = function('y', t)
f = function('f', t)
# ランチェスターの式を定義します
de1 = diff(x,t) == -a*y + c*dirac_delta(t)
de2 = diff(y,t) == -b*x

微分方程式を解く

desolve_in_japaneseと同様に、以下の手順で微分方程式を解きます。

1. 微分方程式にラプラス変換を施し、微分方程式を変数t領域から変数s領域へ移します
2. 得られた方程式は演算子$s$の代数方程式となるので、代数計算により解を求める
3. 求めたs 関数の解を逆ラプラス変換を用いてt の関数に変換する

var('s')
# 微分方程式をラプラス変換します
lde1 = laplace(de1, t, s); show(lde1)
lde2 = laplace(de2, t, s); show(lde2)

s \mathcal{L}\left(x\left(t\right), t, s\right) - x\left(0\right) = -a \mathcal{L}\left(y\left(t\right), t, s\right) + c

s \mathcal{L}\left(y\left(t\right), t, s\right) - y\left(0\right) = -b \mathcal{L}\left(x\left(t\right), t, s\right)

# x, yのラブラス変換をX, Yに、x(0), y(0)をそれぞれx0, y0に置き換えます
var('X Y x0 y0')
sde1 = lde1.subs_expr(laplace(x(t), t, s) == X, laplace(y(t), t, s) == Y, x(0) == x0, y(0) == y0); show(sde1)
sde2 = lde2.subs_expr(laplace(x(t), t, s) == X, laplace(y(t), t, s) == Y, x(0) == x0, y(0) == y0); show(sde2)
# 以下の式が求める代数方程式になります

X s - x_{0} = -Y a + c

Y s - y_{0} = -X b

# solveを使ってX,Yの解を求める
sol = solve([sde1, sde2], X, Y); sol

[[X == (a*y0 - c*s - s*x0)/(a*b - s^2), Y == (b*c + b*x0 - s*y0)/(a*b - s^2)]]

assume(a > 0)
assume(b > 0)
# solveの解からX,Yの右辺を取り出します
r1 = sol[0][0].rhs(); show(r1)
r2 = sol[0][1].rhs(); show(r2)
# 取り出したX, Yを逆ラプラス変換し、時間領域の解を求めます
s1 = inverse_laplace(r1, s, t); show(s1)
s2 = inverse_laplace(r2, s, t); show(s2)

\frac{{(a y_{0} - c s - s x_{0})}}{{(a b - s^{2})}}

\frac{{(b c + b x_{0} - s y_{0})}}{{(a b - s^{2})}}

{(c + x_{0})} \cosh\left(\sqrt{a} \sqrt{b} t\right) - \frac{\sqrt{a} y_{0} \sinh\left(\sqrt{a} \sqrt{b} t\right)}{\sqrt{b}}

y_{0} \cosh\left(\sqrt{a} \sqrt{b} t\right) - \frac{{(b c + b x_{0})} \sinh\left(\sqrt{a} \sqrt{b} t\right)}{\sqrt{a} \sqrt{b}}

係数の決定

論文では、a, bを以下のように求めています。

最初は、bを求めます。 $$\begin{array}{l c l} \frac{dy}{dt} & = & -b x \\ \int^{y(t_f)}_{y(t_i)} dy & = & -b \int^{t_f}_{t_i} x(t)dt \\ \end{array}$$ から $$\begin{array}{l c l} b & = & \frac{y(t_i) - y(t_f)}{\int^{t_f}_{t_i} x(t)dt} & \approx & \frac{21,500 - 0}{2,024,829} \end{array}$$ となります。

戦闘開始時の日本軍の数は21,500人でほぼ全滅したことから、分子は、21,500とし、 米軍には、戦闘開始から終了までの兵士の数の記録が残っており、分母の$\int^{t_f}_{t_i} x(t)dt = \sum^{36}_{i=1}x(i) = 2,024,829$であることが分かっています。

次にaですが、開戦から28日目で戦況が決まったことから、bの時と同様に $$\begin{array}{r c l} \frac{dx}{dt} & = & f(t) - b x \\ \int^{x(t_f)}_{x(t_i)} dx & = & \int^{t_f}_{t_i} (f(t) - a y(t)) dt \\ x(t_f) - x(t_i) & = & \int^{t_f}_{t_i} f(t)dt - a \int^{t_f}_{t_i} y(t) dt \\ a & = & \frac{\int^{t_f}_{t_i} f(t)dt + x(t_i) - x(t_f)}{\int^{t_f}_{t_i} y(t) dt} \end{array}$$ ここで、$\int^{28}_0 ft(t)dt = 54, 000 + 6, 000 + 13, 000 = 73, 000$であり、日本軍の日々の残留兵士の数の記録はありませんが、先のbを使って $$y(i) \approx 21,500 -b \sum^i_{j=1}{x(j)}$$ から求めることができます。 $$\begin{array}{r c l} a & = & \frac{\int^{28}_{0} f(t)dt + x(0)- x(28)}{\int^{28}_{0} y(t) dt} \\ & \approx & \frac{73,000 + 0 - 52,735}{\sum^{28}_{i=1}{(21,500 - 0.0106*\sum^i_{j=1} x(j))} } \\ & \approx & \frac{20265}{351,172} \\ & \approx & 0.0577 \\ \end{array}$$ となります。

結果の出力

微分方程式の解をプロットしてみます。

米軍の兵士の数をf1(t, a, b, c, x0, y0)の関数にセットし、プロットします。（このようにセットすることで設定を変えた結果をプロットすることができます）

# 論文からa=0.0577, b=0.0106, 米軍の増員が54000人, 日本軍の滞在人数が21500人から米軍の滞在人数をプロットします
f1(t, a, b, c, x0, y0) = s1
plot(f1(t, 0.0577, 0.0106, 54000, 0, 21500), [t, 0, 36])

米軍の兵士の数を実際の記録と重ね合わせてプロットしてみます。f1のtをずらし、x0, y0を再計算しながら計算します。

驚くほどよい一致が見られます。

# 米軍の増員が、初日54000, 2日目6000, 5日目13000を含めてプロット
calcRange = [[0, 2, 54000], [2, 5, 6000], [5, 36, 13000]]
list = [0, 52839 , 50945 , 56031, 56031 , 53749 , 66155 , 65250 , 64378 , 62874 , 62339 , 61405 , 60667 , 59549 , 59345 , 59081 , 58779 , 58196 , 57259 , 56641 , 54792 , 55308 , 54796 , 54038 , 53938 , 53347 , 53072 , 52804 , 52735]
g = [list_plot(list)]
x0 = 21500
y0 = 0
for (ts, te, dx) in calcRange:
    p = plot(f1(t-ts, 0.0577, 0.0106, dx, y0, x0), [ts, te]);
    g.append(p)
    tmp = f2(te-ts, 0.0577, 0.0106, dx, y0, x0)
    y0 = f1(te-ts, 0.0577, 0.0106, dx, y0, x0)
    x0 = tmp
sum(g).show(ymin=0, ymax=70000)

参考までに、日本軍の兵士の推移も表示してみます。

# 同様に日本軍の滞在人数をプロットします
f2(t, a, b, c, x0, y0) = s2
g = []
x0 = 0
y0 = 21500
for (ts, te, dx) in calcRange:
    g.append(plot(f2(t-ts, 0.0544, 0.0106, dx, x0, y0), [ts, te]))
    tmp = f2(te-ts, 0.0577, 0.0106, dx, x0, y0)
    x0 = f1(te-ts, 0.0577, 0.0106, dx, x0, y0)
    y0 = tmp
sum(g).show(ymin=0, ymax=25000)

# 武器の性能（兵力が同じ場合）
x,y = var('x y')
a=1
b=1
plot_vector_field((-a*y, -b*x), [x, 0, 1], [y, 0, 1])

# 硫黄島の戦いの場合
a=0.0577
b=0.0106
plot_vector_field((-a*y, -b*x), [x, 0, 1], [y, 0, 1])