共鳴する振り子
例の如く PostScript で書いて ImageMagick で変換。
「どんな時に共鳴が起こるかな? Let's パラメータをいろいろ変えてみよう!」
実際に下のプログラムを実行してみてくれる人なんているのだろうか。
以下ソース
%!PS gsave % 水平方向に l1 だけ離れた 2 個の点に、長さ l2 の棒がそれぞれ繋がれている。 % 棒の先端は長さ l1 の棒で繋がれている。 % 長さ l2 の棒の先には長さ l3 の棒が繋がれ、その先には質量 m3 の重りがある。 % もう片方の長さ l2 の棒の先には、長さ l4 の棒と質量 m4 の重りがある。 % 棒のつなぎ目は滑らかに動く。 % 長さ l1 の棒の質量は m1、その他の棒の質量は無視できる。 /l1 60 /l2 55 /l3 75 /l4 75 /m1 100 /m3 10 /m4 10 7 {def} repeat % 重力加速度を g とする。 /g 40 def % 棒 l2,l3,l4 の鉛直方向に対する傾きを th2,th3,th4 とする。 % 初期値 /th2 0 /th3 -0.3 /th4 0 3 {def} repeat % th2,th3,th4 の時間微分を w2,w3,w4 とする。 % 初期値 /w2 0 /w3 0 /w4 0 3 {def} repeat % 重りの半径 /r 10 def % 時間ステップ dt とステップの数 stepnum /dt 0.2 def /stepnum 420 def % ラジアン単位の三角関数 /rsin {3.141592653589793 div 180 mul sin} def /rcos {3.141592653589793 div 180 mul cos} def % 棒と重りの位置は % m1(左端): (x,z)=( l2*sin(th2), -l2*cos(th2) ) % (右端): (x,z)=( l1+l2*sin(th2), -l2*cos(th2) ) % m3: (x,z)=( l2*sin(th2)+l3*sin(th3), -l2*cos(th2)-l3*cos(th3) ) % m4: (x,z)=( l1+l2*sin(th2)+l4*sin(th4), -l2*cos(th2)-l4*cos(th4) ) /m1xl {l2 th2 rsin mul} def /m1xr {m1xl l1 add} def /m1z {l2 th2 rcos mul neg} def /m1xzl {m1xl m1z} def /m1xzr {m1xr m1z} def /m3x {l2 th2 rsin mul l3 th3 rsin mul add} def /m3z {l2 th2 rcos mul l3 th3 rcos mul add neg} def /m3xz {m3x m3z} def /m4x {l2 th2 rsin mul l4 th4 rsin mul add l1 add} def /m4z {l2 th2 rcos mul l4 th4 rcos mul add neg} def /m4xz {m4x m4z} def % 質量の速度は % m1: (vx, vz)=( l2*w2*cos(th2), l2*w2*sin(th2) ) % m3: (vx, vz)=( l2*w2*cos(th2)+l3*w3*cos(th3), % l2*w2*sin(th2)+l3*w3*sin(th3) ) % m4: (vx, vz)=( l2*w2*cos(th2)+l4*w4*cos(th4), % l2*w2*sin(th2)+l4*w4*sin(th4) ) % 運動エネルギー K と位置エネルギー V % K=1/2*(m1+m3+m4)*l2^2*w2^2 % +1/2*m3*(l3^2*w3^2+2*l2*l3*w2*w3*cos(th2-th3)) % +1/2*m4*(l4^2*w4^2+2*l2*l4*w2*w4*cos(th2-th4)) % V=-(m1+m3+m4)*g*l2*cos(th2)-m3*g*l3*cos(th3)-m4*g*l4*cos(th4) % Lagrangian % L=K-V % 運動方程式は d/dt(dL/dwi)-dL/dthi=0 (i=2,3,4) より % A: % d/dt[(m1+m3+m4)*l2^2*w2+m3*l2*l3*w3*cos(th2-th3)+m4*l2*l4*w4*cos(th2-th4)] % + m3*l2*l3*w2*w3*sin(th2-th3) + m4*l2*l4*w2*w4*sin(th2-th4) % + (m1+m3+m4)*g*l2*sin(th2) =0 % B: % d/dt( m3*(l3^2*w3+l2*l3*w2*cos(th2-th3)) ) % - m3*l2*l3*w2*w3*sin(th2-th3) + m3*g*l3*sin(th3) =0 % C: % d/dt( m4*(l4^2*w4+l2*l4*w2*cos(th2-th4)) ) % - m4*l2*l4*w2*w4*sin(th2-th4) + m4*g*l4*sin(th4) =0 % 変形する % mu1=m1/(m1+m3+m4), mu3=m3/(m1+m3+m4), mu4=m4/(m1+m3+m4) とすると m1 m3 m4 add add dup dup m1 exch div /mu1 exch def m3 exch div /mu3 exch def m4 exch div /mu4 exch def % A: % l2*dw2/dt +mu3*l3*cos(th2-th3)*dw3/dt +mu4*l4*cos(th2-th4)*dw4/dt % = -mu3*l3*w3^2*sin(th2-th3)-mu4*l4*w4^2*sin(th2-th4) -g*sin(th2) % B: % l3*dw3/dt +l2*cos(th2-th3)*dw2/dt % = l2*w2^2*sin(th2-th3) -g*sin(th3) % C: % l4*dw4/dt +l2*cos(th2-th4)*dw2/dt % = l2*w2^2*sin(th2-th4) -g*sin(th4) % A に B,C を代入 % l2*(mu1+mu3*sin(th2-th3)^2+mu4*sin(th2-th4)^2)*dw2/dt % = -mu3*sin(th2-th3)*(l2*w2^2*cos(th2-th3)+l3*w3^2+g*cos(th3)) % -mu4*sin(th2-th4)*(l2*w2^2*cos(th2-th4)+l4*w4^2+g*cos(th4)) % -g*mu1*sin(th2) % [th2 th3 th4 w2 w3 w4] から [w2 w3 w4 dw2/dt dw3/dt dw4/dt] を % 計算する手続き f /f { 9 dict begin % [/th2 /th3 /th4 /w2 /w3 /w4] exch arraydef aload pop [/w4 /w3 /w2 /th4 /th3 /th2] {exch def} forall /dw2dt g mu1 mul th2 rsin mul neg mu3 th2 th3 sub rsin mul l2 w2 dup mul mul th2 th3 sub rcos mul l3 w3 dup mul mul add g th3 rcos mul add mul sub mu4 th2 th4 sub rsin mul l2 w2 dup mul mul th2 th4 sub rcos mul l4 w4 dup mul mul add g th4 rcos mul add mul sub mu1 th2 th3 sub rsin dup mul mu3 mul add th2 th4 sub rsin dup mul mu4 mul add l2 mul div def /dw3dt w2 dup mul th2 th3 sub rsin mul l2 mul g th3 rsin mul sub th2 th3 sub rcos l2 mul dw2dt mul sub l3 div def /dw4dt w2 dup mul th2 th4 sub rsin mul l2 mul g th4 rsin mul sub th2 th4 sub rcos l2 mul dw2dt mul sub l4 div def [w2 w3 w4 dw2dt dw3dt dw4dt] end } def % 絵を描く手続き /drawpendulum { -0.5 l1 mul 75 add 150 translate 0 0 moveto m1xzl lineto m3xz lineto stroke l1 0 moveto m1xzr lineto m4xz lineto stroke 1 0 0 setrgbcolor 3 setlinewidth m1xzl moveto m1xzr lineto stroke 0 1 0 setrgbcolor m3xz r 0 360 arc fill 0 0 1 setrgbcolor m4xz r 0 360 arc fill showpage } def % 配列の加算、スカラー倍、配列を使ってまとめて代入 /arrayadd { 3 dict begin /a2 exch def /a1 exch def [ 0 1 a1 length 1 sub { /i exch def a1 i get a2 i get add } for ] end } def % /arraymul { % 3 dict begin % /c exch def % /a1 exch def % [ % 0 1 a1 length 1 sub { % /i exch def % a1 i get c mul % } for % ] % end % } def /arraymul { 1 dict begin /c exch def [exch {c mul} forall] end } def % /arraydef { % 4 dict begin % /a2 exch def % /a1 exch def % /len a1 length def % 0 1 len 1 sub { % /i exch def % a1 i get a2 i get % } for % len % end % {def} repeat % } def % 2次の Runge-Kutta /c 0.5 def /b 0.5 c div def /a 1 b sub def 0 1 stepnum { % 2 step 毎に図を表示 2 mod 0 eq { drawpendulum % /skip 0 def } if /th_and_w [th2 th3 th4 w2 w3 w4] def /dth_and_dw1 th_and_w f dt arraymul def /dth_and_dw2 th_and_w dth_and_dw1 c arraymul arrayadd f dt arraymul def th_and_w dth_and_dw1 a arraymul arrayadd dth_and_dw2 b arraymul arrayadd % [/th2 /th3 /th4 /w2 /w3 /w4] exch arraydef aload pop [/w4 /w3 /w2 /th4 /th3 /th2] {exch def} forall } for grestore
追記:forall を使うと arraydef は要らなくなる。ので書き換える。arraymul も簡単になる。
