カテゴリー別アーカイブ: お勉強

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

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

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

一番丁寧な運動方程式の解き方は解析的に解くことですが、ほとんどの場合うまい方法がないことがわかっています。例えば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次精度の手法ですが、計算が簡便でしかも陽解法でありながら長期間でエネルギーが変化しにくいという便利な性質を持っています。

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

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

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

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

確率の問題再考「ある家庭の子供2人のうち，1人は男であることがわかっている．このとき…」

コメントをどうぞ

「ある家庭の子供2人のうち，1人は男であることがわかっている．このとき，もう1人が男である確率は?」

5725091674_9d1f4a5688_o

確率の分野でとても有名な問題のひとつの表現です．確率の入門書に限らず，この手のパズルを紹介する書物などでもよく見かけ，またインターネット上の様々なコミュニティで定期的に話題になります．
Yahoo!知恵袋でも検索をかけると数十件程度同じ内容の質問が並んでいたりします．テレビ番組で取り上げられたこともあるようですね．

なんやかんやと炎上したあげく，答えは1/3だという結論に達するのが定番の流れです．

私はこの問題，いつも間違ってるんじゃないか?と思っています．少なくとも問題文に矛盾があって，答えを1/3とするのは乱暴すぎるような気がしています．

まず挙げられる(間違いとされる)答えは

「性別なんてもう片方の性別には関係していない．だから1/2だ」

というものです．これに対する反論はだいたい以下のようになります．

「想定されるパターンは(男，男)，(男，女)，(女，男)，(女，女)の4種類である．このうち4番目は問題文より除外されるので，残り3パターンから考えて1/3だ」

この説明は一見正しそうに見えます．ただ，上の説明で問題文の「もう1人が男」を考えるとどうにも不穏な空気が漂い始めます．特に(男，男)のときのもう1人ってどっちでしょう?

問題文を見てみると，問題文前半では単に「1人は～である」という表現が使われているだけです．これは(答えを1/3とする文脈では)「～という人が1人以上いる」と言っているだけであり，特定の人のことは何も記述していません．じゃあその後に続く「もう1人」は?と考えると，そんなものは定義できないことがわかります．
(問題文のつじつまをあわせるなら，「1人は男である」を特定の人にすることになり，答えは1/2になるでしょう．この場合は上の誤答とされていた表現もすんなりくると思います．)

この記述のおかしさは男女というどちらも起こりそうな現象で考えると気づきにくいですが，もっと確率の偏った出来事に置き換えて書けば雰囲気が変わります．

「太陽と月のうち少なくとも片方は東から出て西に沈む．もう片方は～」

と書くと，いやいや，もう片方なんてないでしょう，となると思います．(前半の文には特に間違いはありません．)

出題者の意向通り答えを1/3としたいなら，問題文は「両方男である確率は?」としないとおかしくなります．細かい違いのようですが意味が全然変わってきます．

以前Slashdotに載った記事で「子供が二人いて、(少なくとも) 一人は火曜日生まれの男の子である。もう一人の子供も男の子である確率は ?」という問題を扱うものがありました．

問題の原文は “I have two children, one of whom is a son born on a Tuesday. What is the probability that I have two boys?” となっており，問題はなさそうです．
ところが英語版の記事に転載された段階で “What’s the probability that my other child is a boy?“となってしまい，ここから変な話が始まってしまっています．

(ちなみに原文での答えは13/27になります．これは火曜日という一見どうでもよさそうな条件が答えに影響するという点に面白さがある，と紹介されています．考えてみて下さい．)

まあこんなのは言葉遊びだと言われればそうかもしれませんが，一度問題を正しく把握してしまえばあとは適当な手段で計算するだけなのであって，その段階に落としこむまでのほうが大事であることはよくあることではないかと思います．

私が一番気になったのは，こういう質問に答える側の人があまりこういう議論(問題文に矛盾がある，とか)をしないのだなと思ったことです．元が有名な問題なのでたぶんどこかで聞いた/議論した答えを輸入しているのだと思いますが，もし正しい答えとされるものを十分な吟味なしに「正しい答え」として広めているのであれば，それは自分の出した答えが間違っている以上に危険なことではないかと私は思います．

後輩向け: 卒論の研究室の選び方まとめ

1件のフィードバック

どの本をとるべきか?

ちょっとした事情でwebサーバーを止めていましたが，復活させました．

私の大学の工学部機械系2学科は4年生の4月に卒論の研究室配属が行われます．研究室からすれば毎年の行事ですが，当の4年生にしてみれば大事なイベントです．
4年生は例年この時期になるとどの研究室がいいかとか，この研究室は誰が志望しているだろうかとか，そういう話があちこちで聞かれるようになります．

一方で，いったい何を基準に選べばいいのか4年生にはよくわからないところもあります．これは研究室が遠くて雰囲気が伝わりにくい(らしい)ことと，そもそも研究生活をしたことがない以上何を重視すればいいのか知りようがないというのがあると思います．

もちろん研究テーマや研究室の専門は大事ですが，機械系は研究室の数がとても多く，また興味のある内容が広い人もいるため，なかなか決めきれない人も多いと思います．
かなりの人数が修士課程まで進むことを考えると，研究室選びというのは大学時代の後半3年間を決める大事なイベントでもあるわけです．

というわけで，新4年生向けに研究室選びに役に立つかもしれない情報をいくつか集めてみました．

念の為に言っておきますが，ここの内容は「何がある研究室は悪い」とか「こういう性格の研究室が良い」とかそういうことを書くつもりではありません．(なるべくニュートラルな立場を取ろうと努めています．)
選ぶときの目安としてどんな違いがあるのかの例を見せられたらと思っています．

まず身近なところで，私の指導教員がいくつか指標を用意してくれています．

４年生が修士で同じ研究室に進学しているか?
修士の定員が少ない場合や，研究室が外部に有名で外部受験生が多い場合など例外もあります．また機械情報は(というか情報理工学系は?)院試の出願時に修士の指導教員まで決めてしまい時間がないため，基本的に卒論配属とおおよそ同じ志望になっていると思います．
OBが研究室に出入りしているか? OBとの交流が活発か?
OB以外にも，他の研究室の先生やメンバーなどとの交流も聞いてみましょう．つながりをもっている研究室もあります．
教授との相性
人間同士なので，あまり他人の意見に流され過ぎないように．

next49さんのブログの「大学院修士向け研究室情報チェックリスト」にも，研究室選びについて様々なことが述べられています．修士向けとありますが卒論生向けでも問題ありません．この方は大学の教員で，他にも数多くの卒論生向けのエントリーや研究全般についてのエントリーがありますので参考になるかと思います．

いくつか抜粋してみます．

実験系か理論系か。拘束時間や研究のスタイルが大きくことなる
実験系は分野にもよりますが，機械の都合や故障などで一時的に時間が制限される場合があります．
両方扱う研究室もありますが，片方だけの場合はもう片方とどういう関係を持ってるのか調べてみてもいいかもしれません．
そもそも研究室のメンバーは何人か？
人数が多いと研究室内でいくつかの研究グループに分割されていたりします．指導教員との距離感に影響します．もちろん人によって丁度良い距離があると思います．
学生が多ければ教員の人数も増えますが，どちらかに比率が偏っている場合は，どういう理由があるのか考えてみてもいいでしょう．
管理か放任か
機械系でコアタイムのある研究室はそれほどないと思いますが，週次イベントの数や時間はまちまちです．
たいてい週のどこかで研究会があります．それ以外にもグループミーティングがある場合なども多いです．他にも例えば機械情報(正確には知能機械情報学専攻)の研究室は必修科目として研究室輪講があります．
メンバーの多様性はどうか？(出身校、学部は、属していたサークル、人種、出身国、地域など）
留学生の人数や結びつきの強い国などにはけっこうばらつきがあります．
また，修士課程などで複数の枠を持っている研究室があります．機械系に近い分野だと創造情報学専攻，学際情報学府，バイオエンジニアリング専攻などがあります．研究室の多様性もそうですが，院試の入り口が複数になる場合もあります．

他にも，せっかくなのでいくつか見ておくとおもしろいかなと思うものを挙げておきます．

卒論生の面倒を見る人
教授や准教授であったり，助教であったり，あるいはドクターであったり，比率は研究室で様々です．
(研究テーマ設定などで)どれくらいの自由度があるのか?
例えば研究室として大きなプロジェクトに関わっていたりすると，4年生もそれに合流して研究の一部を担当したりすることがよくあります．これもどちらのスタイルがよいかは個人差があります．
テーマと問題点が決まっている方がもちろん研究は進みやすく，おそらく知識の吸収もしやすくなります．一方でテーマと問題点の設定というのは研究においてとても重要なプロセスで(一番大事という人もいる)，避けては通れません．一方でそれには深い知識が必要で…と，鶏と卵みたいな関係になっている気がします．
研究室の入りやすさ
「部外者は基本的に禁止」から比較的開放的な研究室まで様々です．研究室の方針が垣間見えます．
同期のメンバー
後輩の好きな先輩というのも多分たくさんいますが，なんだかんだで困ったときには同期が一番頼りになります．仲のいい人がいるからというのも，必ずしも理由にならないわけではないと思います．
これもきっと性格によりますね．

他にも違いを考えだすときりがありませんが，入ってみないとわからないこともたくさんあるので，あんまり深く悩みすぎないようにしましょう．研究室の形態も毎年どんどん変わっていきます．
先生方は皆とんでもなく優秀で，やりたいと思えばいくらでも力になってもらえるはずです．

研究室選びと4月からの研究生活，肩の力を入れすぎずにがんばってください．

男子トイレの力学

コメントをどうぞ

知性は、方法や道具に対しては鋭い鑑識眼を持っていますが、目的や価値については盲目です。 – アインシュタイン

大学院の夏学期の授業で自由課題が出たのですが，その提出物を授業で発表することになりました．

課題で愉快な結果が出たわけでもないのですが，せっかくなのでそのとき作ったプログラムとその結果を書き留めておきます．
テーマは「モンテカルロ法」．何か身近な問題についてシミュレーションを行ってこいとのことでした．

モンテカルロ法は乱択アルゴリズムの一種です．ざっくり言ってしまえばたくさんやってみて統計をとるという方法ですね．
コンピュータのくせに「答えは厳密ではないが，たぶんだいたい合ってる」という結果の出方が大好きです．

さて，たくさんのものをぼんやり見る…となるとやはり気になるのは人の群れのシミュレーションです．人って群れると不思議な振る舞いをしますよね．
研究となると面倒な分野ですが，授業の課題なら気にせず変なことができます．

それで扱う現象ですが，人って誰かの隣には座らないんですよね．電車などに乗っているとよく見る景色です．
これは映画館とか大学の講義でも同じで，隣が開いてる席があったらそっちを使うというのが人間というものの性質のようです．

男子トイレでも同じことが起きますね．大きなトイレでがら空きなら誰かの隣の便器というのは普通使いません．もしがらがらのトイレで隣に人がやってきたら，何か他の目的があるのではという考えが頭をよぎります．

で，あのふるまいのシミュレーションをしてみようと思ったわけです．

トイレのモデルは以下のようになります．

便器の数をn個として，そこに人が次々と到着します．各個人の便器使用時間はガウス分布，到着時間はポアソン到着としています．

各個人は両隣が空いている便器があればそれを，なければそれ以外の便器を使います．全部埋まっていたら待ち行列に並びます．待ち行列の人は便器が空き次第詰めていきます．
同程度の魅力の便器があればランダムに選びます．

パラメータは便器数n, 便器使用時間の平均μおよび標準偏差σ，それと平均到着時間間隔λとなります．
便器使用時間についてはどこかにデータがあればよかったのですが，見つからなかったので自分で1週間計測しました．
平均は20秒，標準偏差は3秒となりました．今回は両側に3秒追加して26±3秒としてシミュレーションを行いました．

n=5として平均到着時間間隔λを変化させた結果が以下になります．右の番号はトイレの便器ですね．端から順に1から5です．

横軸は到着頻度（λ）になります．左端は人が来ない状況で，右に行くにつれてトイレが賑わっていきます．右端は常にトイレが使われているという相当やばい状態ですね．現実にこんなトイレがあったら早いところ増設したほうが事故を未然に防げていいと思います．

人がとてもまばらな時は，どの便器も同様に使われます．また相当混んでいる時も同様に使われます．
まあ他人がいなければ区別なく使いますし，混んでる時は選ぶ余地もなく全部使われてるので当然とも言えます．

問題はこの中間のときで，ちょうど限界の半分くらいの人が押し寄せるときに便器の使用頻度に差が生じます．
端のトイレが比較的よく使われる一方で，端から2番目のトイレ（2,4）の使用頻度は落ちています．おそらく端っこは隣が片方しかないため，そこそこ混んでいても入りやすいのだと思います．それでその隣は使われないのでしょうね．

以下が便器の数n=31とした結果です．大きなトイレですね．

多くの便器達はおおよそグラフの中央に集まっています．これは便器ごとに使用率が変わらないことを意味しています．
一方で上に大きく曲がった線が2本あります．これが両端の便器です．使用率は他のものに対して有意に高くなっています．
下の2本は端から2番目の便器になります．こちらは使われにくくなっていますね．

現実問題を表現できていますかね?
今度駅とかのトイレに行ったときは便器の使用状況に気をつけてみるのもいいかもしれません．

この結果から何がわかるか?というと，まあせいぜい端の便器は汚れてるかもしれないよ，ということくらいです．
まあトイレに限らずなんでもそうかもしれません．大学の教室の椅子でも5席あったら1,5番目，詰まったらその後3番目から埋まりますね．

これ自体は非常に単純なシミュレーションですが，簡単な仮定をおくだけでこのようなデータが得られるのがシミュレーションの面白いところです．
気軽にコードが書けるとこういった遊びの幅も広がります．

コメントも書いてないものですが一応ソースコードを置いておきます．Visual Studio 2010で作ったものですが単純なC++のコードなのでだいたいの環境で動くと思います．

restroom

発表に使った資料も置いておきます．後半は違う話になっていますがいずれここで書けたら面白いかなと思っています．

マルチスケール課題1

シミュレーションの内容に補足をすると，今回はキネティックモンテカルロ法という手法になっています．これは扱う時間のスケールに対して稀な現象を含めて扱うときに使われる手法のひとつで，系で起きるイベントをリストアップし，次のイベントまでにかかる時間を計算してイベント単位で時間を進めていく方法です．

また，このコードにはトイレの待ち行列が伸びすぎるとクラッシュする機能があります．シミュレーションとはいえあんまり我慢させるのはかわいそうですから．

ちなみにトイレシミュレータの原案は昨年の研究室の先輩方でした．先輩方ありがとうございました．

[昔話] 振り子のシミュレーションをした話

1件のフィードバック

今回は真面目系の話．駒場時代の昔話です．

[追記] プログラムを作り直しました。 [作ってみた] [プログラミング] 剛体3重振り子

シミュレーションの課題をやる元気が出ないので昔話を書いてモチベーションを上げようという魂胆です．

東大の理系1・2年には物理実験という授業があります．毎週1回何かを実験して，その結果をまとめて先生に報告したら終了というものですね．

報告はたいてい他の人と時間が重なるので，待機用に控え室というのが用意されています．
で，この部屋，いったい誰が作ったのかわかりませんが物理っぽいおもちゃがたくさん取り揃えてあります．でっかいジャイロと回転台とか，ニュートンのゆりかごとか，そういう類のものです．

ニュートンのゆりかご(wikipedia)

その中でひときわごつくて目を引いたのが，金属製の振り子でした．もっと丁寧にいうと「剛体3重振り子」というもので，回転する棒の先に棒が，またその先に棒がついた形をした細長い振り子です．

これの遊び方は?
先をつまみ上げて落とすだけです．

じゃあ何がおもしろいのか，というと，振り子の挙動です．ふつう振り子というのは左右にゆっさゆっさと揺れるものですが，この振り子は訳の分からない運動を始めます．全体でぐるんぐるん回ったかと思えば突然暴れだしたり，かと思えば先の棒だけが超高速で回りだしたりします．

以下は剛体2重振り子ですが，雰囲気は伝わると思います．

「連成振り子の運動」というのは大学数学では定番の練習問題で，たいていは複雑に見える運動も実は計算してみると単振動の組み合わせでしたね，という結論が待っています．もう大学の授業でかれこれ6回は見た気がします．

そんな中奇妙な振る舞いをする振り子が登場するわけです．中に奇妙なカラクリがあるとかならわかりますが，眼の前にあるのはただ棒が3本つながっただけの代物です．複雑な現象も単純に表されるんだよと言われていた矢先にこれです．

納得していたものが頭の中から吹き飛ぶ感覚はある意味快感ですらありました．

さて，実験の報告も終わって帰ってきて，まあこのわけのわからないものを理解してやろうじゃないか，理論的に記述してやろうじゃないか，となるわけです．

まあそうやって書けばなんかモチベーション高い感じですが，別に目標があったわけでもなく，遊びみたいなものです．

ただし解析的に解けないのはわかっていたため，そこは数値計算で，となります．コンピュータの出番です．
丁度この頃力学はLagrange形式までやっていましたし，プログラミングも一応授業で最低限には習っていました．

となると，残りは微分方程式の数値計算法です．数値計算の分野にはこの頃はじめて出会ったことになります．
数値計算というのはおもしろい分野で，微分とか積分とかの抽象的な数学の問題を，単純な計算作業に落としこむ方法を与えてくれます．

それだけといえばそれだけですが，運動方程式，つまりこの世界そのものを足し算とかけ算だけの世界に引きずり込むという芸当ができるわけです．神様お試し版です．昔そんな映画もありましたね．

で，あれこれ式を書いていくと，最後に”あの”わけのわからない動きが手元ですっかり再現できるようになるわけです．

お手製のプログラムで自分でもよくわからなかった運動を再現するというのはなかなかおもしろいものです．自分で世界を作る感覚，シムシティと同じです．

当時作ったプログラムが公開できたらよかったのですが，ちょっと今手元にないためまた後で用意したいと思います．
かわりにyoutubeにあったものをはっておきます．

2重，3重振り子はカオスの例として有名なんですね．

こんな日を経験してから，計算力学という分野にある程度関心をもつようになります．なんか新しい物性とか流れとかちょっと複雑なことを考えたいときは理論式だけではなく計算機がボトルネックになるのではないか，なんてことをぼんやり考えていたわけです．

私が機械情報工学科に進学することを決めた背景には，このような考えも影響していたのではないかなーと思います．

さて，この遊びで作ったプログラムは特に使い道もないしそのままお蔵入りするかと思っていたのですが，ところがどっこい，学部に進学してからまた引っ張りだしてきて新しい要素が組み込まれることになります．

その話はまたいつかできたらと思います．

ぞうさんの何でもノート

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