二重振り子のシミュレーション
こんにちは。おととしです。
前回、二重振り子の運動方程式の計算が終わったので、次にプログラムをしていきます。
プログラミングは少しだけやったことがあって、学生の頃にFortranとMTALABを用いていました。ただ、どちらもたくさんの人に使われている印象はないので、今回はPythonでやっていこうと思います。とはいっても、そもそもプログラミングは趣味程度しかやったことがないので、参考程度にご覧いただくと幸いです。
プログラムの流れとしては、
１．定数と変数の宣言
２．初期状態
３．θの更新式の定義
４．数値シミュレーションの更新式の定義（ルンゲ・クッタ法）
５．計算
６．プロット
というのを考えています。前回やったのは、３のθの更新式の定義だけです。
まず、定数と変数の宣言を行います。
import matplotlib.pyplot as plt import matplotlib.animation as anm import numpy as np # Const Definition g=9.8 # m/s2 m1=0.1 # kg m2=0.1 # kg l1=0.55 # m l2=0.45 # m mu=(m1+m2)/m2 dt=0.01 itr=1000 # Iteration Times theta1=[0]*itr theta2=[0]*itr V1=[0]*itr V2=[0]*itr X1=[0]*itr Y1=[0]*itr X2=[0]*itr Y2=[0]*itr
gは重力加速度、m1とm2は2つのおもりのそれぞれの質量、l1とl2は2つの振り子のそれぞれの紐の長さです。
muは次の計算のためにここで定義しておきます。
dtは時間の区切りで、細かいほどシミュレーションした時の動きが滑らかになります。
itrは計算の終点で、多いほど長い時間のシミュレーションを行えますが、計算にも時間がかかります。
theta1からY2までは計算結果を保持するための「箱」を用意しています。
次に、初期状態を定義しておきます。
# Initial State theta1[0]=np.pi/180*120 theta2[0]=np.pi/180*120 V1[0]=0 V2[0]=0 X1[0]=l1*np.sin(theta1[0]) Y1[0]=l1*np.cos(theta1[0]) X2[0]=X1[0]+l2*np.sin(theta2[0]) Y2[0]=Y1[0]+l2*np.cos(theta2[0]) t=0
二重振り子は最初の重りの位置でどのような挙動をするか変化します。今回は角度を120°にしています。
V1とV2はθの時間微分を表しています。
X1からY2までは、計算して得られるθの値を直交座標に変換しています。
次に、θの更新式の定義を行います。
# Function Definition def F_V1(g, mu, l1, l2, theta1, theta2, V1, V2): D_theta=theta1-theta2 F_V=(g*(np.sin(theta2)*np.cos(D_theta)-mu*np.sin(theta1))\ -(l1*V1*V1*np.cos(D_theta)+l2*V2*V2)*np.sin(D_theta))\ /(l1*(mu-(np.cos(D_theta))**2)) return F_V def F_V2(g, mu, l1, l2, theta1, theta2, V1, V2): D_theta=theta1-theta2 F_V=(g*mu*(np.sin(theta1)*np.cos(D_theta)-np.sin(theta2))\ +(mu*l1*V1*V1+l2*V2*V2*np.cos(D_theta))*np.sin(D_theta))\ /(l2*(mu-(np.cos(D_theta))**2)) return F_V
ここは、前回の記事で求めた式を間違いのないように書いたものとなります。
次に、今回のプログラムの肝となる、ルンゲ・クッタ法の記述となります。
# RK Definition def RK4(g, mu, l1, l2, theta1, theta2, V1, V2, dt, func1, func2): k1_x1=V1 k1_x2=V2 k1_v1=func1(g, mu, l1, l2, theta1, theta2, V1, V2) k1_v2=func2(g, mu, l1, l2, theta1, theta2, V1, V2) k2_x1=V1+dt*k1_v1/2 k2_x2=V2+dt*k1_v2/2 k2_v1=func1(g, mu, l1, l2, theta1+dt*k1_x1/2, theta2+dt*k1_x2/2, V1+dt*k1_v1/2, V2+dt*k1_v2/2) k2_v2=func2(g, mu, l1, l2, theta1+dt*k1_x1/2, theta2+dt*k1_x2/2, V1+dt*k1_v1/2, V2+dt*k1_v2/2) k3_x1=V1+dt*k2_v1/2 k3_x2=V2+dt*k2_v2/2 k3_v1=func1(g, mu, l1, l2, theta1+dt*k2_x1/2, theta2+dt*k2_x2/2, V1+dt*k2_v1/2, V2+dt*k2_v2/2) k3_v2=func2(g, mu, l1, l2, theta1+dt*k2_x1/2, theta2+dt*k2_x2/2, V1+dt*k2_v1/2, V2+dt*k2_v2/2) k4_x1=V1+dt*k3_v1 k4_x2=V2+dt*k3_v2 k4_v1=func1(g, mu, l1, l2, theta1+dt*k3_x1, theta2+dt*k3_x2, V1+dt*k3_v1, V2+dt*k3_v2) k4_v2=func2(g, mu, l1, l2, theta1+dt*k3_x1, theta2+dt*k3_x2, V1+dt*k3_v1, V2+dt*k3_v2) theta1+=(k1_x1+2*k2_x1+2*k3_x1+k4_x1)*dt/6 theta2+=(k1_x2+2*k2_x2+2*k3_x2+k4_x2)*dt/6 V1+=(k1_v1+2*k2_v1+2*k3_v1+k4_v1)*dt/6 V2+=(k1_v2+2*k2_v2+2*k3_v2+k4_v2)*dt/6 return theta1, theta2, V1, V2
数値シミュレーションでは、ルンゲ・クッタ法と呼ばれる方法がよく用いられているようです。運動方程式のように、微分方程式の形で表されているものの、解析解を求めるのが難しい場合に用いられます。漸化式のようになっていて、前のステップの値を用いて運動の傾きを近似し、次の位置を近似的に求めています。私も方法として用いているだけで、原理や精度については完璧には理解していません。
次に、実際に計算を行います。
# Calculation for i in range(itr-1): theta1[i+1], theta2[i+1], V1[i+1], V2[i+1]=RK4(g, mu, l1, l2, theta1[i], theta2[i], V1[i], V2[i], dt, F_V1, F_V2) # Convert to XY X1[i+1]=l1*np.sin(theta1[i+1]) Y1[i+1]=l1*np.cos(theta1[i+1]) X2[i+1]=X1[i+1]+l2*np.sin(theta2[i+1]) Y2[i+1]=Y1[i+1]+l2*np.cos(theta2[i+1])
ここまでに関数を定義していたので、計算の記述はシンプルになります。
最後に、プロットしていきます。
# Final Plot fig=plt.figure(figsize=(6,6)) def update(i, dt, X1, Y1, X2, Y2): if i !=0: plt.cla() if i<150: len=i else: len=150 T_x1=[0]*len T_y1=[0]*len T_x2=[0]*len T_y2=[0]*len for k in range(len): T_x1[k]=X1[i-len+k] T_y1[k]=Y1[i-len+k] T_x2[k]=X2[i-len+k] T_y2[k]=Y2[i-len+k] t=i*dt plt.plot([0, X1[i], X2[i]],[0, Y1[i], Y2[i]], color='k') # String plt.plot(0,0,marker='.', markersize=10, color='k') # Origin plt.plot(T_x1,T_y1, color='r') # m1 Trace plt.plot(T_x2,T_y2, color='b') # m2 Trace plt.plot(X1[i],Y1[i], marker='.', markersize=20, color='r') # m1 plt.plot(X2[i],Y2[i], marker='.', markersize=20, color='b') # m2 plt.xlim(-1.1,1.1) plt.ylim(1.2,-1.0) plt.title('t= %.1f s' %t) # Animation def animation_play(itr, dt, X1, Y1, X2, Y2): params={ 'fig':fig, 'func':update, 'fargs':(dt, X1, Y1, X2, Y2), 'interval':100, 'frames':itr, 'repeat':False, } ani=anm.FuncAnimation(**params) plt.show() animation_play(itr, dt, X1, Y1, X2, Y2)
アニメーションまでできるというので、調べてやってみました。また、質点の軌跡を描いていた例があって、キレイだなと思ったのでマネしてみました。
ここまでやってみて、この程度のシミュレーションであれば、思ったより簡単にできたという印象です。ただ、世の中に存在する物理現象は二重振り子だけではありません。それらを網羅していてはどれだけリソースがあっても足りないと感じました。いちいち運動方程式を立てて計算するのではなく、もっと違うアプローチがありそうだと考えています。今度はそっちの方法でいろいろシミュレーションできたらいいな、と思います。
コード全文：
import matplotlib.pyplot as plt import matplotlib.animation as anm import numpy as np # Const Definition g=9.8 # m/s2 m1=0.1 # kg m2=0.1 # kg l1=0.55 # m l2=0.45 # m mu=(m1+m2)/m2 dt=0.01 itr=1000 # Iteration Times theta1=[0]*itr theta2=[0]*itr V1=[0]*itr V2=[0]*itr X1=[0]*itr Y1=[0]*itr X2=[0]*itr Y2=[0]*itr # Initial State theta1[0]=np.pi/180*120 theta2[0]=np.pi/180*120 V1[0]=0 V2[0]=0 X1[0]=l1*np.sin(theta1[0]) Y1[0]=l1*np.cos(theta1[0]) X2[0]=X1[0]+l2*np.sin(theta2[0]) Y2[0]=Y1[0]+l2*np.cos(theta2[0]) t=0 # Function Definition def F_V1(g, mu, l1, l2, theta1, theta2, V1, V2): D_theta=theta1-theta2 F_V=(g*(np.sin(theta2)*np.cos(D_theta)-mu*np.sin(theta1))\ -(l1*V1*V1*np.cos(D_theta)+l2*V2*V2)*np.sin(D_theta))\ /(l1*(mu-(np.cos(D_theta))**2)) return F_V def F_V2(g, mu, l1, l2, theta1, theta2, V1, V2): D_theta=theta1-theta2 F_V=(g*mu*(np.sin(theta1)*np.cos(D_theta)-np.sin(theta2))\ +(mu*l1*V1*V1+l2*V2*V2*np.cos(D_theta))*np.sin(D_theta))\ /(l2*(mu-(np.cos(D_theta))**2)) return F_V # RK Definition def RK4(g, mu, l1, l2, theta1, theta2, V1, V2, dt, func1, func2): k1_x1=V1 k1_x2=V2 k1_v1=func1(g, mu, l1, l2, theta1, theta2, V1, V2) k1_v2=func2(g, mu, l1, l2, theta1, theta2, V1, V2) k2_x1=V1+dt*k1_v1/2 k2_x2=V2+dt*k1_v2/2 k2_v1=func1(g, mu, l1, l2, theta1+dt*k1_x1/2, theta2+dt*k1_x2/2, V1+dt*k1_v1/2, V2+dt*k1_v2/2) k2_v2=func2(g, mu, l1, l2, theta1+dt*k1_x1/2, theta2+dt*k1_x2/2, V1+dt*k1_v1/2, V2+dt*k1_v2/2) k3_x1=V1+dt*k2_v1/2 k3_x2=V2+dt*k2_v2/2 k3_v1=func1(g, mu, l1, l2, theta1+dt*k2_x1/2, theta2+dt*k2_x2/2, V1+dt*k2_v1/2, V2+dt*k2_v2/2) k3_v2=func2(g, mu, l1, l2, theta1+dt*k2_x1/2, theta2+dt*k2_x2/2, V1+dt*k2_v1/2, V2+dt*k2_v2/2) k4_x1=V1+dt*k3_v1 k4_x2=V2+dt*k3_v2 k4_v1=func1(g, mu, l1, l2, theta1+dt*k3_x1, theta2+dt*k3_x2, V1+dt*k3_v1, V2+dt*k3_v2) k4_v2=func2(g, mu, l1, l2, theta1+dt*k3_x1, theta2+dt*k3_x2, V1+dt*k3_v1, V2+dt*k3_v2) theta1+=(k1_x1+2*k2_x1+2*k3_x1+k4_x1)*dt/6 theta2+=(k1_x2+2*k2_x2+2*k3_x2+k4_x2)*dt/6 V1+=(k1_v1+2*k2_v1+2*k3_v1+k4_v1)*dt/6 V2+=(k1_v2+2*k2_v2+2*k3_v2+k4_v2)*dt/6 return theta1, theta2, V1, V2 # Calculation for i in range(itr-1): theta1[i+1], theta2[i+1], V1[i+1], V2[i+1]=RK4(g, mu, l1, l2, theta1[i], theta2[i], V1[i], V2[i], dt, F_V1, F_V2) # Convert to XY X1[i+1]=l1*np.sin(theta1[i+1]) Y1[i+1]=l1*np.cos(theta1[i+1]) X2[i+1]=X1[i+1]+l2*np.sin(theta2[i+1]) Y2[i+1]=Y1[i+1]+l2*np.cos(theta2[i+1]) # Final Plot fig=plt.figure(figsize=(6,6)) def update(i, dt, X1, Y1, X2, Y2): if i !=0: plt.cla() if i<150: len=i else: len=150 T_x1=[0]*len T_y1=[0]*len T_x2=[0]*len T_y2=[0]*len for k in range(len): T_x1[k]=X1[i-len+k] T_y1[k]=Y1[i-len+k] T_x2[k]=X2[i-len+k] T_y2[k]=Y2[i-len+k] t=i*dt plt.plot([0, X1[i], X2[i]],[0, Y1[i], Y2[i]], color='k') # String plt.plot(0,0,marker='.', markersize=10, color='k') # Origin plt.plot(T_x1,T_y1, color='r') # m1 Trace plt.plot(T_x2,T_y2, color='b') # m2 Trace plt.plot(X1[i],Y1[i], marker='.', markersize=20, color='r') # m1 plt.plot(X2[i],Y2[i], marker='.', markersize=20, color='b') # m2 plt.xlim(-1.1,1.1) plt.ylim(1.2,-1.0) plt.title('t= %.1f s' %t) # Animation def animation_play(itr, dt, X1, Y1, X2, Y2): params={ 'fig':fig, 'func':update, 'fargs':(dt, X1, Y1, X2, Y2), 'interval':100, 'frames':itr, 'repeat':False, } ani=anm.FuncAnimation(**params) plt.show() animation_play(itr, dt, X1, Y1, X2, Y2)
