1. はじめに

1
2018 年 11 月 12 日

多重剛体振り子の方程式
新1

潟 工科 大学 基 礎教 育・教 養系 竹野茂治

はじめに以前、ある問題を考察するために、「n 重剛体振り子」の方程式を考えたことがあった。また最近 その問題を思い出して、ついでにその数値計算もしてみたのだが、あまりうまくいかなかった。もしや方程式の導出が間 違えているのかと思いネット上を探してみたのだが、「n 重 振り子」の方程式はあるものの、「n重剛体振り子」の方程式はまともなものがあまりないようだった。
一方で、「n 重 剛体振り子」のシミュレーション例はネット上にもあるし、カオスとの関 係 で 2 重や 3 重の剛体振り子の実験 例などは結 構目にする。となれば、私が導出した n 重剛体振り子の方程式を置いておくことにも意味がないわ
けではなさそうなので、それをここにまとめておくことにする。

2

設定

通常の振り子は、天井からぶらさがる軽 い糸の先に重りをひとつつけるが、「剛体振り子」は、その糸にあたる部分が剛体、すなわち金 属や木など、形 が変形せずに質 量があるものになっていて、この場合は先端の重りは必要なく、物理振り子と呼 ばれることもあるようである。
剛体振り子の形 状は、必ずしも棒状のものとは限らないが、本稿ではまっすぐな棒状のもののみを考える (1)。
また、多重振り子とは、振り子の先にまた振り子がつく、という構造になっているもので、n 重は n 個の振り子が端同士でつないであるものを指す (細かくいえば、本稿では「直列」な n 重振り子を考える)。
なお、通常の振り子は、重りの先に次 の振り子の糸を結 べば済むが、剛体振り子の場合は、剛体同士が摩擦なく回 転できるような (質 量のない) ジョイントでつなぐ必要がある。また、大きな回 転角 も想定する場合は、その回 転によって剛体同士、あるいは剛体が天井にぶつかることはない (すりぬけ条件 ) とか、何 回 転でも回 転可 能であるといった設定でシミュレーションをしていることも多い。本稿でも、特に回 転角 の制限やぶつかるような設定はせずに考察する。
各 剛体同士は、そのジョイントを通して互いに力を与えるが、それは作用反作用の関 係 になる。

3. 記号

2

図 1: 多重振り子 (左は通常の振り子、右は剛体振り子)

また本稿では、剛体の運動は 1 平面内で起 きるとし、よって回 転もその平面内での回 転とする。n 個の剛体は、同じ材質 、同じ太さのまっすぐな棒状のものとするが、長さは揃っている必要はないとする。
このような方程式の導出には、ラグランジュの運動方程式が良く用いられるようであるが、本稿では、より煩雑ではあるが古典的な方法、すなわちニュートンの運動方程式を用いることにする。

3

記号

1 平面内の剛体の運動は、各 剛体の向きを示す n 個の角 を決 めればすべての状態が決 定するので、その n 個の角 が求 めるべき未知関 数となる。

剛体の棒には、上から順 に A1 , A2 , . . . , An と名前をつけ、各 Aj の長さ、質 量をそれぞれ lj , mj とするが、剛体棒の線密度を ρ とすれば mj = ρlj となる。鉛直下 向きの方向に対して、Aj が右側に回 転した角 を θj とし、Aj の中心 (= 重心) の位置ベクトルを rj , Aj と Aj−1 のつなぎ目の位置ベクトルを r̂j−1 とする(図 2)。
鉛 直下 向きから右に θ 回 転させた単位ベクトルを p(θ)、それをさらに π/2 回 転させたベクトルを q(θ) とするとp(θ) = (sin θ, − cos θ),q(θ) = (cos θ, sin θ)

であり、A1 は原 点からぶらさげるとすれば、 Aj は p(θj ) の方向を向くので、r̂ 0 = 0,r̂j =jX

k=1

lk p(θk )(1)

4. 運動方程式
PSfrag replacements

3r1
A1

θ1r̂ 1

l1r2
θ2r̂ 2l2r3A2

θ3l3r̂ 3A3

図 2: lj , r j , r̂ j , θjとなり、よって rj は
r j = r̂ j−1 +j−1
X
lj
ljp(θj ) =lk p(θk ) + p(θj )22k=1(2)となる。これで、全体が角 θ1 ∼ θn で表されることになる。

4運動方程式

運動方程式は、各 剛体毎にその重心を中心とする方程式を考えればよく、そのために各 剛体にかかる力を考える。
剛体 Aj の重心 rj には、Aj の質 量に対する重力 −mj gey (ey = (0, 1), g は重力加速度) がかかり、つなぎ目 r̂j−1 , r̂j の場所では、 Aj−1 , Aj+1 からの力がかかる。
Aj+1 が r̂ j の場所で Aj に及 ぼす力を T j とすると、その反作用により、Aj はr̂ j の場所で Aj+1 に −T j の力を及 ぼすことになる。
よって、Aj の運動方程式は、
mj r̈ j = T j − T j−1 − mj gey(1 ≤ j ≤ n)(3)となる。ここで、「˙」は時 刻 t での微分を意味するものとする。なお、−T 0 は
A1 を固 定した始点からの反作用であり、また T n = 0 であることに注意する。

4. 運動方程式4

次 は rj を中心とする回 転運動の方程式を考える。
Aj は単位ベクトル q(θj ) に垂直で、r̂j では q(θj ) 方向の力によって θ が増える
回 転をする (3)。よって、T j が r̂j で、rj を中心に Aj を θj の増える方向に回PSfrag replacements−T j−1r̂ j−1T j・q(θj )lj2−q(θj )rj−T j−1・(−q(θj ))θjlj2q(θj )Tjr̂ j

図 3: Aj に対する回転力のモーメント

そうとするモーメントは、
ljT j・q(θj )2の内積によって与えられる。同様に −T j−1 が r̂j−1 で、rj を中心に Aj を回 そうとするモーメントは、ljlj(−T j−1 )・(−q(θj )) = T j−1 )・q(θj )22となる。よって、Aj の重心での慣 性モーメントを Ij とすると、 Aj の rj 中心での回 転の方程式は次 のようになる。
Ij θ̈j =lj(T j + T j−1 )・q(θj )2(1 ≤ j ≤ n)(4)

慣 性モーメント Ij は、長さ lj と線密度 ρ により2ρ
x dx =
Ij = ρ
3−lj /2
Z lj /22lj
2!3=

mj lj2ρlj3=
1212(5)5. 相互力の消去

5となるので、(4) は
1m j ljθ̈j = (T j + T j−1 )・q(θj )122(1 ≤ j ≤ n)(6)となる。
(3) の運動方程式は、2 次 元 ベクトルなので、成分で考えれば 2n 本あり、(6)は n 本あるので合計 3n 本、未知関 数は T 0 ∼ T n−1 を成分で考えれば 2n 個、
θ1 ∼ θn が n 個なので、(3) と (6) から T 0 ∼ T n−1 を消去 すれば、θ1 ∼ θn に対する丁度 n 本の方程式が得られることになる。

5相互力の消去次 に、ニュートンの運動方程式 (3) と回 転の方程式 (6) から T 0 ∼ T n−1 を消去する。
(3) より、T j−1 − T j = −mj (r̈ j + gey )(1 ≤ j ≤ n)(7)となるが、これを j = k + 1 から j = n まで加 えると、 T n = 0 より
(T k − T k+1 ) + (T k+1 − T k+2 ) + · · · + (T n−1 − T n ) = T kとなるので、よって T k は
Tk = −nXmj (r̈ j + gey )(0 ≤ k ≤ n − 1)(8)

j=k+1と表される。そして、これにより
T j + T j−1 = −mj (r̈ j + gey ) − 2n
Xmi (r̈ i + gey )i=j+1
2= −nXnXd mi  gey
mj r j + 2
mi r i  − mj + 2
dt2
i=j+1
i=j+1(9)となり、これを (6) に代入すれば T j が消去 できることになる。その前に式 (9)を p(θk ) を用いて表しておく。そのために、
rj =n
Xk=1

ajk p(θk )(1 ≤ j ≤ n)(10)

5. 相互力の消去6と書くことにする。 (2) より、この ajk はlk l(k < j)jajk = (11)(k = j)20

(k > j)となる。また、
Lj =nXlj+lm
2 i=j+1(1 ≤ j ≤ n)(12)とすれば、(9) の右辺 2 項目の係 数は
mj + 2n
Xi=j+1

m i = ρ  lj + 2n
Xi=j+1li  = 2ρLj(13)と書ける。
(9) の右辺 1 項目の部分は、(10) により、
mj r j + 2n
X

mi r i = m ji=j+1=n
X

k=1n
X

ajk p(θk ) + 2p(θk ) mj ajk + 2n
Xi=j+1

mii=j+1

k=1nXmi aik  =n
Xn
Xaik p(θk )k=1

bjk p(θk )(14)

k=1の形 となるが、bjk は
bjk = mj ajk + 2nX

mi aiki=j+1

であり、k < j のときは (11) より
bjk = mj lk + 2nXmi lk = 2ρLj lk ,i=j+1

k = j のときは、n
Xn
X
lj
lj
bjj = mj ajj + 2
mi aij = mj + 2
mi lj = 2ρ Lj −2
4
i=j+1
i=j+1

k > j のときは、
bjk = 2mk akk + 2

n
X

i=k+1

mi aik = mk lk + 2n
X

i=k+1

mi lk = 2ρLk lk!

lj ,5. 相互力の消去

7となり、よって、
 2ρLj lk(k < j)bjk =  2ρ Lj − lj4!

!

lj= 2ρ Lj∨k − δjk lk4

lj (k = j)2ρLk lk(15)(k > j)となる。ここで、a ∨ b = max{a, b}, δjk = 0 (j 6= k), δjj = 1 とする。
(15), および (13) により、(9) は、
nX!

lj2ρ Lj∨k − δjk lk p̈(θk ) − 2ρLj gey
T j + T j−1 = −4
k=1(16)となる。ここで、
p̈(θk ) =d  0p (θk )θ˙k = p00 (θk )(θ˙k )2 + p0 (θk )θ¨k = −p(θk )(θ˙k )2 + q(θk )θ¨kdtであり、p(θk )・q(θj ) = (sin θk , − cos θk )・(cos θj , sin θj ) = − sin(θj − θk ),q(θk )・q(θj ) = (cos θk , sin θk )・(cos θj , sin θj ) = cos(θj − θk ),ey・q(θj ) = sin θj

なので、(16) より (6) は、
m j lj
1(T j + T j−1 )・q(θj )θ̈j =122!
nXlj= −ρLj∨k − δjk lk {(θ˙k )2 sin(θj − θk ) + θ¨k cos(θj − θk )}4k=1− ρLj g sin θjとなり、これを ρ で割 れば、
nX
lj2lj
θ̈j +
Lj∨k − δjk lk {θ¨k cos(θj − θk ) + (θ˙k )2 sin(θj − θk )}124
k=1+ Lj g sin θj = 0!となるが、この式の θj の微分の項をみると、
lj2ljθ̈j + Lj −124!

ljlj θ¨j = Lj −6!

lj θ¨j(17)5. 相互力の消去

8であり、また (17) の (θ̇j )2 の係 数は 0 になるので、結 局 (17) はnXk=1!

lj
Lj∨k − δjk lk {θ¨k cos(θj − θk ) + (θ˙k )2 sin(θj − θk )} + Lj g sin θj = 06(18)と書くことができる。これが、陽ではないが、θ1 ∼ θn に対する 2 階 の連立常微分方程式となる。
例えば、lj = l (1 ≤ j ≤ n)、すなわちすべての Aj の長さが等しい場合は、

Lj∨k −

1llj
δjk = l n − j ∨ k +− δjk =626

 

1

 l n−j∨k+ 21
 l n−j+3(j 6= k)(j = k)なので、この場合で n = 2 のときに方程式を書き下 すと、

3g4θ¨1 1 ¨˙ 2 3 + 2 {θ2 cos(θ1 − θ2 ) + (θ2 ) sin(θ1 − θ2 )} + 2l sin θ1gθ¨1
 {θ¨1 cos(θ1 − θ2 ) − (θ˙1 )2 sin(θ1 − θ2 )} + 2 +sin θ2232l=0(19)=0となる。

最後 に、(18) を行列の形 で書いておく。簡 単のため、
αjk = αjk (l1 , . . . , ln ) = Lj∨k −ljδjk
6(20)と書くことにすると、
A(θ) ="

lk αjk cos(θj − θk )#,B(θ) =j,k"lk αjk sin(θj − θk )#(21)
j,kに対して、(18) はθ¨(θ˙1 )2L sin θ1 1  1 .. 
.
..
A(θ)  .  + B(θ)  ..  + g 
=0.


2
θ¨n(θ˙n )Ln sin θn









(22)の形 となる。なお、αkj = αjk なので、lj = l (1 ≤ j ≤ n) の場合は、A(θ) は対称行列、B(θ) は交代行列となる。

6. 最後に69

最後に

今回 、摩擦のない多重剛体振り子の運動方程式の導出を、ラグランジュの運動方程式ではなく、ニュートンの運動方程式を用いて行った。 1 節にも書いたように、ネット上にはその方程式がほとんどなかったので、これを公開 することは意味がないわけでもないだろう。
なお、現 在 (2018 年 11 月 12 日) の日本語版の Wikipedia に「二重振り子」の項目があり ([1])、そこに n = 2 の場合の二重剛体振り子 (物理振り子) の方程式が書いてあるのだが、それは少し (19) の方程式とは係 数が違っている。 g の係 数と対角 成分の係 数が違っているのであるが、前者は記 載ミス、後 者は剛体棒の慣 性モーメントを考慮していない値ではないかと思われる。実際、大学 の図書館 で力学 の演 習 書を見 てみたところ、[2] p171、および [3] p177 には、同じ問題の解 として (19) と同じ式が書かれている。
なお、方程式 (22) の解 であるが、n = 2, 3 で |A(θ)| などを計算してみたところからすると、一般の n に対しても多分 |A(θ)| 6= 0 であると思われ、よって (22) は正規 形 の方程式になるので局 所解 の存在は当然言える。しかし、方程式に (θ˙j )2が含 まれるため、解 の大域的な存在を証明することはそれほど易しくはないと思われる。ただし、既 に論文等で方程式 (22) の解 の大域的な存在が否定されている (反例がある) のか、それとも大域的な存在が示されているのかは私は知らない。
ちなみに、私は全体の長さ nl = L を一定にして、n → ∞ にしたような物理モデ
ルである、いわゆる「くさり」の偏微分方程式に興 味があり、その方程式の解析、あるいは数値計算のための近 似 として n 重剛体振り子が使えないかと思って考え始めたのであるが、今のところその解 析も数値計算もあまりうまくいっていない。それは、n 重剛体振り子の方程式 (22) 自 身が陽の形 ではないために数値計算自 体が面倒なこともあるが、 n が多くなると先の方 (An ) の動きが激しくなってしまうのがその原 因ではないかと思う。そのような現 象は、「くさり」の方程式を単純に差分化 した場合も現 れるので、その意味では n 重剛体振り子が期 待通りにくさりの近 似 になっているのかもしれないが、結 局 どちらでもうまくいかないことになっていて、まだこちらの研 究 は先は長そうである。

参考文献[1] Wikipedia、「二重振り子」 : https://ja.wikipedia.org/wiki/二重振り子[2] 野上茂吉 郎、「基 礎物理学 選書 力学 演 習 」(1982)、裳華 房[3] 小西岳 、納繁男、「共 立物理学 講座、力学 演 習 」(1984)、共 立出版