数式処理ソフト DERIVE(デライブ) de ドライブ

81.連立常微分方程式の数値解法(2重振り子)

(1)2重振り子の運動(数値計算)

「数式処理ソフト DERIVE(デライブ)の第80回では、単振り子の運動を数値計算と解析的計算の両方で調べて見た。
 今回は、2重振り子の運動を考えてみよう」

「わかりやした。
 単振り子の場合の図に似せて描くとこんな感じになる。

 支持点 A(原点Oからの高さ= r+s)から垂直につり下げられた2重振り子の振動を考えます。
 振り子の重り自体は、その大きさを無視できるものとして、ま、要するに質点と考えるわけですね。
 また、糸は、質量を無視できる細い糸とします。
 糸の長さは、図のように、r とs、重りの質量は、2つとも同じとし、m とします。
 位置エネルギー U は、O点の位置でのエネルギーを、ゼロ とすると、
 U=mg(r+s-r cosθ)+mg(r+s-r cosθ-s cosφ)、
 Pの運動エネルギーは、(m/2)r^2θ^2、
 Qの位置は、直交座標による表現は、(r sinθ+s sinφ,r+s-r cosθ-s cosφ)、
 この位置座標のx成分とy成分をそれぞれ時間で微分して、2乗和を作ると、
 r^2θ’^2+s^2φ’^2+2 r sθ’φ’cos(θ-φ)、となる。
 よって、運動エネルギー T は、
 T=(m/2)r^2 θ'^2+(m/2)(r^2θ’^2+s^2φ’^2+2 r sθ’φ’cos(θ-φ))、
 ラグラジアンをLとすると、L=T-U、 オイラーの方程式は、
 d/dt(∂(L,θ'))-∂(L,θ)=0、
 d/dt(∂(L,φ'))-∂(L,φ)=0、となるので、
 まとめると、
 kφ''(t)·COS(θ(t) -φ(t)) + kφ'(t)^2·SIN(θ(t) -φ(t)) + 2ω^2·SIN(θ(t)) + 2θ''(t) = 0、
 θ''(t)·COS(θ(t) -φ(t)) -θ'(t)^2·SIN(θ(t) -φ(t)) +ω^2·SIN(φ(t)) + kφ''(t) = 0、
 となります。
 ここで、k≡s/r、ω^2≡g / r、と書いている」

「ご苦労さんじゃった。
 具体例で、数値計算を行ってみよう。
 ともちゃんが導いてくれた2つの式から、θ' ' とφ' '、を解くと、
 p=dθ/dt、q=dφ/dt、と書くことにすれば、
 dp/dt=((COS(2·θ - 2·φ)·(k·q^2·SIN(θ - φ) - 2·ω^2·SIN(θ)) - SIN(2·θ - 2·φ)·(k·q^2·COS(θ - φ) + 2·p^2) + 4·ω^2·SIN(φ)·COS(θ - φ) - 4·ω^2·SIN(θ)·SIN(θ - φ)^2 - 3·k·q^2·SIN(θ - φ) - 6·ω^2·SIN(θ))/(4·(SIN(θ - φ)^2 + 1)))、
 dq/dt=((k·q^2·SIN(2·θ - 2·φ) + 4·ω^2·SIN(θ)·COS(θ - φ) + 4·p^2·SIN(θ - φ) - 4·ω^2·SIN(φ))/(2·k·(SIN(θ - φ)^2 + 1)))、
 じゃ。
 m次までのオイラー展開法を適用するためには、
 k=1(上下の糸の長さが等しい)、ω^2=98、(r=0.1 m、重力加速度g=9.8 m/sec^2)、として、
 初期位置は、t=0の時に、θ=φ=π/6、初期速度はゼロ、p=0、q=0 、刻み幅をh=0.01として、
 Eulerm_user_R([1,p,q,pの右辺,qの右辺],[t,θ,φ,p,q],[0,π/6,π/6,0,0],0.01,4,1000)、とすればよい。
 しかし、p ' 及びq ' の表式に三角関数が含まれる複雑な式であるため、テイラー展開項も複雑になるようで、このままでは、計算が終わらないようじゃ。
 そこで、ルンゲ・クッタ法に頼るとしよう。
 DERIVEのルンゲ・クッタ法関数は、組み込み関数、RK( )であるが、特別なユーティリティ関数を呼び出す必要は無く、
 RK([p,q,pの右辺,qの右辺],[t,θ,φ,p,q],[0,π/6,π/6,0,0],0.01,1000)のように引数を与えればよい。
 なお、わしらが作った、Eulerm_user_R()と異なり、第1引数と第2、3引数の次元が1つ異なるので注意するようにな。
 k=1、ω^2=98、、初期角度 θ=φ=π/6、初期速度はゼロのグラフは下図の通りじゃ。


 これを見ると、上の振り子の波形(θ=赤)と下の振り子(φ=青)の周期は、同じであるが、波形は、純粋のコサインカーブから少し崩れている。
 また、振幅は、下の振り子の方が上の振り子に比べて、変動が大きいようじゃ」

「じゃ、わたしは、下の振り子の糸の長さが長い、k=2、(r=0.1、s=0.2)の場合を計算してみるよ。
 その他のパラメーターは、同じとします。


 この場合は、上の振り子の波形(赤)は、k=1の場合より、崩れている感じね。
 一方、下の振り子(青)は、多少の崩れはあるもののコサインカーブに近づいてきた。

 では、k=3の場合がこれ。
 
 下の振り子は、ほぼ、コサインカカーブ。
 上の振り子の波形は、崩れているなりに、左右対称になっている。

 最後は、下の振り子の糸が短い場合、k=0.5、(r=0.1、s=0.05)その他の諸元は、同じ。
 

 今度は、上、下とも振り子の波形は、整ってきている感じだわ」

「では、わしは、下の振り子の糸を、k=0.1、(r=0.1、s=0.01)と、さらに短くしてみよう。
 

 上、下の振り子の振幅、周期ともほぼ等しいものとなった。波形は、きれいじゃ。
 次節で、解析的な扱いを行って、みよう」

(2)2重振り子の運動(線形近似)

「2重振り子の運動方程式は、p=θ'、q=φ'、と書いたとき、
 dp/dt=((COS(2·θ - 2·φ)·(k·q^2·SIN(θ - φ) - 2·ω^2·SIN(θ)) - SIN(2·θ - 2·φ)·(k·q^2·COS(θ - φ) + 2·p^2) + 4·ω^2·SIN(φ)·COS(θ - φ) - 4·ω^2·SIN(θ)·SIN(θ - φ)^2 - 3·k·q^2·SIN(θ - φ) - 6·ω^2·SIN(θ))/(4·(SIN(θ - φ)^2 + 1)))、
 dq/dt=((k·q^2·SIN(2·θ - 2·φ) + 4·ω^2·SIN(θ)·COS(θ - φ) + 4·p^2·SIN(θ - φ) - 4·ω^2·SIN(φ))/(2·k·(SIN(θ - φ)^2 + 1)))、じゃった。
 (1)節の数値計算では、上の振り子の糸の長さ r を0.1 m、重力加速度gを9.8m/sec^2、として、ω^2=g/r、から、ω^2=98と置いて計算した。
 初期角度は、θ=φ=π/6≒0.5235・・、であるので、θ、φの一次までとることにする。
 同様に、p、q の2乗等の項もゼロと見なす。
 このような線形近似で、前節の結果をどの程度、説明できるかを調べて見よう。
 上掲の複雑な2つの式は、近似的に、
 2·θ''(t) + k·φ''(t) + 2·ω^2·θ(t) = 0、
 θ''(t) + k·φ''(t) + ω^2·φ(t) = 0、となることが分かる。
 初期条件は、θ(0)=φ(0)=a、θ' (0)=φ' (0)=0、じゃな。
 ともちゃん、ラプラス変換法により、この2つの方程式の初期条件を満たす解(特解)を求めてご覧」

「待ってました。
 θ(t)のラプラス変換をf(s)、φ(t)のそれを、g(s)としますね。
 ここで、sは、ラプラス変換の変数とする。(2つめの振り子の糸の長さではないので注意してね)
 初期条件を勘案して、(ωの代わりにwを使っている)
 s^2·f - s·a + w^2·f + k/2·(s^2·g - s·a) = 0、
 s^2·f - s·a + w^2·g + k·(s^2·g - s·a) = 0、
 整理すると、
 (s^2+w^2)f+(k/2)s^2 g=s*a(1+(k/2))、
 s^2 f+k s^2 g=s*a(1+k)、
 となる。
 この連立方程式から、f、g、を解くと、
 f = a·s·(k·(s^2 + w^2) + 2·w^2)/(k·s^2·(s^2 + 2·w^2) + 2·w^2·(s^2 + w^2))、
 g = a·s·(k·(s^2 + 2·w^2) + 2·w^2)/(k·s^2·(s^2 + 2·w^2) + 2·w^2·(s^2 + w^2))、
 この逆変換を求めれば、特解が得られる。
 一見すると、分母は、sの4次式なんだけど、s^2の2次式に過ぎないので因数分解が可能。
 その上で、部分部数展開を行う、というのが定番ね。
 と、まあ、原理は簡単なんだけど、係数がややこしいので、DERIVEの「代数」メニューの「展開」で、変数をs、「分母を実数の範囲で因数分解」を指示すると、
 f=a·k·s·(√(k^2 + 1) + 1)/(2·√(k^2 + 1)·(k·s^2 - w^2·(√(k^2 + 1) - k - 1)))+a·k·s·(√(k^2 + 1) - 1)/(2·√(k^2 + 1)·(k·s^2 + w^2·(√(k^2 + 1) + k + 1)))、
 が得られる。
 ここで、一般に、cos(αt)のラプラス変換が、s/(s^2-α^2)、であることを使うと、
 θ(t)=a·k·(√(k^2 + 1) + 1)/(2·k·√(k^2 + 1))·COS(w·t·√(k + 1 - √(k^2 + 1))/√k) + a·k·(√(k^2 + 1) - 1)/(2·k·√(k^2 + 1))·COS(w·t·√(√(k^2 + 1) + k + 1)/√k)、
 が得られる。
 同様に、φについても、
 φ(t)=a·k·(√(k^2 + 1) + k + 1)/(2·k·√(k^2 + 1))·COS(w·t·√(k + 1 - √(k^2 + 1))/√k) + a·k·(√(k^2 + 1) - k - 1)/(2·k·√(k^2 + 1))·COS(w·t·√(√(k^2 + 1) + k + 1)/√k)、
 と求まったわ」

「どれどれ、検算してみるかの。
 θとφを、線形近似式に代入して計算させると、2つの式のどちらも、満足していることが分かった。
 ついでに言うと、DERIVEでのラプラス逆変換は、ユーザー定義ファイルが下記に置かれている。
 C:\Program Files (x86)\TI Education\Derive 6\Users\Transforms、内にある「LaplaceTransforms.dfw」を使うこともできる。
 逆変換表を目で引くのが面倒なときに使ってよいじゃろう。
 たとえば、InvLaplace(f(s))、と計算させれば、f(s)の逆変換 θ(t)が求められる。
 近藤先生の本(「演算子法」(近藤次郎 著:培風館:1970年6月30日第13刷)の逆変換表ほどの結果は、得られないかもしれんがの。
 さてと、こうして、線形近似解が求まったので、(1)節の数値計算結果をどの程度、再現できているかを検証してみよう。
 初期角度 a=π/6、w=√(98)、は、共通なので、最初に置いてしまう。
 残されたパラメーターは、k (第2の振り子の糸の長さ/第1の振り子の糸の長さ)のみになる。
 k=1の結果は、次のようになる。
 θ=·(√2/24 + 1/12)·COS(7·t·√(4 - 2·√2)) + ·(1/12 - √2/24)·COS(7·t·√(2·√2 + 4))、
 φ=·(√2/12 + 1/12)·COS(7·t·√(4 - 2·√2)) + ·(1/12 - √2/12)·COS(7·t·√(2·√2 + 4))、
 θを赤、φを青の線で表すと、グラフは下図の通り。
 
 
 一方、、k=1のRK法による計算結果を再掲したものが下図じゃ。
 

 よく似ているようじゃが、時間 t が大きくなると、線形近似は、(1)節の計算結果とずれてくる傾向がある。
 その点を検証するために、下図では、線形近似値と前節のRK法による計算値との比較じゃ。
 k=1、ω^2=98、初期角度π/6、初期速度ゼロ、線形近似値-RK計算値と時間との関係を示している。
 
 θ-RK計算値を水色、φーRK計算値をグレーで示す。
 どちらも、t が大きくなるにつれて、線形近似値は、RK法による計算値から外れ、その差が徐々に拡大している。

 ついでに言うておくと、
 線形近似解は、θ=A cos(λt)+B cos(μt)、φ=C cos(λt)+D cos(μt)、の形をしている。
 すなわち、
 A=a·k·(√(k^2 + 1) + 1)/(2·k·√(k^2 + 1))、B=a·k·(√(k^2 + 1) - 1)/(2·k·√(k^2 + 1))、
 C=a·k·(√(k^2 + 1) + k + 1)/(2·k·√(k^2 + 1))、D=a·k·(√(k^2 + 1) - k - 1)/(2·k·√(k^2 + 1))、
 λ=w√(k + 1 - √(k^2 + 1))/√k、
 μ=w√(√(k^2 + 1) + k + 1)/√k、
 λとμの振動は、「基準振動」と呼ばれる。
 これから、A/B=(√(k^2 + 1) + 1)^2/k^2、
 また、C/D=- (√(k^2 + 1) + k + 1)^2/(2·k)、
 グラフで示すと下図の通り。
 
 kが小さいときは、θ、φとも、λの角周波数の振幅が支配的になる。
 k=0.1の場合がそれじゃな。
 逆にkが大きくなると、θでは、λとμの寄与がほぼ等しくなる。k=3の場合が、そのようなケースじゃ。
 φについては、kの全域でλの振動が支配的じゃ。(k=1のとき、C/Dの絶対値が最小となるが、その値は、約5.82)」

「θとφの図形を描いてみたよ。[θ,φ]ね。
 
 横軸がθ、縦軸がφ、t=0~10まで。
 λとμの比率は、無理数なので、厳密には、θとφの値が同じ組み合わせは、起きないのよね。
 もし、λ=1、μ=3、とすれば、下図のようになる。(t=0~100)
 この場合は、tをいくら増やしても、このカーブ以外にずれることがない。
 
 
 λ=1、μ=1/3、の場合は、下図のよう。(t=0~100)
 この場合も、これ以上、tを増やしても変わらない。
 
 
 とはいえ、有理数の比の場合も、λ=5/13、μ=3/11、などでは、下図のよう。(t=0~100)。
 やや複雑になってくる。
 
 
 この場合は、tを1000までとると、下のように密になる。
 
 しかし、tをこれ以上増やしても、四角形の領域内が、これ以上、塗りつぶされることはない。

 一方、λ=1、μ=√2 では、t=100までで、下図のように、上の有理数比の場合と変わらないように見える。
 
 
 けれども、t=1000まで増やしてみると、下図のようになる。
 
 t=10000まで経過するとこの解像度では、白い部分がなくなってしまう」

「角周波数の比率が有理数か、無理数かにより、同じ軌道を通るか、否かが決まっていることがよく分かるのう。
 比率が無理数の場合は、許される範囲内をまんべんなく、通り過ぎるので、t を増やしていくと、べた塗りになってしまうのじゃな。
 線形近似解についての検討は、不十分かもしれんが、これくらいにしておこう。
 (1)節の数値計算では、(初めて)わしらのオイラー展開法では、うまいかなかったな。
 何故うまくなかったかを調べて見ると、テイラー展開項を求める「Difm_user_R」関数で、2次までは、求められるんじゃが、3次以上では、途方もなく時間がかかってしまう。
 ルンゲ・クッタ法(4次)は、DERIVEに組み込み関数としてあるので、次回は、多段階法の「(陽的)アダムス公式」を取り上げてみよう」 

最終更新日 2016/4/25