sage/硫黄島の戦い
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2009/11/09からのアクセス回数 &counter;
このページのsage notebookを以下のURLで公開しています。
http://www.sagenb.org/pub/945/
注意：このページは、「ランチェスターのモデル」がどの程度...
** 微分方程式で解く硫黄島の戦い [#z3d265c0]
大学時代数学セミナーに微分方程式で硫黄島の戦いがシミュレ...
そこで、Googleで「硫黄島の戦い」と「微分方程式」で検索す...
また、Lanchester iwo jimaをキーワードとすると[[Verifying ...
- 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で表現する [#z3cc1664]
\(f(t)\)を任意の関数にすることは難しいので、\(\delta(t)\)...
sage入力：
#pre{{
# 微分方程式で解く硫黄島の戦い（ランチェスターのモデル）
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
}}
*** 微分方程式を解く [#bb3c25f5]
[[desolve_in_japanese>http://www.sagenb.org/home/pub/906/...
+ 微分方程式にラプラス変換を施し、微分方程式を変数t領域か...
+ 得られた方程式は演算子sの代数方程式となるので、代数計算...
+ 求めたs 関数の解を逆ラプラス変換を用いてt の関数に変換...
sage入力：
#pre{{
var('s')
# 微分方程式をラプラス変換します
lde1 = laplace(de1, t, s); show(lde1)
lde2 = laplace(de2, t, s); show(lde2)
}}
&color(blue){\(s \mathcal{L}\left(x\left(t\right), t, s\r...
&color(blue){\(s \mathcal{L}\left(y\left(t\right), t, s\r...
sage入力：
#pre{{
# 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...
sde2 = lde2.subs_expr(laplace(x(t), t, s) == X, laplace(y...
# 以下の式が求める代数方程式になります
}}
&color(blue){\(X s - x_{0} = -Y a + c\)};
&color(blue){\(Y s - y_{0} = -X b\)};
sage入力：
#pre{{
# solveを使ってX,Yの解を求める
sol = solve([sde1, sde2], X, Y); sol
}}
&color(blue){[ [X == (a*y0 - c*s - s*x0)/(a*b - s^2), Y =...
sage入力：
#pre{{
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)
}}
&color(blue){\( \frac{{(a y_{0} - c s - s x_{0})}}{{(a b ...
&color(blue){\( \frac{{(b c + b x_{0} - s y_{0})}}{{(a b ...
&color(blue){\( {(c + x_{0})} \cosh\left(\sqrt{a} \sqrt{b...
&color(blue){\( y_{0} \cosh\left(\sqrt{a} \sqrt{b} t\righ...
** 係数の決定 [#z99d610a]
論文では、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)d...
\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人でほぼ全滅したことから、...
米軍には、戦闘開始から終了までの兵士の数の記録が残ってお...
次に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) - ...
x(t_f) - x(t_i) & = & \int^{t_f}_{t_i} f(t)dt - a \int^{t...
a & = & \frac{\int^{t_f}_{t_i} f(t)dt + x(t_i) - x(t_f)}{...
\end{array}
$$
ここで、\(\int^{28}_0 ft(t)dt = 54, 000 + 6, 000 + 13, 00...
$$
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^{2...
& \approx & \frac{73,000 + 0 - 52,735}{\sum^{28}_{i=1}{(2...
& \approx & \frac{20265}{351,172} \\
& \approx & 0.0577 \\
\end{array}
$$
となります。
** 結果の出力 [#o151938d]
微分方程式の解をプロットしてみます。
米軍の兵士の数をf1(t, a, b, c, x0, y0)の関数にセットし、...
sage入力：
#pre{{
# 論文からa=0.0577, b=0.0106, 米軍の増員が54000人, 日本軍...
f1(t, a, b, c, x0, y0) = s1
plot(f1(t, 0.0577, 0.0106, 54000, 0, 21500), [t, 0, 36])
}}
&ref(sage0.png);
米軍の兵士の数を実際の記録と重ね合わせてプロットしてみま...
驚くほどよい一致が見られます。
sage入力：
#pre{{
# 米軍の増員が、初日54000, 2日目6000, 5日目13000を含めて...
calcRange = [[0, 2, 54000], [2, 5, 6000], [5, 36, 13000]]
list = [0, 52839 , 50945 , 56031, 56031 , 53749 , 66155 ,...
62339 , 61405 , 60667 , 59549 , 59345 , 59081 , 58779 , ...
54792 , 55308 , 54796 , 54038 , 53938 , 53347 , 53072 , ...
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, t...
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)
}}
&ref(sage0-1.png);
参考までに、日本軍の兵士の推移も表示してみます。
#pre{{
# 同様に日本軍の滞在人数をプロットします
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), [...
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)
}}
&ref(sage0-2.png);
*** 武器の性能の比較 [#facb9a1f]
武器の性能が同じ場合には、両者のベクトル図は、以下のよう...
sage入力:
#pre{{
# 武器の性能（兵力が同じ場合）
x,y = var('x y')
a=1
b=1
plot_vector_field((-a*y, -b*x), [x, 0, 1], [y, 0, 1])
}}
&ref(sage0-3.png);
計算で求まったa, bを使うと以下のようになります。
sage入力:
#pre{{
# 硫黄島の戦いの場合
a=0.0577
b=0.0106
plot_vector_field((-a*y, -b*x), [x, 0, 1], [y, 0, 1])
}}
&ref(sage0-4.png);
** コメント [#g98a33c0]
#vote(おもしろかった[12],そうでもない[0],わかりずらい[2])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2009/11/09からのアクセス回数 &counter;
このページのsage notebookを以下のURLで公開しています。
http://www.sagenb.org/pub/945/
注意：このページは、「ランチェスターのモデル」がどの程度...
** 微分方程式で解く硫黄島の戦い [#z3d265c0]
大学時代数学セミナーに微分方程式で硫黄島の戦いがシミュレ...
そこで、Googleで「硫黄島の戦い」と「微分方程式」で検索す...
また、Lanchester iwo jimaをキーワードとすると[[Verifying ...
- 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で表現する [#z3cc1664]
\(f(t)\)を任意の関数にすることは難しいので、\(\delta(t)\)...
sage入力：
#pre{{
# 微分方程式で解く硫黄島の戦い（ランチェスターのモデル）
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
}}
*** 微分方程式を解く [#bb3c25f5]
[[desolve_in_japanese>http://www.sagenb.org/home/pub/906/...
+ 微分方程式にラプラス変換を施し、微分方程式を変数t領域か...
+ 得られた方程式は演算子sの代数方程式となるので、代数計算...
+ 求めたs 関数の解を逆ラプラス変換を用いてt の関数に変換...
sage入力：
#pre{{
var('s')
# 微分方程式をラプラス変換します
lde1 = laplace(de1, t, s); show(lde1)
lde2 = laplace(de2, t, s); show(lde2)
}}
&color(blue){\(s \mathcal{L}\left(x\left(t\right), t, s\r...
&color(blue){\(s \mathcal{L}\left(y\left(t\right), t, s\r...
sage入力：
#pre{{
# 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...
sde2 = lde2.subs_expr(laplace(x(t), t, s) == X, laplace(y...
# 以下の式が求める代数方程式になります
}}
&color(blue){\(X s - x_{0} = -Y a + c\)};
&color(blue){\(Y s - y_{0} = -X b\)};
sage入力：
#pre{{
# solveを使ってX,Yの解を求める
sol = solve([sde1, sde2], X, Y); sol
}}
&color(blue){[ [X == (a*y0 - c*s - s*x0)/(a*b - s^2), Y =...
sage入力：
#pre{{
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)
}}
&color(blue){\( \frac{{(a y_{0} - c s - s x_{0})}}{{(a b ...
&color(blue){\( \frac{{(b c + b x_{0} - s y_{0})}}{{(a b ...
&color(blue){\( {(c + x_{0})} \cosh\left(\sqrt{a} \sqrt{b...
&color(blue){\( y_{0} \cosh\left(\sqrt{a} \sqrt{b} t\righ...
** 係数の決定 [#z99d610a]
論文では、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)d...
\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人でほぼ全滅したことから、...
米軍には、戦闘開始から終了までの兵士の数の記録が残ってお...
次に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) - ...
x(t_f) - x(t_i) & = & \int^{t_f}_{t_i} f(t)dt - a \int^{t...
a & = & \frac{\int^{t_f}_{t_i} f(t)dt + x(t_i) - x(t_f)}{...
\end{array}
$$
ここで、\(\int^{28}_0 ft(t)dt = 54, 000 + 6, 000 + 13, 00...
$$
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^{2...
& \approx & \frac{73,000 + 0 - 52,735}{\sum^{28}_{i=1}{(2...
& \approx & \frac{20265}{351,172} \\
& \approx & 0.0577 \\
\end{array}
$$
となります。
** 結果の出力 [#o151938d]
微分方程式の解をプロットしてみます。
米軍の兵士の数をf1(t, a, b, c, x0, y0)の関数にセットし、...
sage入力：
#pre{{
# 論文からa=0.0577, b=0.0106, 米軍の増員が54000人, 日本軍...
f1(t, a, b, c, x0, y0) = s1
plot(f1(t, 0.0577, 0.0106, 54000, 0, 21500), [t, 0, 36])
}}
&ref(sage0.png);
米軍の兵士の数を実際の記録と重ね合わせてプロットしてみま...
驚くほどよい一致が見られます。
sage入力：
#pre{{
# 米軍の増員が、初日54000, 2日目6000, 5日目13000を含めて...
calcRange = [[0, 2, 54000], [2, 5, 6000], [5, 36, 13000]]
list = [0, 52839 , 50945 , 56031, 56031 , 53749 , 66155 ,...
62339 , 61405 , 60667 , 59549 , 59345 , 59081 , 58779 , ...
54792 , 55308 , 54796 , 54038 , 53938 , 53347 , 53072 , ...
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, t...
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)
}}
&ref(sage0-1.png);
参考までに、日本軍の兵士の推移も表示してみます。
#pre{{
# 同様に日本軍の滞在人数をプロットします
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), [...
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)
}}
&ref(sage0-2.png);
*** 武器の性能の比較 [#facb9a1f]
武器の性能が同じ場合には、両者のベクトル図は、以下のよう...
sage入力:
#pre{{
# 武器の性能（兵力が同じ場合）
x,y = var('x y')
a=1
b=1
plot_vector_field((-a*y, -b*x), [x, 0, 1], [y, 0, 1])
}}
&ref(sage0-3.png);
計算で求まったa, bを使うと以下のようになります。
sage入力:
#pre{{
# 硫黄島の戦いの場合
a=0.0577
b=0.0106
plot_vector_field((-a*y, -b*x), [x, 0, 1], [y, 0, 1])
}}
&ref(sage0-4.png);
** コメント [#g98a33c0]
#vote(おもしろかった[12],そうでもない[0],わかりずらい[2])
皆様のご意見、ご希望をお待ちしております。
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
