Lanchester_model_in_japanese -- Sage

Hiroshi TAKEMOTO(take@pwv.co.jp)

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

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

そこで、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と同様に、以下の手順で微分方程式を解きます。

微分方程式にラプラス変換を施し、微分方程式を変数t領域から変数s領域へ移します
得られた方程式は演算子$s$の代数方程式となるので、代数計算により解を求める
求めた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)

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

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)]]

[[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}}

\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} = 0.0106 \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])

増兵関数を任意の関数にした場合

これまでは、増兵関数$f(t)$を$\delta(t)$を使いましたが、任意の関数にした場合を見てみましょう。

$\delta(t)$から$f(t)$に変更し、$f(t)$のラプラス変換を$F(s)$として方程式を解きます。

# 微分方程式で解く硫黄島の戦い（ランチェスターのモデル） var('t a b') x = function('x', t) y = function('y', t) f = function('f', t) # ランチェスターの式を定義します de1 = diff(x,t) == -a*y + f de2 = diff(y,t) == -b*x var('s') # 微分方程式をラプラス変換します lde1 = laplace(de1, t, s); lde2 = laplace(de2, t, s); # x, y, fのラブラス変換をX, Y, Fに、x(0), y(0)をそれぞれx0, y0に置き換えます var('X Y x0 y0') F = function('F', s) sde1 = lde1.subs_expr(laplace(x(t), t, s) == X, laplace(f(t), t, s) == F, 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(f(t), t, s) == F, laplace(y(t), t, s) == Y, x(0) == x0, y(0) == y0); show(sde2) # 以下の式が求める代数方程式になります

X s - x_{0} = -Y a + F\left(s\right)
Y s - y_{0} = -X b

X s - x_{0} = -Y a + F\left(s\right)
Y s - y_{0} = -X b

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

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

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

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

\mathcal{L}^{-1}\left(\frac{{(a y_{0} - s x_{0} - s F\left(s\right))}}{{(a b - s^{2})}}, s, t\right)
\mathcal{L}^{-1}\left(\frac{{(b x_{0} + b F\left(s\right) - s y_{0})}}{{(a b - s^{2})}}, s, t\right)

\mathcal{L}^{-1}\left(\frac{{(a y_{0} - s x_{0} - s F\left(s\right))}}{{(a b - s^{2})}}, s, t\right)
\mathcal{L}^{-1}\left(\frac{{(b x_{0} + b F\left(s\right) - s y_{0})}}{{(a b - s^{2})}}, s, t\right)

任意の増兵関数の場合の解

$f(t)$が入っただけで解けなくなってしまいましたが、最初の式を分解して見てみると、 $$ x(t) = \frac{a y0}{(a b - s^2)} + \frac{- s x0}{(a b - s^2)} + \frac{- s F(s)}{(a b - s^2)} $$ となっています。 $$ \frac{- s F(s)}{(ab - s^2) } $$ 以外は、$\delta(t)$の時の解に似ていることが分かります。

ラブラス変換の積

ラブラス変換の積で検索すると、以下のページが見つかりました。ラプラス変換の合成法則

これによると、ラプラス変換関数の積は、 $$ \begin{array}{r c l} F(s)G(s) & = & \mathcal{L}\left( f(t)*g(t) \right) & = & \mathcal{L}\left( \int_{0}^{t} f(t - \tau)g(\tau)d\tau \right) \end{array} $$ ここで、$f(t)*g(t)$の$*$は、かけ算ではなく、「たたみ込み計算」の記号です。

最終的な解

$sinh(\omega t), cosh(\omega t)$のラプラス変換の式から、$\omega = \sqrt{ab}$となり、 $$ \begin{array}{r c l} x(t) & = & -\frac{\sqrt{a} y_0 sinh(\sqrt{ab}t)} {\sqrt{b}} + x_0 cosh(\sqrt{ab}t) + \int_{0}^{t} cosh(\sqrt{ab}(t - \tau))f(\tau)d\tau \\ y(t) & = & y_0 cosh(\sqrt{ab}t) - \frac{bx_0 sinh(\sqrt{ab}t)}{\sqrt{ab}} \\ \end{array} $$

# sinh(ωt), cosh(ωt)のラプラス変換 var('omega') laplace(sinh(omega*t), t, s).show() laplace(cosh(omega*t), t, s).show()

-\frac{\omega}{{(\omega^{2} - s^{2})}}
-\frac{s}{{(\omega^{2} - s^{2})}}

-\frac{\omega}{{(\omega^{2} - s^{2})}}
-\frac{s}{{(\omega^{2} - s^{2})}}

Lanchester_model_in_japanese

5031 days ago by takepwave