機械で「心」を作る～「AIの父」ミンスキー氏が早稲田大学で講演
コンピュータサイエンスは重要だがコンピュータそのものが重要なわけではない、とミンスキー氏は強調した。それは解けない方程式があっても、コンピュータを使うことで何が起こるかを見る事ができるからだ。それがコンピュータがなく数学しかなかった時代との違いだとミンスキー氏はいわゆる「3体問題」など力学の問題を例に出して説明した。
「3体問題」の有名なものに、ピタゴラス3体問題があります。図のような3:4:5の直角三角形の頂点上に静止させた、質量3, 4, 5の質点の運動を調べるというものです。
この問題は1893年には知られていましたが、数値的にでさえ、解決は1966年になってからです。ディアク，ホームズ『天体力学のパイオニアたち』によれば、粒子の2つが連星を形成し、第3体が高速度ではじき飛ばされるという結論は、驚くべきものだったそうです。
これを数値的に解くのは、コンピュータ無しではまあ無理ですよね。「この問題をちゃんと解くのは大変ですよ。（よく知られた微分方程式の数値解法である）ルンゲ・クッタじゃぜんぜんだめですよ」と言っていた天体力学の教授の言葉が思い出されますが、ちょうど、Mathematicaならケプラー問題の数値解を簡単に求められるという記事に対して、「でも、計算精度は大丈夫でしょうか」と訊かれたところだったので、コンピュータがあれば簡単だということを確認してみましょう。
必要な精度は状況によるので、絶対大丈夫とは言えませんが、たいていの場合には、十分よい精度で計算すると思います。有名なピタゴラス3体問題で調べてみましょう。
運動エネルギーとポテンシャルエネルギーからラグランジアンを、ラグランジアンから運動方程式を導き、解きます。解く時は、保存量（全エネルギーと運動量、角運動量）が保存されるようにするといいでしょう。
(* Mass *)
m = {3, 4, 5};
(* Coordination *)
q = {{x1[t], y1[t]}, {x2[t], y2[t]}, {x3[t], y3[t]}};
(* Kinetic energy *)
V = Sum[m[[i]] Sum[D[q[[i, j]], t]^2,
{j, 1, Length@q[[i]]}], {i, 1, Length@q}]/2;
(* Potential energy *)
U = -Total[
Sum[Sum[m[[i]] m[[j]]/Sqrt[Total[(q[[i]] - q[[j]])^2]],
{j, i + 1, Length@q}], {i, 1, Length@q - 1}]];
(* Lagrangian *)
L = V - U;
(* Momentum *)
p = Table[D[V, D[q[[i, j]], t]], {i, 1, 3}, {j, 1, 2}];
(* Angular momentum *)
Mz = Table[
Cross[Append[q[[i]], 0], Append[p[[i]], 0]][[3]], {i, 1, 3}];
(* Equations of motion *)
eqns = Table[D[D[L, D[q[[i, j]], t]], t] == D[L, q[[i, j]]],
{i, 1, 3}, {j, 1, 2}];
initialCondition = {{
x1[0] == 1, y1[0] == 3,
x2[0] == -2, y2[0] == -1,
x3[0] == 1, y3[0] == -1},
Table[D[q[[i, j]], t] == 0 /. t -> 0,
{i, 1, 3}, {j, 1, 2}]} // Flatten;
tmax = 70;
s = NDSolve[
{eqns, initialCondition} // Flatten,
{q, D[q, t]} // Flatten,
{t, 0, tmax},
Method -> {"Projection",
Method -> "ExplicitRungeKutta",
"Invariants" -> {V + U, Total[p][[1]], Total[p][[2]], Total[Mz]}},
StartingStepSize -> 0.01];
plots = Table[
Show[
ParametricPlot[Evaluate[q[[1]] /. s], {t, 0, time},
PlotStyle -> Red, PlotRange -> {{-5, 5}, {-5, 5}}],
ParametricPlot[Evaluate[q[[2]] /. s], {t, 0, time},
PlotStyle -> Blue, PlotRange -> {{-5, 5}, {-5, 5}}],
ParametricPlot[Evaluate[q[[3]] /. s], {t, 0, time},
PlotStyle -> Green, PlotRange -> {{-5, 5}, {-5, 5}}]
], {time, 0.01, tmax}];
ListAnimate[plots]
Wolfram CDF Playerがインストールされていれば動きます（UMMでも動きます）。
この結果は、Burrau’s problem of three bodiesで紹介されているものとだいたいあっています。（リンク先にはこの問題を数値的に解いたSzebehelyの論文があります。）
最後に、保存量が実際に保存されているかどうかを見ておきましょう。
Related posts: