︠42001062-bb73-48b0-b201-eb0479a9aa1bi︠ %html Hiroshi TAKEMOTO(take@pwv.co.jp)

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

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

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

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

とし、硫黄島の戦いを以下のような微分方程式で表すことができます。 $$ \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を使って定義したのが以下の式です。 ︡cee2fe4b-c893-448a-bcc3-a3fc7e2bc12d︡{"html": "Hiroshi TAKEMOTO(take@pwv.co.jp)\r\n\r\n

\u5fae\u5206\u65b9\u7a0b\u5f0f\u3067\u89e3\u304f\u786b\u9ec4\u5cf6\u306e\u6226\u3044

\r\n\u5927\u5b66\u6642\u4ee3\u6570\u5b66\u30bb\u30df\u30ca\u30fc\u306b\u5fae\u5206\u65b9\u7a0b\u5f0f\u3067\u786b\u9ec4\u5cf6\u306e\u6226\u3044\u304c\u30b7\u30df\u30e5\u30ec\u30fc\u30b7\u30e7\u30f3\u3067\u304d\u308b\u3068\u3044\u3046\u8a18\u4e8b\u304c\u3042\u308a\u307e\u3057\u305f\u3002

\r\n\r\n\u305d\u3053\u3067\u3001Google\u3067\u300c\u786b\u9ec4\u5cf6\u306e\u6226\u3044\u300d\u3068\u300c\u5fae\u5206\u65b9\u7a0b\u5f0f\u300d\u3067\u691c\u7d22\u3059\u308b\u3068\u3001\u30e9\u30f3\u30c1\u30a7\u30b9\u30bf\u30fc\u306e\u30e2\u30c7\u30eb\u3092\u4f7f\u3063\u3066\r\n\u786b\u9ec4\u5cf6\u306e\u6226\u3044\u3092\u8868\u3057\u305f\u672c\uff1a\u300c\u5fae\u5206\u65b9\u7a0b\u5f0f \u305d\u306e\u6570\u5b66\u3068\u5fdc\u7528\u300d\u3042\u308b\u3053\u3068\u304c\u5206\u304b\u308a\u307e\u3057\u305f\u3002

\r\n\r\n\u307e\u305f\u3001Lanchester iwo jima\u3092\u30ad\u30fc\u30ef\u30fc\u30c9\u3068\u3059\u308b\u3068Verifying Lanchester's Combat Model Battle of Iwo Jima\u304c\u898b\u3064\u304b\u308a\u3001\u30af\u30a4\u30c3\u30af\u30fb\u30d3\u30e5\u30fc\u3067\u305d\u306e\u8ad6\u6587\u3092\u898b\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3057\u305f\u3002\r\n\r\n\u3053\u306e\u8ad6\u6587\u3067\u306f\u3001\r\n

\r\n\u3068\u3057\u3001\u786b\u9ec4\u5cf6\u306e\u6226\u3044\u3092\u4ee5\u4e0b\u306e\u3088\u3046\u306a\u5fae\u5206\u65b9\u7a0b\u5f0f\u3067\u8868\u3059\u3053\u3068\u304c\u3067\u304d\u307e\u3059\u3002\r\n$$\r\n\\begin{array}{l l l}\r\n\\frac{dx}{dt} &=& - a y + f(t) \\\\\r\n\\frac{dy}{dt} &=& - b x \\\\\r\n\\end{array}\r\n$$\r\na, b\u306f\u7c73\u8ecd\u3001\u65e5\u672c\u8ecd\u306e\u6226\u95d8\u52b9\u7387\u4fc2\u6570\u3068\u3057\u3001$f(t)$\u306f\u7c73\u56fd\u8ecd\u306e\u88dc\u5f37\u95a2\u6570\u3067\u3001\u4ee5\u4e0b\u306e\u3088\u3046\u306b\u306a\u3063\u3066\u3044\u307e\u3059\u3002\r\n$$\r\nf(t) = \\left\\{\r\n \\begin{array}{l l}\r\n 54,000 & 0 \\le t \\lt 1 \\\\\r\n 0 & 1 \\le t \\lt 2 \\\\\r\n 6,000 & 2 \\le t \\lt 3 \\\\\r\n 0 & 3 \\le t \\lt 5 \\\\\r\n 13,000 & 5 \\le t \\lt 6 \\\\\r\n 0 & t \\ge 6 \\\\\r\n \\end{array}\r\n\\right.\r\n$$\r\n

\r\n\r\n

\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092sage\u3067\u8868\u73fe\u3059\u308b

\r\n$f(t)$\u3092\u4efb\u610f\u306e\u95a2\u6570\u306b\u3059\u308b\u3053\u3068\u306f\u96e3\u3057\u3044\u306e\u3067\u3001$\\delta(t)$\u3068\u3057\u3066\u3001\u4e00\u77ac\u3067\u88dc\u5f37\u304c\u5b8c\u4e86\u3059\u308b\u3088\u3046\u306a\u5f0f\u3092\r\nsage\u3092\u4f7f\u3063\u3066\u5b9a\u7fa9\u3057\u305f\u306e\u304c\u4ee5\u4e0b\u306e\u5f0f\u3067\u3059\u3002\r\n"}︡ ︠826c5bda-ce4b-42ec-981b-2f2d46792ac8︠ # 微分方程式で解く硫黄島の戦い（ランチェスターのモデル） 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 ︡7fbd443c-9ef8-4cd9-930d-c2f3289f857a︡︡ ︠9f9b9054-d7b5-4164-ba4f-77d34fe45e52i︠ %html

微分方程式を解く

desolve_in_japaneseと同様に、以下の手順で微分方程式を解きます。
  1. 微分方程式にラプラス変換を施し、微分方程式を変数t領域から変数s領域へ移します
  2. 得られた方程式は演算子$s$の代数方程式となるので、代数計算により解を求める
  3. 求めたs 関数の解を逆ラプラス変換を用いてt の関数に変換する
︡1fe2e8f6-656e-4187-9662-3124102b81fb︡{"html": "\r\n

\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092\u89e3\u304f

\r\ndesolve_in_japanese\u3068\u540c\u69d8\u306b\u3001\u4ee5\u4e0b\u306e\u624b\u9806\u3067\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092\u89e3\u304d\u307e\u3059\u3002\r\n
    \r\n
  1. \u5fae\u5206\u65b9\u7a0b\u5f0f\u306b\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3092\u65bd\u3057\u3001\u5fae\u5206\u65b9\u7a0b\u5f0f\u3092\u5909\u6570t\u9818\u57df\u304b\u3089\u5909\u6570s\u9818\u57df\u3078\u79fb\u3057\u307e\u3059
    2. \r\n
  2. \u5f97\u3089\u308c\u305f\u65b9\u7a0b\u5f0f\u306f\u6f14\u7b97\u5b50$s$\u306e\u4ee3\u6570\u65b9\u7a0b\u5f0f\u3068\u306a\u308b\u306e\u3067\u3001\u4ee3\u6570\u8a08\u7b97\u306b\u3088\u308a\u89e3\u3092\u6c42\u3081\u308b
    3. \r\n
  3. \u6c42\u3081\u305fs \u95a2\u6570\u306e\u89e3\u3092\u9006\u30e9\u30d7\u30e9\u30b9\u5909\u63db\u3092\u7528\u3044\u3066t \u306e\u95a2\u6570\u306b\u5909\u63db\u3059\u308b
    4. \r\n
\r\n\r\n"}︡ ︠cd98c597-949d-485f-98f9-a15834ac06be︠ var('s') # 微分方程式をラプラス変換します lde1 = laplace(de1, t, s); show(lde1) lde2 = laplace(de2, t, s); show(lde2) ︡2b94fcf8-2ec1-44aa-b5f6-a3a46e04299b︡{"html": "
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
"}︡{"html": "
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)
"}︡ ︠ac6ecdeb-96d9-41de-b2c6-a9c5c706d784︠ # 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) # 以下の式が求める代数方程式になります ︡1171bb6a-1fe9-41cc-b264-8126c70afac2︡{"html": "
X s - x_{0} = -Y a + c
"}︡{"html": "
Y s - y_{0} = -X b
"}︡ ︠d2a36b61-edca-460a-94ff-b22f511a6094︠ # solveを使ってX,Yの解を求める sol = solve([sde1, sde2], X, Y); sol ︡3fadf12f-4f4b-4f80-b911-a70af5e2b34e︡{"stdout": "[[X == (a*y0 - c*s - s*x0)/(a*b - s^2), Y == (b*c + b*x0 - s*y0)/(a*b - s^2)]]"}︡ ︠3b0e922b-e01e-47ae-ab2a-4d40870b0d14︠ 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) ︡b0bced87-0ec2-4aa4-8e50-640380ed445b︡{"html": "
\\frac{{(a y_{0} - c s - s x_{0})}}{{(a b - s^{2})}}
"}︡{"html": "
\\frac{{(b c + b x_{0} - s y_{0})}}{{(a b - s^{2})}}
"}︡{"html": "
{(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}}
"}︡{"html": "
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}}
"}︡ ︠fca84b50-de8a-4d74-9f65-fcb43fb4b3b3i︠ %html

係数の決定

論文では、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)の関数にセットし、プロットします。（このようにセットすることで設定を変えた結果をプロットすることができます） ︡1d6ece5a-f8ce-4e44-b579-18f81001d092︡{"html": "\r\n

\u4fc2\u6570\u306e\u6c7a\u5b9a

\r\n\u8ad6\u6587\u3067\u306f\u3001a, b\u3092\u4ee5\u4e0b\u306e\u3088\u3046\u306b\u6c42\u3081\u3066\u3044\u307e\u3059\u3002

\r\n\r\n\u6700\u521d\u306f\u3001b\u3092\u6c42\u3081\u307e\u3059\u3002\r\n$$\r\n\\begin{array}{l c l}\r\n\\frac{dy}{dt} & = & -b x \\\\\r\n\\int^{y(t_f)}_{y(t_i)} dy & = & -b \\int^{t_f}_{t_i} x(t)dt \\\\\r\n\\end{array}\r\n$$\r\n\u304b\u3089\r\n$$\r\n\\begin{array}{l c l}\r\nb & = & \\frac{y(t_i) - y(t_f)}{\\int^{t_f}_{t_i} x(t)dt}\r\n & \\approx & \\frac{21,500 - 0}{2,024,829}\r\n\\end{array}\r\n$$\r\n\u3068\u306a\u308a\u307e\u3059\u3002

\r\n\r\n\u6226\u95d8\u958b\u59cb\u6642\u306e\u65e5\u672c\u8ecd\u306e\u6570\u306f21,500\u4eba\u3067\u307b\u307c\u5168\u6ec5\u3057\u305f\u3053\u3068\u304b\u3089\u3001\u5206\u5b50\u306f\u300121,500\u3068\u3057\u3001\r\n\u7c73\u8ecd\u306b\u306f\u3001\u6226\u95d8\u958b\u59cb\u304b\u3089\u7d42\u4e86\u307e\u3067\u306e\u5175\u58eb\u306e\u6570\u306e\u8a18\u9332\u304c\u6b8b\u3063\u3066\u304a\u308a\u3001\u5206\u6bcd\u306e$\\int^{t_f}_{t_i} x(t)dt = \\sum^{36}_{i=1}x(i) = 2,024,829$\u3067\u3042\u308b\u3053\u3068\u304c\u5206\u304b\u3063\u3066\u3044\u307e\u3059\u3002

\r\n\r\n\u6b21\u306ba\u3067\u3059\u304c\u3001\u958b\u6226\u304b\u308928\u65e5\u76ee\u3067\u6226\u6cc1\u304c\u6c7a\u307e\u3063\u305f\u3053\u3068\u304b\u3089\u3001b\u306e\u6642\u3068\u540c\u69d8\u306b\r\n$$\r\n\\begin{array}{r c l}\r\n\\frac{dx}{dt} & = & f(t) - b x \\\\\r\n\\int^{x(t_f)}_{x(t_i)} dx & = & \\int^{t_f}_{t_i} (f(t) - a y(t)) dt \\\\\r\nx(t_f) - x(t_i) & = & \\int^{t_f}_{t_i} f(t)dt - a \\int^{t_f}_{t_i} y(t) dt \\\\\r\na & = & \\frac{\\int^{t_f}_{t_i} f(t)dt + x(t_i) - x(t_f)}{\\int^{t_f}_{t_i} y(t) dt} \r\n\\end{array}\r\n$$\r\n\u3053\u3053\u3067\u3001$\\int^{28}_0 ft(t)dt = 54, 000 + 6, 000 + 13, 000 = 73, 000 $\u3067\u3042\u308a\u3001\u65e5\u672c\u8ecd\u306e\u65e5\u3005\u306e\u6b8b\u7559\u5175\u58eb\u306e\u6570\u306e\u8a18\u9332\u306f\u3042\u308a\u307e\u305b\u3093\u304c\u3001\u5148\u306eb\u3092\u4f7f\u3063\u3066\r\n$$\r\ny(i) \\approx 21,500 -b \\sum^i_{j=1}{x(j)}\r\n$$\r\n\u304b\u3089\u6c42\u3081\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3059\u3002\r\n$$\r\n\\begin{array}{r c l}\r\na & = & \\frac{\\int^{28}_{0} f(t)dt + x(0)- x(28)}{\\int^{28}_{0} y(t) dt} \\\\\r\n & \\approx & \\frac{73,000 + 0 - 52,735}{\\sum^{28}_{i=1}{(21,500 - 0.0106*\\sum^i_{j=1} x(j))} } \\\\\r\n & \\approx & \\frac{20265}{351,172} \\\\\r\n & \\approx & 0.0577 \\\\\r\n\\end{array}\r\n$$\r\n\u3068\u306a\u308a\u307e\u3059\u3002\r\n\r\n\r\n\r\n\r\n

\u7d50\u679c\u306e\u51fa\u529b

\r\n\u5fae\u5206\u65b9\u7a0b\u5f0f\u306e\u89e3\u3092\u30d7\u30ed\u30c3\u30c8\u3057\u3066\u307f\u307e\u3059\u3002

\r\n\u7c73\u8ecd\u306e\u5175\u58eb\u306e\u6570\u3092f1(t, a, b, c, x0, y0)\u306e\u95a2\u6570\u306b\u30bb\u30c3\u30c8\u3057\u3001\u30d7\u30ed\u30c3\u30c8\u3057\u307e\u3059\u3002\uff08\u3053\u306e\u3088\u3046\u306b\u30bb\u30c3\u30c8\u3059\u308b\u3053\u3068\u3067\u8a2d\u5b9a\u3092\u5909\u3048\u305f\u7d50\u679c\u3092\u30d7\u30ed\u30c3\u30c8\u3059\u308b\u3053\u3068\u304c\u3067\u304d\u307e\u3059\uff09\r\n"}︡ ︠693fa3eb-d25f-45b5-b8a8-242ff978779a︠ # 論文から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]) ︡b5af7e5d-0dd7-47c2-8068-b15e6296b450︡{"html": ""}︡ ︠54172d05-d9b2-42d3-b8cd-48c6dfc1de6bi︠ %html 米軍の兵士の数を実際の記録と重ね合わせてプロットしてみます。f1のtをずらし、x0, y0を再計算しながら計算します。

驚くほどよい一致が見られます。 ︡f584c912-e613-4ea8-b087-03f65ef91339︡{"html": "\r\n\u7c73\u8ecd\u306e\u5175\u58eb\u306e\u6570\u3092\u5b9f\u969b\u306e\u8a18\u9332\u3068\u91cd\u306d\u5408\u308f\u305b\u3066\u30d7\u30ed\u30c3\u30c8\u3057\u3066\u307f\u307e\u3059\u3002f1\u306et\u3092\u305a\u3089\u3057\u3001x0, y0\u3092\u518d\u8a08\u7b97\u3057\u306a\u304c\u3089\u8a08\u7b97\u3057\u307e\u3059\u3002

\r\n\r\n\u9a5a\u304f\u307b\u3069\u3088\u3044\u4e00\u81f4\u304c\u898b\u3089\u308c\u307e\u3059\u3002\r\n"}︡ ︠729d4994-cb93-47ce-8b31-b1d153899e57︠ # 米軍の増員が、初日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) ︡c7837073-e0e1-4767-b598-c2e20dccbabe︡{"html": ""}︡ ︠38834219-cb4b-435c-a8d6-2aa36145a52ai︠ %html 参考までに、日本軍の兵士の推移も表示してみます。 ︡fbd501f6-4f67-47ac-9a38-5db66b307642︡{"html": "\r\n\u53c2\u8003\u307e\u3067\u306b\u3001\u65e5\u672c\u8ecd\u306e\u5175\u58eb\u306e\u63a8\u79fb\u3082\u8868\u793a\u3057\u3066\u307f\u307e\u3059\u3002\r\n"}︡ ︠a3bdf73c-dda4-4003-9616-4662a7f14ebf︠ # 同様に日本軍の滞在人数をプロットします 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) ︡c3d5b1e9-3b3c-42d7-9a16-3db72ec0d125︡{"html": ""}︡ ︠28af3e95-673d-48e3-b696-910b50e4feca︠ # 武器の性能（兵力が同じ場合） x,y = var('x y') a=1 b=1 plot_vector_field((-a*y, -b*x), [x, 0, 1], [y, 0, 1]) ︡43dbc758-95e2-4aed-b1f9-856a298c4fa1︡{"html": ""}︡ ︠afcbe8c2-c560-4671-8311-05fcdcd4c23d︠ # 硫黄島の戦いの場合 a=0.0577 b=0.0106 plot_vector_field((-a*y, -b*x), [x, 0, 1], [y, 0, 1]) ︡19050810-3751-4302-b75e-66cceafb4e90︡{"html": ""}︡