HOME > natural science Laboratory > コンピュータ・シミュレーション講座 > 仮想物理実験室

N重振り子のアルゴリズムとシミュレータ

文責:藤原 亮 (2012年4月10日) カテゴリ:仮想物理実験室(247)計算物理学(126)

これまで,以下の振子の運動シミュレーションを行ってきました.
ラグランジュ運動方程式1:極座標を用いた単振子
ラグランジュ運動方程式2:極座標を用いた球面振子
ラグランジュ未定乗数法を用いた球面振子のシミュレーション
ラグランジュ未定乗数法を用いた2重振子のシミュレーション
これまでのシミュレーションを見てきて,おもりの数が増えたらどんな動きをするのだろうか? という疑問が湧いた人もいると思います.
ここでは,N重振り子の運動方程式を定式化し,HTML5(canvas要素)+WebGL(Three.jsライブラリ)を用いた, プラグインなしでお手軽に実行できるシミュレータで,その運動を再現します.

3次元空間中にある,質量のない剛体棒でつながれたN個の質点で構成された振り子の運動方程式を求めてみよう!

ここでは,シミュレーションを行うことを前提として,N重振り子の運動方程式を記述していく. また,一般的に振り子の運動を記述するときは,剛体棒の角度	heta_iなどを使うが, この方法は座標の数が減り解析がしやすくなる反面, 質点が剛体棒から受ける拘束力などといった情報がカットされているのだ. したがって今回は,質点のデカルト座標系における位置質点が剛体棒から受ける拘束力 で振り子の運動を記述してみようと思う.

運動方程式を求める前に,振り子を構成する要素に名前を付ける. 振り子が固定されている点を\mathbf{r}_0, 固定点から数えてi番目の質点の位置を\mathbf{r}_i, その質点の質量をm_iとする. また,振り子の質点の数をN,重力加速度を\mathbf{g}=(0 , 0 , -g)^	ranとする. i番目の質点には,重力m_i \mathbf{g}と, 前の質点との間に働く張力- \alpha_{i} (\mathbf{r}_{i} - \mathbf{r}_{i-1}), 後の質点との間に働く張力\alpha_{i+1} (\mathbf{r}_{i+1} - \mathbf{r}_{i})が作用する.


運動方程式と,その解き方の方針

これを運動方程式に書き表すと式(1)のとおりとなる.

(1)

上式における,位置\{\mathbf{r}_i\},速度\{\dot{\mathbf{r}}_i\},張力に含まれる係数\{\alpha_i\}という 7N個の変数を求めれば,振り子を再現できるのだ!! さて,求めたい変数は7N個あるのだが,式(1)は3N本しかない. どうするのか?

まず,「一階化」というテクニックや「拘束条件」を使って式の本数を7N本とし,未知変数と式の数を一致させる. この操作をすると,位置\{\mathbf{r}_i\},速度\{\dot{\mathbf{r}}_i\}をRunge-Kutta法などの差分スキームを使って求めることができる. ただ,\{ \alpha_i \}が決まらなければ,位置\{\mathbf{r}_i\},速度\{\dot{\mathbf{r}}_i\}は決めることができないのだ. ということで,\{ \alpha_i \}をいかに決定するかということが,ここでのヤマ場となる. 結論からいうと,\{ \alpha_i \}に関する連立方程式を立式するのだが, いろいろなベクトル演算のテクニックを存分に活用していくので,頑張ってついてきてくれ!!


一階化

まずは「一階化」を行う. 式(1)中にある加速度\{\ddot{\mathbf{r}}_i\}は求めたい変数ではないので,これを消去したい. ということで式(1)を,\{\mathbf{r}_i\}\{ \mathbf{v}_i \} = \{ \dot{\mathbf{r}}_i \}の一階微分だけで表現する. これが「一階化」なのだ!

(2)

境界条件

しかし,一階化しても式(2)は6N本しかないので,まだ全ての変数を決定することはできない. そこで拘束条件を考えてみる. 各質点のあいだには質点間の距離を一定に保つ拘束条件がある. この条件式は式(3)となる.

(3)

張力に含まれる係数\{\alpha_i\}を求める

さて,これが今回のヤマ場だ!! 運動方程式(1)と拘束条件式(3)から, 位置\{\mathbf{r}_i\}と速度\{ \mathbf{v}_i \} = \{ \dot{\mathbf{r}}_i \}で表現された, 張力に含まれる係数\{ \alpha_i \}に関する連立方程式を立式する. この\{ \alpha_i \}さえ求められれば式(2)をRunge-Kutta法などの差分スキームで存分にブン回せるのだ!

まずは下ごしらえとして,式(3)を時刻tで微分する

(4)

式(4)の両辺を2で割り,さらにもう1回tで微分する.

(5)

相対位置\{ ({\mathbf{r}}_{i} - {\mathbf{r}}_{i-1}) \}の二階微分\{ (\ddot{\mathbf{r}}_{i} - \ddot{\mathbf{r}}_{i-1}) \}が出てきた! 運動方程式(1)からも,うまく\{ (\ddot{\mathbf{r}}_{i} - \ddot{\mathbf{r}}_{i-1}) \}を出して消去すれば,\{\alpha_i\}を求められそうだぞ! ということで,式(1)の両辺を質量\{ m_i \}で割ってみる.

(6)

さてここで\{ (\ddot{\mathbf{r}}_{i} - \ddot{\mathbf{r}}_{i-1}) \}を出すために,i番目の式から(i-1)番目の式を辺々引く.

(7)

ここで\ddot{\mathbf{r}_0}はどうするんだ?と思った人がいるかもしれない. {\mathbf{r}_0}は固定点なので力がつり合っており\ddot{\mathbf{r}_0}=\mathbf{0}になっている,としておこう. ここで式(7)の両辺と\{ ({\mathbf{r}}_{i} - {\mathbf{r}}_{i-1}) \}の内積をとる.

(8)

式(5)から導かれる関係式(\ddot{\mathbf{r}}_{i} - \ddot{\mathbf{r}}_{i-1}) \cdot (\mathbf{r}_{i} - \mathbf{r}_{i-1}) = - (\dot{\mathbf{r}}_{i} - \dot{\mathbf{r}}_{i-1})^2を使って, 式(8)の左辺にある相対位置の二階微分\{ (\ddot{\mathbf{r}}_{i} - \ddot{\mathbf{r}}_{i-1}) \}を消そう. これでやっと加速度を消去できたな!} ついでに,後々のために右辺と左辺を入れ替えたり係数を掛ける順番を少し入れ替えておく.

(9)

さて,これは\{ \alpha_i \}についての連立方程式になりそうだな... とりあえず式(9)を行列化してみよう. 式(9)を下記のとおりに置き換える.

(10)

ただしここで,

(11)

なのだが,\mathbf{R}

となる. あとは,Gauss-Jordanの消去法などを用いて\mathbf{\alpha}を求めれば完了だ! やっと\{ \alpha_i \}を求められたな!!おつかれっす!!

N重振り子シミュレータ

本シミュレータは、HTML5(canvas要素)+WebGL(Three.jsライブラリ)を用いて、プラグインなし実行することができます。 ただし、クライアントサイドのブラウザとグラフィックカードが、HTML5とWebGLに対応している必要があります。

↓デモ画像(「点電荷による電気力線シミュレータ」はこちら



▲このページのトップNPO法人 natural science トップ

関連記事

仮想物理実験室







計算物理学







▲このページのトップNPO法人 natural science トップ