Mathematicaで多体問題を解く

Mathematicaの強力な記号操作能力を使って多体問題を解くプログラムを作ります。
#粒子のリストを作成する
粒子は質量、電荷、位置、運動量の変数名、初期位置の変数名、初期位置、初期運動量の情報を持つとする。

##`putindex`
添字つきの変数を作成する関数`putindex`を作ります。

```mathematica
putindex[x_,i_] :=ToExpression[ToString[x] <> ToString[i]];
```
実行結果
<img width="155" alt="putindex.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/db155e0a-36ca-cacf-ee0b-7e3dd9349e1f.png">
##`makeparticles`
粒子のリストを作成する関数`makeparticles`を作成する。`Map`を使って`putindex`をそれぞれの項目に作用させ、添字つきの変数名を作成している。`Map[f,list,{-1}]`はlistの最下層の要素にfを作用させるという意味です。

```mathematica
makeparticles[n_] :=
 Table[
  Map[
   putindex[#, i]&,
   {
    m, (*質量*)
    q, (*電荷*)
    {x, y}, (*位置*)
    {px, py}, (*運動量*)
    {x0, y0}, (*初期運動量*)
    {px0, py0} (*初期運動量*)
   }, {-1}
  ],{i, n}
];
```
実行結果
<img width="440" alt="makeparticles.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/d1af0cc6-f575-49b2-b2d4-387efd858fe3.png">
##`makeparams`
粒子の質量や初期位置の値を入れるルールのリストを作成する関数`makeparams`を作成する。MapAt[f,list,{All,1}]はリストの第二レベルの1番目にfを作用させる。
`Flatten`は入れ子のリストを開封します。例:`Flatten@{a,{b,c}}`->`{a,b,c}` `f@a`は`f[a]`と同じ意味です。

```mathematica
makeparams[n_] := Flatten@Table[
 MapAt[
  putindex[#, i] &,
   {
    m -> RandomReal[{0.1, 1}],
    q -> RandomReal[{-0.1, 0.1}],
    x0 -> RandomReal[{-1, 1}],
    y0 -> RandomReal[{-1, 1}],
    px0 -> RandomReal[{-0.1, 0.1}],
    py0 -> RandomReal[{-0.1, 0.1}]
   }, {All, 1}
  ],
  {i, n}     
];
```

実行結果
<img width="786" alt="makeparams.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/b1ef17f4-6cf6-b977-c0f6-dfbc6c8cc0d4.png">

`->`は`Rule`を表す記号で、`Replace`:`/.`を使ってこのルールを好きな式に作用させれます。例：`f[a]/.a->b`->`f[b]`。`makeparams`で作ったルールを`makeparticles`で作ったリストに作用させるとさっき作ったリストの`m1`,`q1`などに数値が入ったリストができます。
最初から数値を入れていないのは、後でこの変数を組み合わせた式を作ったとき、記号式の状態で見れるようにするためです。こうするとデバッグなどで便利です。

実行結果
<img width="776" alt="particles:params.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/4f7e5987-ba82-9a48-5d9d-aa8c53e01ee0.png">

#ハミルトニアンを記述する。
ハミルトニアンとは系の全エネルギーを粒子の位置と運動量で表したものです。ハミルトニアンが記述できれば、正準方程式として機械的に解くべき微分方程式が得られます。
準備としてダサいですが、`makeparticles`で作ったリストの何番目に何の項目が入ってるかを添字に直接名前を着けることで分かりやすくします。

```
mass = 1; charge = 2; position = 3; momentum = 4; position0 = 5;
momentum0 = 6;
```
##運動エネルギー関数`T`
```mathematica
T[particles_] := 
Sum[Total@
 ((particles[[i, momentum]]^2)/(2 particles[[i, mass]])),
 {i, Length@particles}
]
```
実行結果

<img width="334" alt="T.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/901cf5eb-889e-d668-34dc-bc73ab29f3aa.png">

##相互ポテンシャル関数`V`
ここでは静電力による相互作用を考えます。（プラズマを考える）

```mathematica
V[particles_,i_, j_] := 
(particles[[i, charge]]*particles[[j, charge]])
/Sqrt[Total@((particles[[i, position]] - particles[[j, position]])^2) + 0.01];
```
実行結果

<img width="277" alt="V12.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/0f0b4918-6fef-e318-5891-4acef068e57a.png">

粒子どうしが完全に重なったときにエネルギーが無限大になる問題を防ぐため、分数の分母のルートの中に0.01を足しています。こうすると穏やかなポテンシャルになります。
<img width="488" alt="gentle_v.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/2027421b-7cd2-1375-509a-cbca84403457.png">
##ハミルトニアン`H`
```mathematica
H[particles_] :=
 T[particles] +
  Sum[If[i > j, V[particles,i,j], 0], 
  {i, Length@particles}, {j,Length@particles}
]
```
実行結果

<img width="1093" alt="H.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/cb034573-bfad-d782-5e42-d85c26c3f692.png">

#正準方程式を記述する。
まず位置と運動量の変数のリストを作ります。
<img width="557" alt="ps.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/1f2c3697-74bc-0636-ce4d-dccffbb02e94.png">
<img width="734" alt="pst.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/aa813854-f5bf-ba02-236c-aad37491383c.png">
ハミルトニアンもこれに合わせます。
<img width="841" alt="Ht.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/3b792cc8-dcf1-18e2-fae8-949d89ddff46.png">
さて、正準方程式はこういう形をしています。

```math
\begin{align*}
\dot{q}_i &= \frac{\partial H}{\partial p_i} \\
\dot{p}_i &= -\frac{\partial H}{\partial q_i}
\end{align*}
```
まずこの式の左辺は、単純に`pst`全体を微分すれば作れます。右辺に関しては、まずハミルトニアンを`pst`で微分したリストを作り、これを半分に分割し、前後を入れ替え、位置による微分には全てマイナス1を掛けることで作れます。
これはMathematicaの強力なリスト操作関数群を使ってこう書けます。

```mathematica
Thread[
 D[pst, t] == 
  Flatten[Reverse@Partition[D[Ht, {pst}], Length@ps/2]*{1, -1}]
 ]
```
この部分を理解したい人は偏微分を通常表記にしてくれる関数(下記リンクのpdconv)を使ってこの式の各部分を調べてみてください。
http://blog.wolfram.com/2011/12/15/mathematica-qa-series-converting-to-conventional-mathematical-typesetting/

```
pdconv[f_] :=
 TraditionalForm[
  f /. Derivative[inds__][g_][vars__] :>
    Apply[Defer[D[g[vars], ##]] &, 
     Transpose[{{vars}, {inds}}] /.
      {{var_, 0} :> 
        Sequence[], {var_, 1} :> {var}}]
  ]
```
<img width="579" alt="Canonical.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/2cbc6415-7d97-72f7-c4ef-c33a84620b49.png">
#`NDSolve`で微分方程式を解く。
正準方程式を立式して微分方程式を解いてくれる関数`hamiltonSolve`を作ります。
Methodに関してですが、"SymplecticPartitionedRungeKutta"というのが正準方程式に特化したソルバらしいのですが、Mathematica 8ではバグ（下記URL）があるらしく使えませんでした。https://mathematica.stackexchange.com/questions/90066/symplecticpartitionedrungekutta-shows-strange-error?noredirect=1&lq=1

```mathematica
hamiltonSolve[H_, particles_, tcond_] :=
  
  Module[{ps, pst, ips, t, te, Ht}, {t, te} = tcond;
   ps = Flatten[{particles[[All, position]], 
      particles[[All, momentum]]}];
   ips = Flatten[{particles[[All, position0]], 
      particles[[All, momentum0]]}];
   pst = Through[ps[t]];
   Ht = H /. Thread[ps -> pst];
   NDSolve[
    Join[ 
     Thread[
      D[pst, t] == 
       Flatten[Reverse@Partition[D[Ht, {pst}], Length@ps/2]*{1, -1}]
      ],
     Thread[(pst /. t -> 0) == ips]
     ], ps, {t, 0, te},
    (*Method->{
    "SymplecticPartitionedRungeKutta","DifferenceOrder"->2,
    "PositionVariables"\[Rule]Drop[pst,Length@pst/2]}*)
    (*Method->{"ImplicitRungeKutta",
    "DifferenceOrder"->4},*)
    (*Method->{"ExplicitRungeKutta","DifferenceOrder"->4},*)
    Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 4},
    StartingStepSize -> 0.01,
    MaxStepSize -> 0.1,
    MaxSteps -> Infinity
    ]
   ];
```
実行結果
<img width="781" alt="hamiltonsolve.png" src="https://qiita-image-store.s3.amazonaws.com/0/178902/70301370-c53d-09aa-2c90-4c04d6d08026.png">
#可視化
電荷を色で、質量を円の大きさで表しました。

```mathematica
DrawParticles[particles_, sol_, t_] :=
 Graphics[
  {
   {
       Blend[{Red, Blue}, Rescale[#[[charge]], {-0.1, 0.1}]],
       Circle[
        {#[[position, 1]][t], #[[position, 2]][t]}, #[[mass]]
        ]
       } & /@ particles /. sol
   , Circle[{0, 0}, 10]
   }
  , PlotRange -> {{-15, 15}, {-15, 15}}
  ]
```
```
Manipulate[
 DrawParticles[particles /. params, sol, t],
 {t, 0, 100}
 ]
```
#結果
このままだと粒子が遠くにいってしまうので、原点から半径10の位置に急激に増加するポテンシャル`V2`をつけて壁にします。壁は`Tanh`関数を使って作りました。

```mathematica
V2[particles_] := 
  Sum[Tanh[(Total@((particles[[i, position]])^2) - 100)/0.1], {i, 
    Length@particles}];
```


実行結果（10体問題）
![10body.gif](https://qiita-image-store.s3.amazonaws.com/0/178902/1da6b8c5-96a8-6ba1-7cd6-9e94d2cf20eb.gif)
青はマイナスの電荷、赤はプラスの電荷を表します。赤と青の丸がちゃんと引き合っているのが分かります。

読んでくれた方ありがとうございました。
参考ページ
https://mathematica.stackexchange.com/questions/102890/solve-motion-from-hamiltons-equations
http://eman-physics.net/analytic/hamilton.html
#コード全体
```mathematica
putindex[x_, i_] := ToExpression[ToString[x] <> ToString[i]]
makeparticles[n_] :=
  Table[
   Map[
    putindex[#, i] &,
    {
     m,
     q,
     {x, y},
     {px, py},
     {x0, y0},
     {px0, py0}
     }, {-1}
    ],
   {i, n}
   ];
makeparams[n_] := Flatten@Table[
    MapAt[
     putindex[#, i] &,
     {
      m -> RandomReal[{0.1, 1}],
      q -> RandomReal[{-0.1, 0.1}],
      x0 -> RandomReal[{-1, 1}],
      y0 -> RandomReal[{-1, 1}],
      px0 -> RandomReal[{-0.1, 0.1}],
      py0 -> RandomReal[{-0.1, 0.1}]
      }, {All, 1}
     ],
    {i, n}
   ];
mass = 1; charge = 2; position = 3; momentum = 4; position0 = 5; momentum0 = 6;
V[particles_, i_, j_] := (
  particles[[i, charge]]*particles[[j, charge]])/Sqrt[
  Total@((particles[[i, position]] - particles[[j, position]])^2) + 
   0.01];
V2[particles_] := 
  Sum[Tanh[(Total@((particles[[i, position]])^2) - 100)/0.1], {i, 
    Length@particles}];

T[particles_] := 
 Sum[Total@((particles[[i, momentum]]^2)/(
   2 particles[[i, mass]])), {i, Length@particles}]
H[particles_] :=
  T[particles] +
   Sum[If[i > j, V[particles, i, j], 0], {i, Length@particles}, {j, 
     Length@particles}] + V2[particles];
hamiltonSolve[H_, particles_, tcond_] :=
  
  Module[{ps, pst, ips, t, te, Ht}, {t, te} = tcond;
   ps = Flatten[{particles[[All, position]], 
      particles[[All, momentum]]}];
   ips = Flatten[{particles[[All, position0]], 
      particles[[All, momentum0]]}];
   pst = Through[ps[t]];
   Ht = H /. Thread[ps -> pst];
   NDSolve[
    Join[ 
     Thread[
      D[pst, t] == 
       Flatten[Reverse@Partition[D[Ht, {pst}], Length@ps/2]*{1, -1}]
      ],
     Thread[(pst /. t -> 0) == ips]
     ], ps, {t, 0, te},
    Method -> {"ImplicitRungeKutta", "DifferenceOrder" -> 4},
    StartingStepSize -> 0.01,
    MaxStepSize -> 0.1,
    MaxSteps -> Infinity
    ]
   ];
DrawParticles[particles_, sol_, t_] :=
  Graphics[
   {
    {
        Blend[{Red, Blue}, Rescale[#[[charge]], {-0.1, 0.1}]],
        Circle[
         {#[[position, 1]][t], #[[position, 2]][t]}, #[[mass]]
         ]
        } & /@ particles /. sol
    , Circle[{0, 0}, 10]
    }
   , PlotRange -> {{-15, 15}, {-15, 15}}
   ];
```
```mathematica
particles = makeparticles[10];
params = makeparams[10];
sol = hamiltonSolve[H[particles] /. params, particles /. params, {t, 100}]
Manipulate[
 DrawParticles[particles /. params, sol, t],
 {t, 0, 100}
]
```