【python】GIFアニメーションを作ってみよう★N体シミュレーションを例に添えて★

2019 9/03
【python】GIFアニメーションを作ってみよう★N体シミュレーションを例に添えて★

こんにちは、おはようございます、こんばんは、理系夫婦の夫たんたん(@rikeifufu_otto)です!

今日は、2019年の超大型連休GW(10連休!)の3日目。

折角の長期休みなのでちょっと凝ったことでもしようかと、pythonで動画(gif)を作ってみました。

まほちゃんのアイコン画像まほちゃん

それより、どっか連れってよ!


たんたんのアイコン画像たんたん

勉強が先だ!and めんど〜い

まほちゃんのアイコン画像まほちゃん

ひ、ひどい親だ・・・

たんたんのアイコン画像たんたん

・・・

目次

pythonでアニメーションを作るには

今回はpythonで動くグラフ(アニメーション)を作ってみたいと思います。

pythonのグラフ作成モジュール「matplotlib」には、アニメーションを作成するツールが2つ用意されている。

  • ArtistAnimation: あらかじめ計算結果を格納しておき、plot時には静止グラフを繋げる
  • FuncAnimation: plot時に、その都度計算を行う

今回は計算量が毎回大きそうなので、「ArtistAnimation」で作ってみました。

簡単な例

まずは非常に簡単な例を示したいと思います。

上記の結果が下記のgifです。

このプログラムでの肝の部分は17行目からのfor文のところです。

このループでは、座標(x, y)が毎回更新されます。xはループ文の指標iと同じで、yはiの2倍の値が代入されます。

この座標でplotを20行目で行い、21行目でアニメーションの元になるリストへ格納されます。

作られたリストgrapha_listは最終的に、23行目でアニメーションに加工されて、表示(24行目)やファイル保存(25行目)されます。

(なお最後のgif化する部分は「imagemagick」をインストールしておく必要がある?)

スポンサーリンク

N体シミュレーションの例

上記のアニメーション化の道具を利用して、重力相互作用のみが働く天体運動をシミュレーションしてみようと思いました。

天体が2次元空間に、初期位置と初期速度をもってばらまかれ、時間発展とともに各天体が重力相互作用のためにどのように動くかを計算、アニメーション化するのが目的です。

これが実はなかなか難しかったです。。。

ものすごく簡単な物理モデルを仮定しました。

ある1つの天体が、次の時間ステップの時の位置と速度が次の様に決定されるモデルにしました。

x方向は、

$$ x_{i+1} = x_i + v_{x,i}\Delta t + a_{x,i}\frac{(\Delta t)^2}{2} \\ v_{x,i+1} = v_{x,i} + a_{x,i}\Delta t $$

y方向は、

$$ y_{i+1} = y_i + v_{y,i}\Delta t + a_{y,i}\frac{(\Delta t)^2}{2} \\ v_{y,i+1} = v_{y,i} + a_{y,i}\Delta t $$

ここで、\( a_{x,i}, a_{y,i} \)は万有引力による他の天体からの合計力によるそれぞれの方向での加速度です。式にすると下記です。(注目している天体を\( \eta \))だとして、

$$ a_{x,i} = \frac{1}{m} \Sigma_{\zeta \neq \eta}G\frac{m\cdot m}{r_\zeta^2} \frac{x_\zeta – x_\eta}{r_\zeta} \\ a_{y,i} = \frac{1}{m} \Sigma_{\zeta \neq \eta}G\frac{m\cdot m}{r_\zeta^2} \frac{y_\zeta – y_\eta}{r_\zeta} $$

ここで\( r_\zeta \)は天体\( \zeta \)と今注目している天体\( \eta \)の距離である。また天体の質量は全て同じ\( m \)とした。

$$ r_\zeta = \sqrt{(x_\zeta – x_\eta)^2 + (y_\zeta – y_\eta)^2} $$

イメージ図:

イメージ

この計算を適当な時間間隔(\( dt \))で全天体に対して行うことで、天体の運動がシミュレーションできるだろうと考えたわけです。

コードが下記。

まずは59行目から62行目を有効にして、適当な位置と初速度を持った2天体の場合を計算してみました。結果が下記です。

2体問題

これはそれっぽい結果に見えた。

2つの天体がお互いの周り(実際には2体の重心)を楕円運動している。

3体問題も解いてみました。

3 bodies

折角なので、天体を100個に設定してみた。

すごい、みんな宇宙の彼方に飛んでいく。

この原因を調べてみるとどうやら、次ステップの計算を行う近似モデルが悪かった様だ。

今回の計算モデルは「リープフロッグ法」と呼ばれる計算手法をもっと簡単にしたので、その時の近似が良くなく、この様なことになったと思います。

もっと計算の時間間隔を細かくすれば多体でもいい感じにできかもしれません。

上記のコード内にはエネルギーの推移を示す部分も含まれるので、みてもらえばわかりますが、エネルギーが増加していきます。

(もしも根本的におかしな点を見つけましたら、ご一報ください。)

まほちゃんのアイコン画像まほちゃん

中途半端な結果を載せんなよ!

たんたんのアイコン画像たんたん

すみません、今回は「アニメーション」作成がメインですので、天体シミュレーションは多めにみてください。
「リープフロッグ法」を真面目に組めばいいのだと思いますが、君のお世話があるので、今日はここまでです!
今度必ず直します!

まほちゃんのアイコン画像まほちゃん

必ずね!

参照

[進研ゼミ]等加速度運動

N体シミュレーションの基礎

matplotlibでアニメーションを作る

関連記事

応援よろしくお願いします☆

この記事を書いた人

天文の博士号をもつ理系パパ。
3歳の娘を子育て中。
最近はダイエットに挑戦中!

コメント

コメントする

目次
閉じる