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

文責：藤原亮（2012年4月10日）カテゴリ：仮想物理実験室(325)｜計算物理学(165)

これまで，以下の振子の運動シミュレーションを行ってきました．
・ラグランジュ運動方程式１：極座標を用いた単振子
・ラグランジュ運動方程式２：極座標を用いた球面振子
・ラグランジュ未定乗数法を用いた球面振子のシミュレーション
・ラグランジュ未定乗数法を用いた２重振子のシミュレーション
これまでのシミュレーションを見てきて，おもりの数が増えたらどんな動きをするのだろうか？という疑問が湧いた人もいると思います．
ここでは，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に対応している必要があります。

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

	【受講生募集】『natural science 科学・技術講座』新講座のご案内 2024.01.09 【大草芳江｜TOPICS】
	サイエンス・デイオブザイヤー2023贈賞式（文部科学大臣賞表彰等）を開催しました 2023.11.07 【大草芳江｜TOPICS】
	■「natural science 科学・技術講座」夏休み短期講座のご案内 2023.07.14 【遠藤理平｜TOPICS】
	「国際イノベーションコンテスト2022世界大会」でnatural science『科学・技術講座』チームが世界3等入賞 2023.07.12 【大草芳江｜TOPICS｜パブリシティ】
	サイエンス・デイ限定『学都「仙台・宮城」サイエンスマップ～光編～』第7版プレゼント申込開始（先着500名）！ 2023.07.08 【遠藤理平｜TOPICS】
	学都「仙台・宮城」サイエンス・デイ2023　7月16日（日）開催！ 2023.07.08 【大草芳江｜TOPICS】
	人材募集のご案内 2022.08.19 【大草芳江｜TOPICS】
	学都「仙台・宮城」サイエンス・デイ2022　7月17日（日）開催！ 2022.07.01 【大草芳江｜TOPICS】
	「第13回国際イノベーションコンテスト2022」国内予選大会　natural science チームが第３位入賞、４年連続通算８回目の世界大会出場へ 2022.05.23 【大草芳江｜TOPICS】
	「キャンパスベンチャーグランプリ2021」東北大会で最優秀賞、全国大会で日刊工業新聞社賞を受賞 2022.03.09 【大草芳江｜TOPICS】
	【プレスリリース】「第12回国際イノベーションコンテスト2021」世界大会　natural science チームが世界３等入賞 2022.01.17 【大草芳江｜TOPICS】
	【第８話】有限の高さの障壁へ照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.29 【遠藤理平｜仮想物理実験室】
	【第７話】無限に高い障壁に向けた電子パルスの照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.24 【遠藤理平｜仮想物理実験室】
	【第６話】無限に高い障壁へ照射アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.23 【遠藤理平｜仮想物理実験室】
	【第５話】電子パルスの運動アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.21 【遠藤理平｜仮想物理実験室】
	【第４話】電子パルスの作り方【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.20 【遠藤理平｜仮想物理実験室】
	【第３話】自由粒子の運動アニメーション【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.19 【遠藤理平｜仮想物理実験室】
	【第２話】自由粒子のスナップショット【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.18 【遠藤理平｜仮想物理実験室】
	【第１話】プログラムの動作確認【Pythonコピペで量子力学完全攻略マニュアル】 2021.09.17 【遠藤理平｜仮想物理実験室】
	「natural science 科学・技術講座」講師３人で設立した株式会社weCAN。アイディアの商品化第一弾「美姿勢メガネ」を朝日新聞さんに取材いただきました。 2021.08.23 【遠藤理平｜パブリシティ】

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

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

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

一階化

境界条件

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

N重振り子シミュレータ

関連記事

仮想物理実験室

計算物理学

Ranking アクセスランキング

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

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

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

一階化

境界条件

張力に含まれる係数を求める

N重振り子シミュレータ

関連記事

仮想物理実験室

計算物理学

Ranking アクセスランキング

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