二重振り子

2013/10/31追記：循環参照と行列演算を使った二重振り子シミュレーションExcelファイルup

学生で言う夏休みあたりの時期に、放送大学で「集中講義期間」と称して、1日2回のペースで全15回の講義をぶっ続けで行う期間があったんですよ。
そのうちごく一部を録画しまして
僕はどちらかというと物理寄りの人間なので
録画した講義の中には「力と運動の物理」っちゅうやつがあったんですね。


最初は簡単な内容だったんですよ。
実際に大学で習った程度の内容だったんです。
ところが中盤以降から急に見知らぬ内容に突入しましてね
ニュートンの運動方程式を拡張する必要があるとか、画面の先生が言い始めるんです。

それでラグランジュの運動方程式だとか、ハミルトンの運動方程式だとかの話になったんですけど
録画を見てた時点での僕は

「かえってややこしくしてるだけじゃん」とか「それやってなんになんのよ」
みたいな感想しか出てこなくてですね
それはもう退屈な、物理にしては退屈な授業だったんです。

「束縛力や慣性力もひっくるめた～」
なんてよくわからなくてですね


ラグランジアンやラグランジュのナンチャラっていうのは言葉としては聞いたことがありましたが
それ以上の何たるかをまったく習っていませんで

ハミルトニアンはかろうじて知ってましたよ。
量子力学に出てきますからね。
でもあれもわかりづらかったなぁ
シュレディンガー方程式がHψ＝Eψってたったのこれだけであらわされた日には
Hの中身が知りたいのになんで具体的な式かいてないんだよ！！！
って身悶えたもんでした。


ラグランジュ関係は人名と･･･数学？微分方程式近辺だったと思います
それからラグランジュ方程式っていう名前と数式を「知っていた」だけで、
理解はしてなかったような･･･


ニュートンの運動方程式をラグランジュの方程式に拡張する際にやった
運動エネルギーとポテンシャルの差？
が、解せなくてですね
なぜ同じ物理現象の力学を考えるのに和になったり差になったりするのかっていう感じですよ


あとで何度も授業を見返して、じっくりやらねばって思ってます。



まあそんな状態で
ふと、「二重振り子の問題」を解いてみようと思ったんですよ。
3体問題みたいに、「決定論的だけどもカオスになる」っていうのはわかってたんです。
ただ、「振り子を2つ連結するだけでカオスな動きを見せる」ように
運動方程式も、それ自体は作るのが簡単で、あとは解くのがめんどくさい
程度のことだと思ったんですね。

しかしいざ面と向かうと全然そんなことはなく
1つ目の振り子が2つ目の振り子に慣性で引っ張られるわ
張力はどうとか考えなきゃならんわで、立てることすら叶わぬ
ニュートンの運動方程式の非力さを痛感しましたよ･･･ええ。


ただ、それでラグランジュ方程式のありがたみがわかったところでラグランジュ方程式を使いこなせる自信もなく
とりあえずネットで解き方を見てみよう
ってなったんですが
30過ぎの初老にあのごちゃごちゃした数式はちょっと目が疲れます＞＜
合ってるかどうかもわかりませんし･･・


まあそんなわけで、途中過程を行列で解いているサイトを見つけたので
どうせ連立微分方程式を解いたところで解析的な解が出るわけでもなし
ここはもう最初から具体的に数値を入力して数値計算してしまえって感じですよ。



幸い、プログラムだとバグが全然見えなくてエクセルしか使えない僕でも
循環参照と行列計算という武器は持っているので
連立微分方程式は行列差分方程式にしてしまえば見える化がぐんと進んでミスも少なくすむわけです。


どこぞのサイトを頼りに
なんか加速度2つの列ベクトルP（α1、α2）があって、D（（ａ，ｂ）、（ｃ、ｄ））って2行2列の行列とR(ｐ、ｑ)って列ベクトルとの関係
DP=R
の関係があるそうなんですよ。

ただし
a=(m1+m2)L1^2
b=m2L1L2cos(s1-s2)
c=m2L1L2cos(s1-s2)
d=m2L2^2
p=-m2L1L2w2^2sin(s1-s2)-(m1+m2)gL1sins1
q=m2L1L2w1^2sin(s1-s2)-m2gL2sins2


で、mとLはそれぞれ振り子の重りの質量と振り子の長さ
添え字はそれぞれ1つ目と2つ目の振り子のものって意味で
g=9.8[m/s2]は重力加速度です。
例：m1は1つ目(根っこ)の振り子の質量、L2は2つ目（枝）の振り子の長さって感じです。


DP=RのPを求めたいんですから、両辺の左側からDの逆行列をかければPの中身である加速度α1とα2は求まります。
P=R/Dって感じなんですが、行列なんでP=invD*Rですね。


あ、加速度は角加速度です。単位はrad/s2

角速度[rad/s]をw1、w2とおいて
振り子の振れた角度(鉛直真下からの角度：rad)をs1、s2
として、○の時間微分を○’と記述すると

s'=w
w'=α
になります。

この微分方程式の微分を差分にして
ds/dt=Δs/Δt=(s(+)-s(-))/dt=w
ってすると
s(+)=s(-)+w・dt
ってなるじゃないすか。

こう計算すると循環参照しますよね？

なので、

初期値をそれぞれs10、s20、w10、w20として

s=if(初期化スイッチ=0,初期値,s自身+w・dt)
とかって感じで式4本を連結させると

s1=if(sw=0,s10,s1+w1*dt)
s2=if(sw=0,s20,s2+w2*dt)
w1=if(sw=0,w10,w1+α1*dt)
w2=if(sw=0,w20,w2+α2*dt)

αは行列計算から求めます。

みたいな感じで、あとは初期値とdtとスイッチswをテキトーに設定すれば
再計算のためにdelボタンでも押し続けると
逐次カオスな挙動不審な動きを再計算してくれるわけです。


グラフにしたければ
1つ目の振り子のx軸にx=L1sins1、y軸にy=L1coss1
2つ目の振り子のx軸にx=L1sins1+L2sins2、y軸にy=L1coss1＋L2coss2
とかってやって
振り子の質量と長さL1=L2=m1=m2=1、（dtは0.05とか適当に）って設定して
初期値にw2=1(枝だけつっつく)、ほかは全部0とかって設定してやればdel押しっぱなしで見ていて飽きないかもしれないカオス映像が画面に映りますよ。




ちなみに単純モデル化してるんで、振り子の棒同士の摩擦もないですし
棒の質量やら慣性モーメントも全然考慮してないです。
剛体とはいっても質量ゼロの剛体ですからね。


ただ、太さゼロの棒とか、大きさゼロの質点なんですが
摩擦はいつか考慮したいなぁと思うところだったりして。
やっぱ、二重振り子ってなんとなく、止まるまでが二重振り子のような気がしないでもないというか
動画として表現する際に、繰り返し表現あるいは初めと終わりがある表現のほうが好きって理由もあります。
このまま減衰がないとdel押し続ける限りいつまでも廻り続けちゃいますからねえ
その辺もパラメータに加えたいかもなあとか思ったり。

ぐぐってもいまいち出ないんですよね、摩擦ありのときの二重振り子の方程式。
実体がないからなのか、ラグランジュ方程式を使っても立てるのが難しいのか
あるいはあまり興味がないのか。

そういえば、二重振り子の極限として、m2=L2=0とかってして単振り子モードを表示しようとしたらなんかすぐに発散するんですよね。
速度(の1乗とか2乗)に比例する落下運動もなんかそんな感じでした。
空気抵抗のパラメータをゼロにすると発散するんですよ。
なんとかならないのかなぁこれ。

それはそれとして、摩擦ありの二重振り子の運動方程式
いつかラグランジュ方程式のやり方を理解して、シミュレーションしてみたいなぁ





メガネをかけた巨大ロボ

にほんブログ村