任意多重結合振子運動の計算アルゴリズム（直交座標系）

多重結合振子のモデルと運動方程式

今回は多数の振り子が伸び縮みしない任意の形状で接続された多重振子をシミュレーションする方法を解説するよ。下図は$${i}$$番目のおもりに多数のおもりが接続された模式図だよ。図では３つのおもりが結合されているけれども、個数は何個でも問題ないよ$${i }$$番目のおもりに接続された伸び縮みしないひもの張力を$${\boldsymbol{f}_{ki}}$$と表しているよ。また、i番目のおもりの位置を$${\boldsymbol{r}_i}$$、質量を$${m_i}$$とし、外力（重力など）を$${\boldsymbol{F}^{({\rm ex})}_i}$$とするね。

多数の結合が存在するおもりの模式図

$${i }$$番目のおもりの加速度ベクトルは次ように表すことができるね。

$$
\boldsymbol{a}_i = \frac{ \boldsymbol{F}_i}{m_i} = \frac{1}{m_i}\left[ \sum\limits_{k} \boldsymbol{f}_{ki}+\boldsymbol{F}^{\,({\rm ex})}_i\right]
$$

結合されている振り子の個数がN個の場合、$${i=1}$$から$${i=N}$$まで同様の式が得られることになるね。このN本の運動方程式から張力に関する連立方程式を導いて、数値的に解くことで多重振子の運動をシミュレーションすることができるようになるね。

張力に関する連立方程式の導き方

先の式から導かれるN本の運動方程式はそれぞれ３次元ベクトル量なので3N本の方程式となるね。しかしながら、$${i}$$番目と$${k}$$番目のおもりを結合するひもの張力$${\boldsymbol{f}_{ik}}$$の向きは$${i}$$番目と$${k}$$番目のおもり位置関係で一意に決める事ができるので、実際は大きさのみが未知の値となるね。このことから張力に関する3N本の方程式はN本に集約することができるわけだね。まずはその方法を解説するよ。

番目とj番目のおもりが結合された模式図

上図は$${i }$$番目と$${j}$$番目のおもりが伸び縮みしないひもで結合された模式図です。この２個のおもりに着目して張力に関する関係式を導出するね。

位置ベクトル、速度ベクトル、加速度ベクトルの条件

２つのおもりの距離は一定に固定されているため、各々は自由に運動することができません。２つのおもりが満たす条件式が存在するね。その条件式の出発点は「２つのおもりの距離が一定」となるという

$$
|\boldsymbol{r}_i-\boldsymbol{r}_j|=L_{ij}
$$

だね。この条件式を時間微分することで２つのおもりの運動に関する条件式を導出することができるね。条件式を内積の形

$$
(\boldsymbol{r}_i-\boldsymbol{r}_j)\cdot (\boldsymbol{r}_i-\boldsymbol{r}_j)=L_{ij}^2
$$

と表しておいて両辺を時刻$${t}$$で微分するね。

$$
(\boldsymbol{v}_i-\boldsymbol{v}_j)\cdot (\boldsymbol{r}_i-\boldsymbol{r}_j)=0
$$

この式は速度ベクトルと位置ベクトルが満たすべき関係式で、相対位置ベクトルと相対速度ベクトルは垂直であることを意味しているね。さらに時刻tで微分するね。

$$
(\boldsymbol{a}_i-\boldsymbol{a}_j)\cdot (\boldsymbol{r}_i-\boldsymbol{r}_j) + |\boldsymbol{v}_i-\boldsymbol{v}_j|^2=0
$$

これは加速度ベクトルと速度ベクトル、位置ベクトルが満たすべき関係式だね。このように、元の結合条件を時刻で微分することで運動に関する様々な量の関係式を導くことができるわけだね。

i番目とj番目のおもりの関係の導出

ニュートンの運動方程式に位置、速度、加速度に関する条件を課します。$${i}$$番目と$${j}$$番目のおもりの運動方程式を以下のとおり表しておきます。

$$
\boldsymbol{a}_{i} = \sum\limits_{k} \frac{\boldsymbol{f}_{ki}}{m_i}+\frac{\boldsymbol{F}^{\,({\rm ex})}_i }{m_i}
$$

$$
\boldsymbol{a}_{j} = \sum\limits_{k} \frac{\boldsymbol{f}_{kj}}{m_j}+\frac{\boldsymbol{F}^{\,({\rm ex})}_j }{m_j}
$$

先に導出した加速度が満たすべき条件を課すために両式の片々を引き算して位置ベクトルの差との内積をとるね。

$$
(\boldsymbol{a}_i- \boldsymbol{a}_j)\cdot( \boldsymbol{r}_{i} - \boldsymbol{r}_{j} ) = \left\{\sum\limits_{k} \left[ \frac{ \boldsymbol{f}_{ki}}{m_i} - \frac{\boldsymbol{f}_{kj}}{m_j} \right] + \left[ \frac{\boldsymbol{F}^{\,({\rm ex})}_i}{m_i} -\frac{\boldsymbol{F}^{\,({\rm ex})}_j}{m_j} \right] \right\} \cdot( \boldsymbol{r}_{i} - \boldsymbol{r}_{j} )
$$

この関係式を整理するために、$${i}$$番目と$${j}$$番目のおもりをつなぐ張力$${\boldsymbol{f}_{ij}}$$を、おもりの位置ベクトルを用いて表しておきます。

$$
\boldsymbol{f}_{ij}=f_{ij}\,\hat{\boldsymbol{n}}_{ij} , f_{ji} = f_{ij}
$$

$$
\hat{\boldsymbol{n}}_{ij} \equiv \frac{ \boldsymbol{r}_i- \boldsymbol{r}_j }{|\boldsymbol{r}_i- \boldsymbol{r}_j|}
$$

$${\hat{\boldsymbol{n}}_{ij}}$$はj番目のおもりの位置を基準としたi番目のおもりの方向を表す単位ベクトルだよ。このように定義することで、大きさが同じで向きが正反対となる張力の作用・反作用の法則（$${\hat{\boldsymbol{n}}_{ji}= - \hat{\boldsymbol{n}}_{ij}}$$）が自然な形で表現できるね。さらに、これまでベクトル量であった未知の量をスカラー量$${f_{ij}}$$へと変換することもできるね。これまでの議論を整理すると、先の加速度ベクトルに関する条件式は次のように得られるね。

$$
-|\boldsymbol{v}_i-\boldsymbol{v}_j|^2 = \sum\limits_{k} \left[ \frac{ L_{ij}\hat{\boldsymbol{n}}_{ki} \cdot \hat{\boldsymbol{n}}_{ij} }{m_i} f_{ki} - \frac{L_{ij}\hat{\boldsymbol{n}}_{kj} \cdot \hat{\boldsymbol{n}}_{ij} }{m_j} f_{kj} \right] + \left[ \frac{\boldsymbol{F}^{\,({\rm ex})}_i}{m_i} -\frac{\boldsymbol{F}^{\,({\rm ex})}_j}{m_j} \right] \cdot( \boldsymbol{r}_{i} - \boldsymbol{r}_{j} )
$$

未知の量$${f_{ij}}$$を左辺に、既知の量を右辺に移項した結果が次のとおりです。

【運動方程式】張力（$${f_{ij}}$$）に関する連立方程式

$$
\sum\limits_{k} \left[ \frac{\hat{\boldsymbol{n}}_{ki} \cdot \hat{\boldsymbol{n}}_{ij} }{m_i} f_{ki} - \frac{ \hat{\boldsymbol{n}}_{kj} \cdot \hat{\boldsymbol{n}}_{ij} }{m_j} f_{kj} \right] = \frac{1}{L_{ij}}\left[-|\boldsymbol{v}_i-\boldsymbol{v}_j|^2-\left[ \frac{\boldsymbol{F}^{\,({\rm ex})}_i}{m_i} -\frac{\boldsymbol{F}^{\,({\rm ex})}_j}{m_j} \right] \cdot( \boldsymbol{r}_{i} - \boldsymbol{r}_{j} ) \right]
$$

おもりの位置ベクトル$${\boldsymbol{r}}$$と速度ベクトル$${\boldsymbol{v}}$$、外力$${\boldsymbol{F}^{\,({\rm ex})}}$$、方向ベクトル$${\hat{n}}$$、質量$${m}$$、距離$${L}$$はすべて既知の量で、未知の量は張力の大きさ$${f}$$だね。N個のおもりにそれぞれひもに結合されている場合（結合されていないペアも可）、未知の張力は$${{}_N{\rm C}_2=N(N-1)/2}$$個存在するのに対して、両辺のインデックス$${i}$$と$${j}$$を全ペアに対して与えることで、同数の方程式を与えることができるね。つまり、もとの連立方程式を解くことで全張力を決定することができるわけだね。なお、ひもで結合されていないペアの存在する場合でも、そのペアに対する方向ベクトルを０（$${\hat{\boldsymbol{n}}}$$）と定義することでそのままの表記で問題ないね。

強制振り子運動で検証

先の連立方程式を検証するために、以前に示した支点と振動子が1個の場合の強制振子の計算アルゴリズムを導出してみるね。先の右辺の和は$${k=i}$$と$${j}$$の場合のみとなるね。　

$$
- \left[\frac{ 1 }{m_i}+ \frac{1 }{m_j} \right]f_{ij} = \frac{1}{L_{ij}}\left[-|\boldsymbol{v}_i-\boldsymbol{v}_j|^2-\left[ \frac{\boldsymbol{F}^{\,({\rm ex})}_i}{m_i} -\frac{\boldsymbol{F}^{\,({\rm ex})}_j}{m_j} \right] \cdot( \boldsymbol{r}_{i} - \boldsymbol{r}_{j} )\right]
$$

この式はすでに未知の量は$${f_{ij}}$$ひとつなので実質的には連立方程式にはなってないね。$${f_{ij}}$$について解くことで張力が得られるね。

強制振子運動を想定して$${i}$$番目を振り子のおもり、$${j}$$番目を支点と考えるね。強制振子は支点が外場とは関係なく移動するので、$${m_j=\infty}$$としつつ速度を$${\boldsymbol{v}_j =\boldsymbol{v}^{\,({\rm ex})}_j }$$、加速度を$${\boldsymbol{a}^{\,({\rm ex})}_j = \boldsymbol{F}^{\,({\rm ex})}_j/m_j }$$と考えるね。おもりの$${\boldsymbol{a}_i = \boldsymbol{F}^{\,({\rm ex})}_i/m_i = \boldsymbol{g}}$$とを与えて整理すると

$$
f_{ij}=\frac{m_i}{L_{ij}}\left[ |\boldsymbol{v}_i -\boldsymbol{v}^{\,({\rm ex})}_j|^2 + (\boldsymbol{g}-\boldsymbol{a}^{\,({\rm ex})}_j)\cdot( \boldsymbol{r}_{i} - \boldsymbol{r}^{\,({\rm ex})}_{j} ) \right]
$$

となるね。これが張力$${S}$$で、以前に導出した結果と一致するね。強制振子運動における張力を正しく導出できることが示されました。

次回はこのアルゴリズムを用いた最初の例として、おなじみの２重振子運動シミュレーションを行うよ。