全体は部分の総和に勝る。 – アリストテレス
マンデルブロー集合、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円、となるので…
ここまで見てだまされないぞと思った人は、しっかり悩んでみてください。
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回あたりの獲得金額」は減っていくので右下がりになるという具合です。その代わり「もっとすごい大当たり」をひく確率も出てきます。
モンテカルロ法などでは試行回数を増やすと値が落ち着くのですが、この問題はおもしろいことにいくら試行回数を増やしてもばらつきが収まらないようになっています。よくできています。
儲け額の分散も出してみようと一度書きかけたのですが、ちょっと考えてからこの問題に対して無駄な試みであることがわかりました。
コンピュータはしばしば想像力を働かせる道具として役に立つのですが、今回のように神通力の効かないヤバそうな問題ではちょっとわけが違ってきます。こんな問題でも扱えてしまう数学の力は偉大です。
結果の意味とこの賭けの本当の期待値についてはいろいろ考えてみて下さい。
前半の補足です。前半で紹介したシミュレータ云々について。
「なんで星の運動を模してるのに楕円軌道じゃないんだ?」と思った方はとても勘のいい方です。
これは2次元の運動をどう考えるかによるものです。
3次元空間の平面上に粒子を配置させて運動させた場合、重力は距離の2乗に反比例しますが、
3次元空間に無限長の棒を(平行に)置いて運動させた場合、重力は距離の1乗に反比例します。
これはまあ真面目に積分してもそうなりますし、電磁気学を勉強したことのある方は雰囲気で理解できるかと思います。後者は宇宙空間では仮定が不自然すぎるだろうと思われるかもしれませんが、万有引力が3次元空間で逆2乗の法則に従っているとするなら、2次元では逆1乗になるのが自然だとも考えられます。
まあ、今回はそこまで哲学的な話ではなく、単に逆1乗のほうが見栄えが良かったからです。
どう見栄えに影響するかですが、これは無限遠点での挙動の違いを考えればだいたいわかります。位置エネルギーは力の積分(の-1倍)で与えられますが、これより力が逆2乗、逆1乗の場合のそれぞれについて位置エネルギーは逆1乗、logとなります。つまり前者ならエネルギーが無限遠である値に収束するのに対して、逆1乗なら際限なく増えていきます。
つまり、逆2乗ではひとたび粒子が第三宇宙速度に達してしまうと二度と戻ってこないのに対し、逆1乗なら少なくともいずれは引き戻されてくるという大きな違いがあります。見た目に粒子がどこかに飛んでいってしまうのは寂しいですね。
よく「この世が2次元だったらいいのに」という言葉を聞きますが、そんなことをしたら、上記の理由から世界中のものがある一点に集まってしまってつまらない世界になるだろうなあと思っています。
なお、粒子があまり近づきすぎると計算上の問題からヤバいことが起きるので、引力を打ち消すだけの反発力を与えてあります。
現実だと衝突しますが、無理にエネルギーを保存させようとして星が完全弾性衝突をするのもあまりリアリティがないですね。
これは速度ベルレ法の適用範囲に限界があるからです。(下は剛体3重振り子の動画)
シンプレクティック法は、影のハミルトニアンと呼ばれる量を厳密に保存する運動を記述します。この際、位置
しかし都合のいいことに、単純な場合ハミルトニアン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;
}
速度ベルレ法だけ、位置と速度を明示的に扱っているのがわかります。
前々回の振り子の欄でちょっと触れた運動方程式の数値計算の話の続きです。
運動方程式は立てるだけでなく、実際に方程式を解かないとシミュレーションは行えません。
一番丁寧な運動方程式の解き方は解析的に解くことですが、ほとんどの場合うまい方法がないことがわかっています。例えばN体に万有引力が働いて運動する系ではNが3以上になると解析的に解けなくなることが知られています。太陽と地球と月でもうオーバーです。
もっと一般的に扱いたいときは、コンピュータが強力な武器になります。微分方程式を差分方程式に近似してしまえば計算機で簡単に扱えるので、後は時間の刻み幅を十分に小さくとってあげれば微分方程式に近い解が得られるというわけです。
じゃあそれで話は終わりかというとそうでもなく、この差分方程式にする段階がなかなかうまくいかないわけです。それに伴っていろいろな手法が提案されています。今回はその話です。
今回作ったものは互いに引力を及ぼし合う粒子たちのシミュレーションです。
イメージしているのは銀河系的なアレですが、まあ数値積分の例で作っただけなので本気で似せているわけではありません。私は中学校で習う電子の描像のイメージに似ているなと思いました。
IAEAのマークにもあるアレです。
中央に質量の大きめな粒子をひとつ、あとは小さな粒子をばら撒いています。小さな粒子同士も相互作用します。
大事なのはこの系では力学的エネルギー・運動量・角運動量が保存するということです。シミュレーションでエネルギーが保存していなかったら、それは何かおかしいことが起こっているぞ、というわけです。
可視化の部分は振り子の流用ですね。数字キーでシミュレーション手法を変えたり、初期化したり軌跡のオンオフを行ったりします。
後述するオイラー法、ホイン法、ルンゲ・クッタ法、速度ベルレ法を実装しています。
いつものようにプログラム: particle
ソースコード: GitHub
微分方程式を解く際、おそらく最も単純な差分化はオイラー法です。
ニュートンの運動方程式は2階の微分方程式ですが、速度v自体も変数と置いてしまうことで1階の微分方程式に置き換えることができます。以下の式で左辺をx、右辺をfと置いてしまえばいいわけですね。
オイラー法では以下のような差分化を行います。⊿tは十分に短い時間をとります。
見た目にわかりやすいため学習用として頻繁に登場しますが、正確なシミュレーションをしたい場合向きではありません。誤差がとても大きいからです。
誤差が大きい理由は2点の補間を直線、つまり1次式近似にしているからで、それなら2次や3次の項を追加して補間しようというアイデアがあります。実際これはうまく行きます。
よく使われるものに4次多項式で補間するルンゲ・クッタ法というのがあります。
tからt+⊿tまでの4点で傾きを計算して、適当に混ぜることで4次の精度を出しています。式中の諸係数を求めるのにはそこそこの手間がかかります。
次数を上げるほど多項式補間の精度は上がりますが、そうすると計算の安定性の議論から時間刻み⊿tをかなり小さくとらないといけなくなってしまうこと、また計算コストが上がるのでそれなら次数を上げずに時間刻みを減らしたほうが都合がいい場合が多いことなどがあり、結果としてこの4次のルンゲ・クッタ法はよく使われます。
1次、2次、4次の各手法で時間経過とエネルギーのずれをプロットしたものが以下になります。
1次のオイラー法は目に見えて粒子が遠ざかっていってしまいます。2次(ホイン法)と4次(ルンゲ・クッタ法)はシミュレーション結果を見た雰囲気では似ているものの、エネルギーを見るとホイン法はだんだんとずれていくことがわかります。
もちろんオイラー法にも利点はあって、軽くて実装が楽なので、正確さにあまりこだわりがない場合にはいちばん手軽な解法になります。
陰的解法と言われる手法群もあるのですが、一般に複雑なので今回は省略します。
ここで、今回メインで書きたかった手法「速度ベルレ法」を紹介したいと思います。2次の精度をもつ手法です。
2次ならルンゲ・クッタ法より精度が悪いのかというとそうでもありません。先のグラフの倍率を上げたものに速度ベルレ法を追加してみます。
ルンゲ・クッタ法に対してぶれは大きいですが、ホイン法のようにエネルギーが変化していくことはありません。しかしまだこれではルンゲ・クッタ法のほうが頼りがいがあるように見えます。
面白いのは長時間の振る舞いです。このまま横軸のスケールを100倍にしてみます。
ルンゲ・クッタ法もこの長時間では計算誤差が積もってエネルギーがずれていくのに対して、速度ベルレ法はほぼ初期の振る舞いを保ったままとなっています。短期間のぶれは大きいのに、長期間でずれていかない不思議なグラフになっています。
速度ベルレ法の式は以下のようになります。位置と速度を分けて扱うのが今までの手法と違います。
なぜこのような振る舞いを見せるのかにはわけがあって、この速度ベルレ法はシンプレクティック性というとても重要な性質を持つように作られています。
ここからちょっと難しい話ですが、ニュートンの運動方程式は時代を追ってラグランジュ形式、ハミルトン形式へと再公式化されていきます。このハミルトン形式で表現したとき、運動は位相空間上で体積を保つ性質があることがわかっています。
速度ベルレ法も含まれるシンプレクティック積分法と呼ばれる手法では、この位相空間での体積を保存する性質を有していて、結果的に長時間でのシミュレーションで(ある種の)誤差が積もっていく現象が起きなくなっています。
ハミルトン形式におけるとある保存量を厳密に保存しており、時間反転対称性(終了状態から反対向きにシミュレーションを行うと開始状態に行くことができる)、エネルギーなどの保存量を保存する、といった今までの手法にはなかった性質があります。
速度ベルレ法は、私が研究でよく扱う分子動力学法でもよく使われています。系を長時間観察するときなどにこの性質がとても有効だからです。また分子動力学法では計算コストのネックが力の計算になっており、1回に多段で力の計算を行うルンゲ・クッタ法などは相性があまりよくありません。
式自体は単純なため「分子動力学法で使われる数値計算は単純」という表現をよく耳にします。しかしその構成はよく考えられていて、裏側には先人の知恵がたっぷり詰まっていることは確かです。
運動方程式ひとつをとっても様々な数値計算法がありますが、どの積分法が正しいかといった議論はあまり有益ではなく、私のイメージでは「どれも間違っているが、間違え方が違うので場合に応じてどうしても譲れないものに応じて選ぶのが良い」としか言えないかなと思っています。
例えばゲームの物理演算などでは、そこらへんの物が綺麗な放物線を描かなくてもそんなに気になりませんが、プルプル震えたり突然とんでもない勢いで吹っ飛んでいったりすると見た目からしてヤバい、ということになるわけです。それなら計算精度より安定性と計算コストを重視しようとなるわけです。
では、現実のこの世界は、どんな手法で運動方程式が解かれているのでしょう。これはちょっと危険な問いで、そもそも現実世界が運動方程式に厳密に従っているのか、宇宙は連続体なのか、シミュレーション可能なのか、といったことからあれこれ考えなければいけません。
まあでも、私が宇宙をひとつ作っていいと言われたら、たぶんオイラー法を採用します。一番単純だし、世界がそんなに複雑な式でできているとは考えたくないからです。
後半では、今回のシミュレーションで作った系について(なぜ前回の振り子のままではだめだったのか、星の運動としてこの動きの胡散臭さは何か、ソースコード云々)の説明をしたいと思います。
速度ベルレ法と言っても万能ではないのです。
