先日ツイートした「拍と拍子の概念を数学で表現する」の検証を行いました。一応計算はできたのですが、拍子の概念を定義するのは非常に難しいということが分かりました。 pic.twitter.com/9MZpi5mHBW

— 博士課程よしだドクタ (@jikkenyametatta) January 5, 2014

「音楽的センスが皆無だが数学と物理なら話が通じる。拍と拍子の概念をどう伝えればよいか?」

ふむ。早い話が宇宙人を想定すればいいのですね。
(私の脳内宇宙人はいつも死ぬほど数学と物理に長けている設定になっています。)

先のツイート内では拍子テンプレート関数なるものを用意していました。これなら拍子以外の要素(コード進行とか?)も考えられそうですが、あまり複雑にすると楽曲のテンポに対しては結果が安定しなくなりそうです。

私ならどうするか。拍子だけならフーリエ変換でもいけそうです。
音をエネルギーか何かの時系列データで表現して、それをフーリエ変換して波数空間で見てあげれば、拍や拍子のピークが見えてもよさそうです。拍の半分の周波数でピークが出れば2拍子、1/3倍の周波数でピークが出れば3拍子です。

と、予想を書いてみるのは簡単ですが、本当にこんなことになるのでしょうか。

というわけで試しにやってみました。
関連文献をあたれば何かしら見つかりそうですが…

続きを読む →

「進捗・どう・です・か」をランダムに表示し「進捗どうですか」が完成したら煽ってくるプログラム

コメントをどうぞ

どうしてこんなものを作ってしまったんだろう。

進捗どうどう進捗ですか

というわけで、今回のネタは無限の猿定理です。直接の元ネタはTwitterですが、下ネタ全開なのでリンクを開く際はご注意ください。

概要

タイトルにもあるように「進捗・どう・です・か」をランダムな順序で表示し、「進捗どうですか」が完成したところで適当にクエスチョンマークをつけて去っていくプログラムです。

特に意味もなく文字カウントします。おおよそ数十文字から1000文字程度で目的を成し遂げて去っていきます。

生成された煽り文句たち

煽り文句も煽り以外も結構バリエーションがあるものです。

「進捗どう進捗どう進捗どう」

焦っているんですかね。

「進捗ですです進捗どうです」

激しい進捗があった人のセリフでしょう。

「どう進捗ですか」

進捗を伝えてからの冷酷な迎え撃ち。

「進捗かどうかどうですか」

自信がなくなってきた。

「どうかどうか進捗進捗進捗」

締め切り間近。もうやばそう。

続きを読む →

マンデルブロー集合へ飛び込む

1件のフィードバック

先日のマンデルブロー集合の10の200乗倍で作成したプログラムがあるので、今回はマンデルブロー集合の様々な地点へ飛び込んだ動画を作ってみました。

前回の記事を見ていない方はそちらからどうぞ:

[フラクタル] マンデルブロー集合 10の200乗倍への挑戦

今回訪れる地点は計5箇所です。どの点も拡大していく最中に集合全体と同じ形が現れます。フラクタルらしくて面白いと思います。

具体的には以下の点です。緑は前回の場所。それぞれ個性的な場所だと思います。

Mandelbrot Set with Circles (Picture from Wikimedia)

Mandelbrot Set with Circles (from Wikimedia)

この集合はフラクタルなので自己相似の図形ではあるのですが、ただ単純にどこでも同じ形を持っているわけではありません。拡大する場所を変えると景色が一変してしまいます。

肝心の動画はこちら。ニコニコ動画はこちら。

続きを読む →

院試勉強しながらついTwitterをしてしまうあなたに「院試まで残り何秒bot」作ってみた

1件のフィードバック

お勉強の合間に邪悪なbotを作成。ちょっぴり内輪向けです。

こちらです → @inshi_sister

東大工学系院試まであと85日と42319秒だよ。

— Inshi_Counter (@inshi_sister) June 1, 2013

最近徐々に身の回りで院試勉強をしている人が増えてきています。
お勉強というのは人によっていろいろなスタイルがあります。私は図書館くらいの空気がちょうどいいのですが、人が大勢いる環境で捗る人も意外と多いようです。

そういう方々が集中するために使っているのか、よく見る(煽り合っている?)のが「院試カウントダウンタイマー」やそれに類するカウントダウン系焦らせグッズです。

私がこんなの使いながら勉強したら発狂するんじゃないかと思いますが、まあいろいろな人がいるものです。

一方で図書館などでもよく見るのがTwitterを見ながらのお勉強です。なるほどTwitterというのはひとりで勉強中でも連帯感みたいなものを提供してくれます。一方でこいつはかなりの時間泥棒で、どうやったらTwitterから離れられるか四苦八苦している人も多いようです。

今回はその2つを単純に足して、Twitter上で院試までの残り時間を主張するbotを作ってみました。

続きを読む →

[フラクタル] マンデルブロー集合 10の200乗倍への挑戦

2件の返信

全体は部分の総和に勝る。 – アリストテレス

マンデルブロー集合、10の200乗倍、グーゴル倍のグーゴル倍に挑戦しました。

Youtubeのほう、最後がぶつ切りになってしまっていますが許してください。
ニコニコ動画はこちら。

(追記)もっと画質のよい動画が欲しい場合はGoogle Driveからダウンロードして下さい。1GB弱あります。

現実世界だと陽子や電子の半径が10の-15乗メートル、観測可能な宇宙の範囲でもたったの(?)10の27乗メートルしかないため、これほどの倍率はまず見かけることはありません。

ニコニコ動画で何度も見返した動画というのはそれほど多くないのですが、その一つに「【超高倍率】マンデルブロー集合【10の100乗倍】 (hoho様)」というものがあります。今回はその動画に感謝しつつ、記録更新を狙いました。
先の動画ではニコニコ最高記録と記されているので、それが正しければ今回がニコニコ新記録です。(logスケールで)一気に2倍まで来たことになります。

続編はこちら: マンデルブロー集合へ飛び込む

続きを読む →

[絶対に儲かる賭け?] サンクトペテルブルグのパラドックスを律儀に試す

コメントをどうぞ

「サンクトペテルブルグのパラドックス」は、現実的にはあまり儲からなそうな賭けの期待値が、計算上無限大に発散するとされる有名なパラドックスです。

パラドックスを知らない人のために賭けの要点を書いておくと、以下のようになります。

「コインを投げ続け、表が出たらゲーム終了。表が出たのが1回目なら1円獲得、2回目なら2円、3回目なら4円…というように、裏が1回出るたびにもらえる金額が倍になる。いくらの賭け金ならこの賭けに参加するべきか。」

儲け額の期待値Wは以下のように計算されます。1/2で1円、1/4で2円、となるので…
$W=\sum_{k=0}^{\infty}\left(2^{k-1}\frac{1}{2^k}\right) = \frac{1}{2}+\frac{1}{2}+\frac{1}{2}+\dots = \infty$

ここまで見てだまされないぞと思った人は、しっかり悩んでみてください。

Taken by xJason.Rogersx, Flicker.

先日友人との間で話題になって、実際のところおおよそどんな具合で儲かるのかやってみようということでさくっと書いてみたものです。確率についてなんやかんやと議論するのも楽しいですが、ゴリ押しでデータを集めるのもまた面白いものがあります。あと探しても意外とそういう資料が見当たらなかったので。

ソースコード

いつもどおりC++です。C++11の機能を使っていますので、例えばgccとかでコンパイルしたい場合は -std=++11 とか -std=c++0x とかのオプションをつけてください。

#include <ctime>
#include <iostream>
#include <random>
#include <functional>

int main(int argc, char** argv){
  if(argc != 2){
    std::cerr << "Wrong argument!" << std::endl;
    return -1;
  }
  long long int loopNum = 1;
  for(int i = 0; i < atoi(argv[1]); i++){
    loopNum *= 2;
  }
  std::cout << "Loop: " << loopNum << std::endl;

  time_t timer;
  time(&timer);
  auto rand = std::bind(std::uniform_int_distribution(0,1),
      std::mt19937(static_cast<unsigned int>(timer)));
  long long int total = 0;
  for(long long int trial = 0; trial < loopNum; trial++){
    long long int earn = 1;
    while(rand()){
      earn*= 2;
    }
    total += earn;
  }
  std::cout << "Earned: "
    << static_cast<double>(total)/static_cast<double>(loopNum)
    << std::endl;
  return 0;
}

特記するようなものもありません。コマンドライン引数から2の何乗回計算するかを決めて、律儀に試行を繰り返して平均獲得金額を求めています。

カウントがlong long int型なので2の50乗回程度までは問題ないと思いますが、そのはるか手前の2の30乗前後(数億程度)で計算時間がきつくなってくると思います。
まあ乱数も毎回1ビットしか使っていないしシングルスレッドだしで、高速化する余地はまだ残っています。

結果

大量に試行した結果は以下のようになります。あくまで愚直に。力強く。

横が1セットあたりの試行回数、縦が1回あたりの平均獲得金額です。両方logスケールなのに注意。
(表が出るまでを1回の賭けとし、1セットでn回行うと定義しています。平均獲得金額は儲け額をnで割った値。)

不思議なグラフです。まず、回数を増やすと平均獲得金額は順調に上がっていくことがわかります。1,000,000,000回以上やるつもりなら、1回あたり10円以上儲かることがほぼ確実になります。このままグラフを右に伸ばしていったらどうなるでしょう。

もうひとつ面白いことがあって、それはばらつきが減らないことです。例えば右上にある点、だいたい10,000,000回で平均獲得金額が10,000円を超えていますが、つまりこの1セットで単純にトータル一千億円以上儲けているということです。とんでもない大当たりです。

胴元と挑戦者、どちらの席につきたいですか?

今回は横軸がlogスケールで均等になるようにサンプリングしています。離散的な値しかとらないので左側はすかすかに見えますが、実際は点が重なっているだけです。
また、左上から右下へのぼんやりとした斜線が大量に見られます。これは「1回大当たり」を成し遂げたときの線になっています。回数を増やすと同じ大当たりでも「1回あたりの獲得金額」は減っていくので右下がりになるという具合です。その代わり「もっとすごい大当たり」をひく確率も出てきます。

モンテカルロ法などでは試行回数を増やすと値が落ち着くのですが、この問題はおもしろいことにいくら試行回数を増やしてもばらつきが収まらないようになっています。よくできています。

儲け額の分散も出してみようと一度書きかけたのですが、ちょっと考えてからこの問題に対して無駄な試みであることがわかりました。

コンピュータはしばしば想像力を働かせる道具として役に立つのですが、今回のように神通力の効かないヤバそうな問題ではちょっとわけが違ってきます。こんな問題でも扱えてしまう数学の力は偉大です。

結果の意味とこの賭けの本当の期待値についてはいろいろ考えてみて下さい。

運動方程式の数値計算と速度ベルレ法の威力 (後半)

1件のフィードバック

前半の補足です。前半で紹介したシミュレータ云々について。

この粒子の振る舞いについて

「なんで星の運動を模してるのに楕円軌道じゃないんだ?」と思った方はとても勘のいい方です。

これは2次元の運動をどう考えるかによるものです。
3次元空間の平面上に粒子を配置させて運動させた場合、重力は距離の2乗に反比例しますが、
3次元空間に無限長の棒を(平行に)置いて運動させた場合、重力は距離の1乗に反比例します。

これはまあ真面目に積分してもそうなりますし、電磁気学を勉強したことのある方は雰囲気で理解できるかと思います。後者は宇宙空間では仮定が不自然すぎるだろうと思われるかもしれませんが、万有引力が3次元空間で逆2乗の法則に従っているとするなら、2次元では逆1乗になるのが自然だとも考えられます。

まあ、今回はそこまで哲学的な話ではなく、単に逆1乗のほうが見栄えが良かったからです。

どう見栄えに影響するかですが、これは無限遠点での挙動の違いを考えればだいたいわかります。位置エネルギーは力の積分(の-1倍)で与えられますが、これより力が逆2乗、逆1乗の場合のそれぞれについて位置エネルギーは逆1乗、logとなります。つまり前者ならエネルギーが無限遠である値に収束するのに対して、逆1乗なら際限なく増えていきます。

つまり、逆2乗ではひとたび粒子が第三宇宙速度に達してしまうと二度と戻ってこないのに対し、逆1乗なら少なくともいずれは引き戻されてくるという大きな違いがあります。見た目に粒子がどこかに飛んでいってしまうのは寂しいですね。

よく「この世が2次元だったらいいのに」という言葉を聞きますが、そんなことをしたら、上記の理由から世界中のものがある一点に集まってしまってつまらない世界になるだろうなあと思っています。

なお、粒子があまり近づきすぎると計算上の問題からヤバいことが起きるので、引力を打ち消すだけの反発力を与えてあります。
現実だと衝突しますが、無理にエネルギーを保存させようとして星が完全弾性衝突をするのもあまりリアリティがないですね。

以前の3重振り子を使わなかった理由

これは速度ベルレ法の適用範囲に限界があるからです。(下は剛体3重振り子の動画)

シンプレクティック法は、影のハミルトニアンと呼ばれる量を厳密に保存する運動を記述します。この際、位置 $q_{t+1}$ を計算する際に $q_{t+1}$ のときのハミルトニアン(の偏微分)が必要になります。これは普通に考えると右辺と左辺がからまってさくっと解けません。
しかし都合のいいことに、単純な場合ハミルトニアンHは位置qと運動量pの足し算に分離でき、このときハミルトニアンの偏微分の形からqがなくなります。その場合の解法のひとつが速度ベルレ法となるわけです。

以前の3重振り子では、ハミルトニアン(というよりラグランジアン)の計算に3つのθという一般化座標を使っていました。このためハミルトニアンの分離が行えなくなり、(実際にはできるかもしれないけど私はその数学を持ち合わせていないので、)速度ベルレ法をはじめとした陽解法が使えない状態に陥っていました。

こういった場合には陰解法のガウス・ルジャンドル法などを使って解くことができます。

なんとなく当然のような気もしますが、よくよく考えてみるとこの仮定は結構不思議なことです。分子動力学法などではたいていこの仮定が成り立ってしまいます。

ソースコード

ちょっとだけかいつまんで書いておきます。C++です。

まずは系の状態の定義。数値積分に関しては位置と速度のペアのリストさえあればいいのだけど、肝心な力の計算ではどの値がどの物理量に対応しているかを知らなければならず、もどかしい。

かといって値を参照して保持しておくほどの大げさなものでもない。そこで、どうせ数値積分では四則演算くらいしか使わないので、粒子のvectorに必要な演算子のみ定義した構造体を持っておきます。

	struct partState{
		double x, y;
		double vx, vy;
		partState()	: x(0.0), y(0.0), vx(0.0), vy(0.0){}
		partState(double _x, double _y, double _vx, double _vy)
			: x(_x), y(_y), vx(_vx), vy(_vy)
		{}
	};
	struct partSystem{
		std::vector<partState> pss;
		partSystem(int size){
			pss = std::vector<partState>(size);
		}
	};
	const partState operator+(const partState& lhs, const partState& rhs);
	const partState operator-(const partState& lhs, const partState& rhs);
	const partSystem operator+(const partSystem& lhs, const partSystem& rhs);
	const partSystem operator-(const partSystem& lhs, const partSystem& rhs);
	const partState operator*(const double lhs, const partState& rhs);
	const partSystem operator*(const double lhs, const partSystem& rhs);

力の計算はもうほんとに全ペアで単純に逆1乗で計算してるだけで、高速化のかけらもありません。
ガチで分子動力学法などをやるときには今回のO(n^2)のようにならない工夫が必要です。

		partSystem calcForce(const partSystem& _ps){
			partSystem ans(size);
			for(int i = 0; i < size; i++){
				for(int j = i+1; j < size; j++){
					double dx = _ps.pss.at(i).x - _ps.pss.at(j).x, dy = _ps.pss.at(i).y - _ps.pss.at(j).y;
					double r2 = dx*dx + dy*dy;
					double r = std::sqrt(r2);
					if(r > r_threshold){
						//double r12_e3 = r2/0.0001;
						//r12_e3 = r12_e3*r12_e3*r12_e3*r12_e3*r12_e3*r12_e3;
						double f = -g*masses.at(i)*masses.at(j)/r;// + g2*masses.at(i)*masses.at(j)/r12_e3;
						double fx = f*dx/r, fy = f*dy/r;
						ans.pss.at(i).vx += fx;
						ans.pss.at(i).vy += fy;
						ans.pss.at(j).vx -= fx;
						ans.pss.at(j).vy -= fy;
					}
				}
				ans.pss.at(i).x = _ps.pss.at(i).vx;
				ans.pss.at(i).y = _ps.pss.at(i).vy;
			}
			return ans;
		}

反発力を与える際、分子動力学チックに高次の多項式の反発力を与えようとした名残があります。結局そこまでする必要もなく削除。かわりに距離にしきい値を導入します。
エネルギーの計算(表示用)でも不整合がないようにしておく必要があります。

あとは今回のテーマだった数値積分。ソースはとっても単純。

		void move_euler(const double dt){
			auto psDiff = calcAccel(ps);
			ps = ps + dt * psDiff;
		}
		void move_heun(const double dt){
			auto k1 = calcAccel(ps);
			auto k2 = calcAccel(ps + dt*k1);
			ps = ps + dt/2.0*(k1+k2);
		}
		void move_runge_kutta(const double dt){
			auto k1 = calcAccel(ps);
			auto k2 = calcAccel(ps + dt/2.0*k1);
			auto k3 = calcAccel(ps + dt/2.0*k2);
			auto k4 = calcAccel(ps + dt*k3);
			ps = ps + dt/6.0*(k1+2*k2+2*k3+k4);
		}
		void move_velocity_velret(const double dt){
			auto k1 = calcAccel(ps);
			auto p2 = ps;
			for(int i = 0; i < size; i++){
				p2.pss.at(i).x += dt*ps.pss.at(i).vx + dt*dt/2.0*k1.pss.at(i).vx;
				p2.pss.at(i).y += dt*ps.pss.at(i).vy + dt*dt/2.0*k1.pss.at(i).vy;
			}
			auto k2 = calcAccel(p2);
			for(int i = 0; i < size; i++){
				p2.pss.at(i).vx += dt/2.0*(k2.pss.at(i).vx + k1.pss.at(i).vx);
				p2.pss.at(i).vy += dt/2.0*(k2.pss.at(i).vy + k1.pss.at(i).vy);
			}
			ps = p2;
		}

速度ベルレ法だけ、位置と速度を明示的に扱っているのがわかります。

運動方程式の数値計算と速度ベルレ法の威力 (前半)

3件の返信

前々回の振り子の欄でちょっと触れた運動方程式の数値計算の話の続きです。

運動方程式は立てるだけでなく、実際に方程式を解かないとシミュレーションは行えません。

一番丁寧な運動方程式の解き方は解析的に解くことですが、ほとんどの場合うまい方法がないことがわかっています。例えばN体に万有引力が働いて運動する系ではNが3以上になると解析的に解けなくなることが知られています。太陽と地球と月でもうオーバーです。

もっと一般的に扱いたいときは、コンピュータが強力な武器になります。微分方程式を差分方程式に近似してしまえば計算機で簡単に扱えるので、後は時間の刻み幅を十分に小さくとってあげれば微分方程式に近い解が得られるというわけです。

じゃあそれで話は終わりかというとそうでもなく、この差分方程式にする段階がなかなかうまくいかないわけです。それに伴っていろいろな手法が提案されています。今回はその話です。

今回作ったものは互いに引力を及ぼし合う粒子たちのシミュレーションです。

イメージしているのは銀河系的なアレですが、まあ数値積分の例で作っただけなので本気で似せているわけではありません。私は中学校で習う電子の描像のイメージに似ているなと思いました。
IAEAのマークにもあるアレです。
中央に質量の大きめな粒子をひとつ、あとは小さな粒子をばら撒いています。小さな粒子同士も相互作用します。

大事なのはこの系では力学的エネルギー・運動量・角運動量が保存するということです。シミュレーションでエネルギーが保存していなかったら、それは何かおかしいことが起こっているぞ、というわけです。

可視化の部分は振り子の流用ですね。数字キーでシミュレーション手法を変えたり、初期化したり軌跡のオンオフを行ったりします。

後述するオイラー法、ホイン法、ルンゲ・クッタ法、速度ベルレ法を実装しています。

いつものようにプログラム: particle
ソースコード: GitHub

微分方程式を解く際、おそらく最も単純な差分化はオイラー法です。

ニュートンの運動方程式は2階の微分方程式ですが、速度v自体も変数と置いてしまうことで1階の微分方程式に置き換えることができます。以下の式で左辺をx、右辺をfと置いてしまえばいいわけですね。
$\frac{d}{dt} \left( \begin{array}{cc} x\\ v\\ \end{array} \right) = \left( \begin{array}{cc} v\\ f(t,x,v)\\ \end{array} \right)$

オイラー法では以下のような差分化を行います。⊿tは十分に短い時間をとります。
$x_{t+\Delta t}=x_t+\Delta t\cdot f(t,x)$

見た目にわかりやすいため学習用として頻繁に登場しますが、正確なシミュレーションをしたい場合向きではありません。誤差がとても大きいからです。

誤差が大きい理由は2点の補間を直線、つまり1次式近似にしているからで、それなら2次や3次の項を追加して補間しようというアイデアがあります。実際これはうまく行きます。
よく使われるものに4次多項式で補間するルンゲ・クッタ法というのがあります。
$k_1=f(t,x_t)\\ k_2=f(t+\frac{1}{2}\Delta t,x_t+\frac{1}{2}\Delta t\cdot k_1)\\ k_3=f(t+\frac{1}{2}\Delta t,x_t+\frac{1}{2}\Delta t\cdot k_2)\\ k_4=f(t+\Delta t,x_t+\Delta t\cdot k_3)\\ x_{t+\Delta t} = x_t + \frac{1}{6}\Delta t(k_1+2k_2+2k_3+k_4)$

tからt+⊿tまでの4点で傾きを計算して、適当に混ぜることで4次の精度を出しています。式中の諸係数を求めるのにはそこそこの手間がかかります。
次数を上げるほど多項式補間の精度は上がりますが、そうすると計算の安定性の議論から時間刻み⊿tをかなり小さくとらないといけなくなってしまうこと、また計算コストが上がるのでそれなら次数を上げずに時間刻みを減らしたほうが都合がいい場合が多いことなどがあり、結果としてこの4次のルンゲ・クッタ法はよく使われます。

1次、2次、4次の各手法で時間経過とエネルギーのずれをプロットしたものが以下になります。

1次のオイラー法は目に見えて粒子が遠ざかっていってしまいます。2次(ホイン法)と4次(ルンゲ・クッタ法)はシミュレーション結果を見た雰囲気では似ているものの、エネルギーを見るとホイン法はだんだんとずれていくことがわかります。
もちろんオイラー法にも利点はあって、軽くて実装が楽なので、正確さにあまりこだわりがない場合にはいちばん手軽な解法になります。

陰的解法と言われる手法群もあるのですが、一般に複雑なので今回は省略します。

速度ベルレ法

ここで、今回メインで書きたかった手法「速度ベルレ法」を紹介したいと思います。2次の精度をもつ手法です。
2次ならルンゲ・クッタ法より精度が悪いのかというとそうでもありません。先のグラフの倍率を上げたものに速度ベルレ法を追加してみます。

ルンゲ・クッタ法に対してぶれは大きいですが、ホイン法のようにエネルギーが変化していくことはありません。しかしまだこれではルンゲ・クッタ法のほうが頼りがいがあるように見えます。
面白いのは長時間の振る舞いです。このまま横軸のスケールを100倍にしてみます。

ルンゲ・クッタ法もこの長時間では計算誤差が積もってエネルギーがずれていくのに対して、速度ベルレ法はほぼ初期の振る舞いを保ったままとなっています。短期間のぶれは大きいのに、長期間でずれていかない不思議なグラフになっています。

速度ベルレ法の式は以下のようになります。位置と速度を分けて扱うのが今までの手法と違います。

$x_{t+\Delta t}=x_t+\Delta t\cdot v_t+\frac{1}{2}{\Delta t}^2\cdot f(t,x_t)\\ v_{t+\Delta t}=v_t+\frac{1}{2}\Delta t(f(t,x_t)+f(t+\Delta t,x_{t+\Delta t}))$

なぜこのような振る舞いを見せるのかにはわけがあって、この速度ベルレ法はシンプレクティック性というとても重要な性質を持つように作られています。

ここからちょっと難しい話ですが、ニュートンの運動方程式は時代を追ってラグランジュ形式、ハミルトン形式へと再公式化されていきます。このハミルトン形式で表現したとき、運動は位相空間上で体積を保つ性質があることがわかっています。
速度ベルレ法も含まれるシンプレクティック積分法と呼ばれる手法では、この位相空間での体積を保存する性質を有していて、結果的に長時間でのシミュレーションで(ある種の)誤差が積もっていく現象が起きなくなっています。

ハミルトン形式におけるとある保存量を厳密に保存しており、時間反転対称性(終了状態から反対向きにシミュレーションを行うと開始状態に行くことができる)、エネルギーなどの保存量を保存する、といった今までの手法にはなかった性質があります。

速度ベルレ法は、私が研究でよく扱う分子動力学法でもよく使われています。系を長時間観察するときなどにこの性質がとても有効だからです。また分子動力学法では計算コストのネックが力の計算になっており、1回に多段で力の計算を行うルンゲ・クッタ法などは相性があまりよくありません。

式自体は単純なため「分子動力学法で使われる数値計算は単純」という表現をよく耳にします。しかしその構成はよく考えられていて、裏側には先人の知恵がたっぷり詰まっていることは確かです。

運動方程式ひとつをとっても様々な数値計算法がありますが、どの積分法が正しいかといった議論はあまり有益ではなく、私のイメージでは「どれも間違っているが、間違え方が違うので場合に応じてどうしても譲れないものに応じて選ぶのが良い」としか言えないかなと思っています。

例えばゲームの物理演算などでは、そこらへんの物が綺麗な放物線を描かなくてもそんなに気になりませんが、プルプル震えたり突然とんでもない勢いで吹っ飛んでいったりすると見た目からしてヤバい、ということになるわけです。それなら計算精度より安定性と計算コストを重視しようとなるわけです。

では、現実のこの世界は、どんな手法で運動方程式が解かれているのでしょう。これはちょっと危険な問いで、そもそも現実世界が運動方程式に厳密に従っているのか、宇宙は連続体なのか、シミュレーション可能なのか、といったことからあれこれ考えなければいけません。
まあでも、私が宇宙をひとつ作っていいと言われたら、たぶんオイラー法を採用します。一番単純だし、世界がそんなに複雑な式でできているとは考えたくないからです。

後半では、今回のシミュレーションで作った系について(なぜ前回の振り子のままではだめだったのか、星の運動としてこの動きの胡散臭さは何か、ソースコード云々)の説明をしたいと思います。
速度ベルレ法と言っても万能ではないのです。

[作ってみた] [プログラミング] 剛体3重振り子

3件の返信

以前ここで書きかけた剛体振り子のプログラムですが、いろいろと紛失してしまったので新しく作りました。

ソフトのダウンロードはこちらから pendulums

ソースコードはgithubからどうぞ https://github.com/So-Takamoto/triple_pendulums/
手元の環境はVisual Studio 2010です。

(以前の記事: [昔話] 振り子のシミュレーションをした話)

例によってC++, OpenGL+freeglutで作成しています。
キーの1,2,3で微分方程式を解く方法を変更でき、iで初期化、lで軌跡の表示切り替え、qで終了です。

質量や重力加速度、初期条件、各振り子の長さを変えたい場合はソースコードをいじってみてください。
64bit版でコンパイルを行いたい場合はfreeglutのライブラリを64bit版に変更してください。今回は軽いのであまり必要はないかと思いますが。

例によって動画を用意しておきます。

なんでこんなシミュレーションをしようと思ったかの話は以前書いた(早い話がこういうおもちゃが大学に置いてあった)ので、今回は原理とか内容の紹介メインにしたいと思います。

今回扱う系は平面上の棒3個なので、そのままで言えば3自由度×3個=9自由度の系なのですが、各ジョイントが2自由度分の拘束を持つため結局本来の自由度は3になります。

ひとつのやり方は、これらのジョイントにかかる力を陽に扱う方法です。利点は自由度を実座標と対応付けられるため意味がわかりやすいことと、比較的機械的な処理で立式できることです。
これらジョイントにかかる力は幾何学的な関係から導出することもできますし、今回のような単純なリンク機構であればエンドエフェクタ(一番先の振り子)から芋づる式に数値解析の範囲で導出することもできます。ロボットハンドなどの制御ではこの辺りの話が動力学計算としてたくさん出てきます。

一方で今回のような系ならそこまで面倒な話をしなくても、解析的に運動方程式を求めることができます。こういう系ではラグランジュ形式を出発点にすると便利です。ジョイントに働く力を陽に扱わずに、自由度の数だけ運動方程式を考えればいいからです。

ラグランジュの運動方程式の出発点は、ラグランジアンL(=T-U)の表現です。運動エネルギーTは重心の並進運動と回転運動から、位置エネルギーUは重心の位置から求まります。

3つの剛体を根本から順に0,1,2と置いて(0から始めるのはC言語使いの癖です)、それぞれ棒の角度を $\theta_{i}$ 、長さを $l_{i}$ 、質量を $m_{i}$ 、根元側のジョイントからの重心の相対位置を ${\bf x}_{ri}$ と置きます。3自由度はここでは3つのθに対応します。

このとき、相対位置と相対速度について、

${\bf x}_{ri} = (\frac{1}{2}l_{i}\sin{\theta_{i}},\frac{1}{2}l_{i}\cos \theta_{i})\\ {\bf v}_{ri} = \frac{\partial}{\partial t}{\bf x}_{ri} = (\frac{1}{2}l_{i}\dot{\theta}\cos{\theta_{i}},-\frac{1}{2}l_{i}\dot{\theta}\sin{\theta_{i}})$

また、絶対位置 ${\bf x}_{i}$ について、

${\bf x}_{0}={\bf x}_{r0}\\ {\bf x}_{1}=2{\bf x}_{r0}+{\bf x}_{r1}\\ {\bf x}_{2}=2{\bf x}_{r0}+2{\bf x}_{r1}+{\bf x}_{r2}$

速度 ${\bf v}_{i}$ も上式を微分するだけです。欲しいのは $\frac{1}{2}m_{i}|{\bf v}_{i}|^{2}=\frac{1}{2}m_{i}{\bf v}\cdot{\bf v}$ なので、展開して内積を計算すればそれぞれ導出できます。

回転のエネルギーについては $\dot{\theta}$ と慣性モーメントから計算できます。長さlの棒の慣性モーメントIは $I=\frac{1}{12}ml^2$ なので、回転のエネルギーについては $\frac{1}{24}ml^2\dot{\theta}^2$ と求まります。

なんやかんやで計算すると以下のようになります。

$T=(\frac{1}{6}m_{0}+\frac{1}{2}m_{1}+\frac{1}{2}m_{2})l_0^2\dot{\theta_0}^2\\ +(\frac{1}{6}m_{1}+\frac{1}{2}m_{2})l_1^2\dot{\theta_1}^2\\ +(\frac{1}{6}m_{2})l_2^2\dot{\theta_2}^2\\ +(\frac{1}{2}m_{1}+m_{2})l_0l_1\dot{\theta_0}\dot{\theta_1}\cos(\theta_0-\theta_1)\\ +(\frac{1}{2}m_{2})l_0l_2\dot{\theta_0}\dot{\theta_2}\cos(\theta_0-\theta_2)\\ +(\frac{1}{2}m_{2})l_0l_2\dot{\theta_1}\dot{\theta_2}\cos(\theta_1-\theta_2)$

汚いようですがθに関しては対称的な形をしています。

同様に位置エネルギーに関しても

$U=(\frac{1}{2}m_0+m_1+m_2)gl_0\cos\theta_0+(\frac{1}{2}m_1+m_2)gl_1\cos\theta_1+(\frac{1}{2}m_2)gl_2\cos\theta_2$

ラグランジアンLは先に書いたようにL=T-Uとして求まります。あとはラグランジュの運動方程式

$\frac{d}{dt}\frac{\partial L}{\partial \dot{\theta_i}}-\frac{\partial L}{\partial \theta_i}=0$

を解けば $\ddot{\theta_i}$ について求まるという算段です。ただし、上の式を解こうとすると $\ddot{\theta_i}$ についての3元連立方程式になるので、微分方程式を解く前に適当な方法で変数を独立させておく必要があります。3変数なら一応解析的にも書き下すことができますが、大変面倒なので今回はガウス・ジョルダンの消去法を使いました。

ラグランジュ形式は便利ですが、最後に出てくる式がいまいち意味付けしにくいという不思議なところがあります。もっとも世の中の物理というのは別に人間にわかってもらうために運動しているわけではないのですが。
そもそもラグランジアンL自体なんの値なのかと考えだすと悩みどころはどんどん増えていきます。教科書には「意味のある値ではない」ときっぱり書いてあったりします。清々しいですね。

まあそれでもハミルトン形式まで持っていけばハミルトニアンHに意味を持たせることができるのですが。

運動方程式が解けたら次は微分方程式の数値計算の出番です。

微分方程式を数値計算で解く方法は数多くあります。今回のプログラムでは3種類から選べるようにしています。すべて陽解法で、それぞれ1次、2次、4次の精度を持っています。

おそらく最もシンプルでわかりやすい方法がオイラー法で、離散化した点の間を直線で補間する方法です。直線なので1次の精度がある、と表現されます。
オイラー法は精度が悪いため一般には使われません。左上のエネルギーの値のぶれが大きいことからもわかると思います。

より高い次数まで精度を持たせた方法があり、例えば4次のルンゲ・クッタ法が有名です。今回は2次のホイン法も用意しました。エネルギーの値を見れば違いがわかると思います。

ただし、たとえ4次精度を持っていたとしても長期間のシミュレーションを行うと次第にエネルギーが変化していってしまいます。長期間のシミュレーションでこういう性質があると結構問題で、例えば太陽系のシミュレーションをしていたら計算誤差で惑星が吹っ飛んでいったとか、分子の運動をシミュレーションしていたら計算誤差で勝手に沸騰した、なんていうよくわからないことが起きてしまいます。こういう場合にはエネルギーが長期間で安定するよう工夫した数値計算を行う必要があります。シンプレクティック性を持つ計算などと呼ばれますが、いちいち書いていたらとんでもない量になるのでここでは割愛します。

私の研究分野でよく使われる速度ベルレ法などの計算手法はうまい具合にこのシンプレクティック性を持っており、2次精度の手法ですが、計算が簡便でしかも陽解法でありながら長期間でエネルギーが変化しにくいという便利な性質を持っています。

系としてはとっても単純でありながら、運動はとても不規則です。時には大きくゆったり回ったり、時には先の振り子だけぐるぐる回ったり、それらを予告なく切り替え続けます。

このような運動はカオスと呼ばれますが、単純な法則から複雑な運動が生まれる例は他にもたくさんあって、例えば前回の反応拡散系も支配している方程式は単純ながら模様が複雑になったりもします。

それなら逆に、世の中の一見複雑な現象も実は単純な法則によって成り立っているのではないか?と考えるものの見方があります。では複雑な運動から単純な法則は見つけることができるのか?振り子ではなく、例えば世の中の気候の変化、株価、人の行動のような複雑な現象は実は単純な法則から成り立っているのではないか?そうだとしたらその法則は何か?という疑問に答えるのが、「複雑系」と呼ばれる学問のひとつの目標なのではないか、と私は勝手に考えています。

まあ、私は複雑系の研究者ではないので、好き勝手に妄想しているところが大きいのですが。

ぞうさんの何でもノート

思いついたことを書いていきます

カテゴリー別アーカイブ: 作ってみた

[作ってみた]物理ぷよぷよ

拍と拍子の概念を宇宙人に伝えよう

「進捗・どう・です・か」をランダムに表示し「進捗どうですか」が完成したら煽ってくるプログラム

概要

生成された煽り文句たち

マンデルブロー集合へ飛び込む

[フラクタル] マンデルブロー集合 10の200乗倍への挑戦

院試勉強しながらついTwitterをしてしまうあなたに「院試まで残り何秒bot」作ってみた

[フラクタル] マンデルブロー集合 10の200乗倍への挑戦

[絶対に儲かる賭け?] サンクトペテルブルグのパラドックスを律儀に試す

ソースコード

結果

運動方程式の数値計算と速度ベルレ法の威力 (後半)

この粒子の振る舞いについて

以前の3重振り子を使わなかった理由

ソースコード

運動方程式の数値計算と速度ベルレ法の威力 (前半)

速度ベルレ法

[作ってみた] [プログラミング] 剛体3重振り子