二重振り子の運動
2011SE173 ： 望月大輔
指導教員：杉浦洋

はじめに1二重振り子に関するラグランジュの運動方程式をたてる．それを数値的に解き，運動の二重振り子のシミュレーションを行う．初期条件を変えてカオス現象を起こす条件を調べる．二重振り子の関数 L(t) を計算し，リアプノフ指数を推定する．リアプノフ指数と初期値の関係について調べる．

運動方程式

2

今回の二重振り子の運動方程式は，4.2(θ1 の場合)(m1 + m2 )l1 θ̈1 + m2 l2 θ̈2 cos(θ1 − θ2 )+ m2 l2 θ̇22 sin(θ1 − θ2 ) + (m1 + m2 )g sin θ1 = 0(θ2 の場合)l1 l2 θ̈1 cos(θ1 − θ2 ) + l22 θ̈2 − l1 l2 θ12 sin(θ1 − θ2 )+ gl2 sin θ2 = 0とする．
運動方程式を用いて近似線形方程式は，(m1 + m2 )l1 θ̈1 + m2 l2 θ̈2 + (m1 + m2 )gθ1 = 0l1 θ̈1 + l2 θ̈2 + gθ2 = 0となる．

3

固有振動数

振り角度が小さいときの運動の固有振動数を求めて運動シミュレーションを行う．
近似線形方程式より固有振動数の値は l1 = 2，l2 = 3，m1 = 2，m2 = 1．となる．θ1 (0) = θ2 (0) = 1◦ ，θ̇1 (0) =θ̇2 (0) = 0◦ /s とする．
0.04

0.02

51015

また，天気が予測不可能だとか，上昇する煙の粒子がどう運動するか予測できないとか，あるいは株価の変動は混沌としてるなどとも言う．カオスの数学的な概念の背後にあるものは予測不可能性である．しかし，力学系の不規則で複雑な振る舞いが広く関心を持たれて研究されるようになってからすでに 20 年以上過ぎたが，現在でもなお我々のカオスに対する本質的な理解は十分であるとは言い難い．実際，「カオス」ということに対する定義さえも未だに確立されていない [1]．
カオスの例

30 年近く前にアメリカの気象学者 エドワード・ローレンツは 決定論的乱流 という論文を発表し，天気のパターンは初期値に鋭敏であることをコンピュータ・シミュレーションにより示した．当時，初期値の非常に小さな変化がその後の天気に大きな影響を与えるという考えは驚くべき発見であった．エドワード・ローレンツは大気の状態を単純化して表した方程式をコンピュータで計算した．ある程度計算を進行させ，プリントアウトされた途中の計算結果を見たエドワード・ローレンツは，先程の計算と全く異なる結果になっていることを見て驚いた．当然エドワード・ローレンツもプログラムのバグや機械の故障を疑ったが，調べた結果そのようなものはなかった．エドワード・ローレンツが見た結果の大きな違いこそ，初期値鋭敏性だったのである [1]．
4.3

初期値鋭敏性

カオスとは大雑把に言えば力学系における軌道の「予測困難な」不規則な振る舞いのことであるとされている．
すなわち初期値が必然的に含むであろう誤差が力学系の時間発展と共に急速に増大し，長時間にわたる予測を困難にすることを指している．すなわち，初期値を少し変えただけで全く違った結果を生み出すことである [2]．

20-0.02-0.04

図 1 θ̇1 と θ̇2 を合わせた図 (縦軸が θ̇1 と θ̇2 ，横軸が時間)

4.4

リアプノフ指数

力学系においてごく接近した軌道が離れていく度合いを表す量である．定量的に表すと，位相空間内の２つの軌道の初期の距離を δZ0 としたとき，ある時点 ｔ における距離は次のように表せる [1]，[2]．

シミュレーションによれば，2 つの固有振動数の角速度| δZ(t) | ≈ eλt | δZ0 |(1)2πの周期は ω1 と ω2 とすると，2ω1 = ω2 となり，周期 ω1運動が表れた．
離れていく度合いは，初期のベクトルの向きの違いによっ
結果，特徴的な運動は発見されなかった．振り子時計は， ても変わってくる．従って，位相空間の次元数で分解しちゃんとした周期をもち動いている．
た リアプノフ指数のスペクトル が存在する．そのう

4

カオスについて4.1

カオスとはカオス（混沌）という言葉は日常会話でよく耳にする．
これは普通は秩序ないしは予測不可能の欠如を意味する．

ち最大リアプノフ指数は力学系の予測可能性を決定づける値である．正の最大リアプノフ指数は系がカオス的であることを示唆する基準の一つである．

5Mathematica による数値実験

二重振り子が行う運動のシミュレーションを Mathemat200
400
600
800
1000
ica 上に実現し，数値実験を行った．この問題で得られた，運動方程式に初期値を与える．質量，糸の長さ，固有振動数で求めた初期値を用いて，l1 = 2，l2 = 3，m1 = 2-5，m2 = 1．とし，重力加速度を g=9.8 とし，値を固定する．それ以外の変数 θ1 ，̇
θ1 ，
θ2 ，̇
θ2 のそれぞれに初期値を与え，二重振り子の変化について調べる．この問題の初期-10値の異なる解を θ = (θ1 (t)，
θ2 (t)，̇
θ1 (t)，̇θ2 (t))，̃
θ = (θ̃1 (t)，˜˜θ̃2 (t)，θ̇1 (t)，θ̇2 (t)) とする．初期値基準モデルと初期値に誤差を加えたモデル誤差の変化について詳しく分析する．今図 3 関数 L(T ) のグラフ．横軸は時間 t．
回は微分方程式の数値解法には Mathematica の NDSolveを用いた．
関数 L(t) が 0 に到達した．カオスが発生した．
カオス性の判定
正のリアプノフ指数を持つ．また 1000 秒実験して関数
θ1 (0) = θ2 (0) = 180◦ ，θ̇1 (0) = θ̇2 (0) = 0◦ /s の場合
L(t) が 0 に到達する．
δθ(t) = θ(t) − θ̃(t) とするとき，ある λ > 0 について式(1) より40∥δθ(t)∥ ≈ eλt ∥δθ(0)∥(2)∥δθ(t)∥ は誤差で，∥δθ(0)∥ は初期誤差である．
実験では δθ(0) の成分を区間 10−5 [−1, 1] の乱数として与え，L(t) = log ∥δθ(t)∥(3)の値を観察する．関数 L(t) のグラフの傾き λ > 0 としてリアプノフ指数が観察できる．
ただし振り子の運動は有限の空間とエネルギーにとどまるので，L(t) = 1 に到達した時点で観察を打ち切る．

6数値実験

3020105101520-10

図 4 関数 L(T ) のグラフ．横軸は時間 t．

振り子の振り幅が小さいときは，図 1 から図 3 で見たよ
図 6 のの直線は傾き 2.9 の直線であり，リアプノフ指数う現象は二つの固有振動の和に過ぎず，カオスは起こらな は λ=2.9 で，式 (1) の指数関数因子は e2.9t となる．
い．カオスが起こるとすれば，それは初期値に依存する．それを調べるために，基準解 θ(t) の初期値を θ1 (0) = θ2 (0) 7 おわりに，̇θ1 (0) = θ̇2 (0) = 0 として，実験する．
本研究では，二重振り子のカオスについて観察した．リ
θ1 (0) = θ2 (0) = 69.8◦ ，θ̇1 (0) = θ̇2 (0) = 0◦ /s の場合
アプノフ指数が指数増加していると，初期値鋭敏性が発関数 L(t) が 0 に到達しなかった．カオスが発生しない． 生して，二重振り子はカオスの現象を発生することが分かった．今回の結果より，69.8◦ 以下の場合カオスの現象-4は発見されなった．69.9◦ になるとカオスが発生する．更-6に，180◦ のとき最大リアプノフ指数であった．カオスが発生したとき位置と速度の初期値を決めても未来を予測-8することは出来なかった．
今後の課題は，リアプノフ指数では L(t) が 0 以上にな-10ればこの解では必然的に L(t) の増大が止まり，カオス的-12な状態とそうでない状態を区別できなくなる．そのよう200400600
8001000
な状況において小藤俊幸教授は解の周波数解析を行うよ
うに勧めてくれたが，時間がなくてそこまで進むことが図 2 関数 L(T ) のグラフ．横軸は時間 tできなかった．その方向での研究が今後の課題である．

θ1 (0) = θ2 (0) = 69.9◦ ，θ̇1 (0) = θ̇2 (0) = 0◦ /s の場合

参考文献[1] デニー・グー・リック，『カオスとの遭遇』，産業図書，
1995．

図 5 の直線は傾き 0.0155 の直線であり，リアプノフ指数『力学系の基礎』，朝倉書店，2000．
は λ=0.0155 で，式 (1) の指数関数因子は e0.0155t となる． [2] 國府寛司，