RuntaScience diary

気象系データを扱う学生 旅が好きです

Welcome to my blog

About

Matplotlibで振り子ー単振り子

こんにちは

今日は関数を使って、単振り子をアニメーションを作成したいと思います!

 アニメーションの基礎は前の記事を参照してください!

 

runtascience.hatenablog.com

 

 

 

微分方程式

単振り子を描写するのに関数が必要なので書き方を確認しましょう

docs.scipy.org

モジュール

 

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np

 

解く

それでは、solve_ivpで指数関数的減衰を解いてみます。

 指数関数的減衰は半減期の計算で使われますね

https://www.nli-research.co.jp/report/detail/id=58700?site=nli

 

dy dt = -λ y

 

# 指数関数的減衰
def exponential_decay(t, y):
    return -0.5 * y

t_span = [0, 10] #tの範囲
y0 = [2, 4, 8] #y初期値
t_eval = np.arange(t_span[0], t_span[-1], 0.1) #tの間隔
sol = solve_ivp(exponential_decay, t_span, y0, t_eval=t_eval) #方程式を解く

 

解(sol)の中身を確認します。

print(sol.t)
⇒[0. 0.1 0.2 0.3…9.7 9.8 9.9]
print(sol.y)
⇒[[2. 1.90245885 1.809668…0.01492886 0.01420166]
[4. 3.8049177 3.619336…0.02985772 0.02840332]
[8. 7.6098354 7.238672…0.05971543 0.05680663]]

 tに対するyが作成されていることがわかります。

 

描写

それでは描写してみます

fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0])
ax.plot(sol.t, sol.y[1])
ax.plot(sol.t, sol.y[2])

ax.set_xlabel("T")
ax.set_ylabel("Y-axis")

ax.set_xlim(t_span[0],t_span[-1])
ax.set_ylim(0, 9)

plt.show()

#保存
#fig.savefig("XXX.png", format="png", dpi=330)

 

f:id:RuntaScience:20200621150622p:plain

 

 

単振り子

それでは単振り子のアニメーション作成してみます。

以下を参考に作成しました。

qiita.com

 

モジュール

まずはモジュールをインポートします

from scipy.integrate import solve_ivp
from matplotlib.animation import FuncAnimation
import matplotlib.pyplot as plt
import numpy as np

 

 

単振り子

 

%matplotlib nbagg
G = 9.8  # 重力加速度(m/s^2)
L = 1.0  # 振り子の長さ(m)

# 運動方程式を解く
def eom(t, state):
    dydt = np.zeros_like(state)
    dydt[0] = state[1]
    dydt[1] = -(G/L)*np.sin(state[0])
    return dydt

#設定
t_span = [0,20] #時間 (s)
dt = 0.05# 時間の間隔 (s)
t = np.arange(t_span[0], t_span[-1], dt)# 時間の間隔のarray作成

#初期状態設定
th1 = 45.0# 角度 [deg]
w1 = 0.0 # 初期角速度 [deg/s]
state = np.radians([th1, w1])# 初期状態

#微分方程式を解く
sol = solve_ivp(eom, t_span, state, t_eval=t)
theta = sol.y[0] #y[0]=θ、y[1]=ω

x = L * np.sin(theta)       
y = -L * np.cos(theta)      

#ーーーーーーグラフーーーーーーー
fig, ax = plt.subplots()

line, = ax.plot([], [], "o-", linewidth=2) 
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, " ", transform=ax.transAxes)

def init():
    line.set_data([], [])
    time_text.set_text("")
    return line, time_text

def animate(i):
    thisx = [0, x[i]]
    thisy = [0, y[i]]
    line.set_data(thisx, thisy)
    time_text.set_text(time_template % (i*dt))
    return line, time_text


ani = FuncAnimation(fig, animate, frames=np.arange(0, len(t))
                    , interval=25, blit=True,init_func=init)

ax.set_xlim(-L-0.5,L+0.5)
ax.set_ylim(-L-0.5,L+0.5)
ax.set_aspect("equal")
ax.grid()

plt.show()

#ani.save("XXX.gif", writer="pillow", fps=20)

 

 

youtu.be

 

時間変化

次に単振り子の角度、x,  yの時間変化をプロットしてみます

 

fig = plt.figure(figsize=(10,5))
ax1 = plt.subplot(111)
ax2 = ax1.twinx()

ax1.plot(sol.t, sol.y[0], label="$\\theta$")
ax2.plot(sol.t, x , color="r", linestyle="--",label="x")
ax2.plot(sol.t, y , color="g", linestyle="--",label="y")

ax1.set_yticks(np.array([-np.pi, -(1/3)*np.pi, -(1/4)*np.pi, -(1/6)*np.pi, 0,(1/6)*np.pi,(1/4)*np.pi,(1/3)*np.pi,np.pi]))
ax1.set_yticklabels(["-$\pi$","-1/3$\pi$","-1/4$\pi$","-1/6$\pi$", "0","1/6$\pi$","1/4$\pi$","1/3$\pi$","$\pi$"])
ax1.set_ylabel("$\\theta$")
ax2.set_yticks(np.array([-1,-(1/2)*np.sqrt(3),-(1/2)*np.sqrt(2),-(1/2)*np.sqrt(1),0,(1/2)*np.sqrt(1),(1/2)*np.sqrt(2),(1/2)*np.sqrt(3),1]))
ax2.set_yticklabels(["$-$1","$-\sqrt{\\frac{3}{2}}$","$-\sqrt{\\frac{2}{2}}$","$-{\\frac{1}{2}}$", "0","${\\frac{1}{2}}$","$\sqrt{\\frac{2}{2}}$","$\sqrt{\\frac{3}{2}}$","1"])
ax2.set_ylabel("X, Y")
ax1.set_xlabel("T")

ax1.set_xlim(t_span[0],t_span[-1])
ax1.set_ylim(-(1/3)*np.pi,(1/3)*np.pi)

#凡例
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, ncol=3, fontsize=12)

plt.show()

#保存
#fig.savefig("XXX.png", format="png", dpi=330)

 

f:id:RuntaScience:20200621174250p:plain

 

単振り子の軌跡

 

%matplotlib nbagg
G = 9.8  # 重力加速度(m/s^2)
L = 1.0  # 振り子の長さ(m)

# 運動方程式を解く
def eom(t, state):
    dydt = np.zeros_like(state)
    dydt[0] = state[1]
    dydt[1] = -(G/L)*np.sin(state[0])
    return dydt

#設定
t_span = [0,20] #時間 (s)
dt = 0.05# 時間の間隔 (s)
t = np.arange(t_span[0], t_span[-1], dt)# 時間の間隔のarray作成

#初期状態設定
th1 = 45.0# 角度 [deg]
w1 = 0.0 # 初期角速度 [deg/s]
state = np.radians([th1, w1])# 初期状態

#微分方程式を解く
sol = solve_ivp(eom, t_span, state, t_eval=t)
theta = sol.y[0] #y[0]=θ、y[1]=ω

x = L * np.sin(theta)       
y = -L * np.cos(theta)      

#ーーーーーーグラフーーーーーーー
fig, ax = plt.subplots()

line, = ax.plot([], [], "o-", linewidth=2) 
locus, = ax.plot([], [], 'r-', linewidth=2)
xlocus, ylocus = [], []
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, " ", transform=ax.transAxes)


def init():
    line.set_data([], [])
    time_text.set_text("")
    return line, time_text

def animate(i):
    thisx = [0, x[i]]
    thisy = [0, y[i]]
    line.set_data(thisx, thisy)
    xlocus.append(x[i])
    ylocus.append(y[i])

    locus.set_data(xlocus, ylocus)
    time_text.set_text(time_template % (i*dt))
    return line,xlocus, ylocus, time_text


ani = FuncAnimation(fig, animate, frames=np.arange(0, len(t))
                    , interval=25, blit=True,init_func=init)

ax.set_xlim(-L-0.5,L+0.5)
ax.set_ylim(-L-0.5,L+0.5)
ax.set_aspect("equal")
ax.grid()

plt.show()

#ani.save("XXX.gif", writer="pillow", fps=20)

 

youtu.be

 

複数の単振り子

次に複数の単振り子です。

 

%matplotlib nbagg
G = 9.8  # 重力加速度(m/s^2)
L = 1.0  # 振り子の長さ(m)

# 運動方程式を解く
def eom(t, state):
    dydt = np.zeros_like(state)
    dydt[0] = state[1]
    dydt[1] = -(G/L)*np.sin(state[0])
    return dydt

#設定
t_span = [0,20] #時間 (s)
dt = 0.05# 時間の間隔 (s)
t = np.arange(t_span[0], t_span[-1], dt)# 時間の間隔のarray作成

#初期状態設定

th1s = np.arange(40.0,90.0,10.0)
w1 = 0.0 # 初期角速度 [deg/s]

states=[]
for th1 in th1s:
    state = np.radians([th1, w1])# 初期状態
    states.append(state)

thetas=[]
for state in states:
    sol = solve_ivp(eom, t_span, state, t_eval=t)
    theta = sol.y[0] #y[0]=θ、y[1]=ω
    thetas.append(theta)

xs,ys=[],[]
for theta in thetas:
    x = L * np.sin(theta)  
    y = -L * np.cos(theta)  
    xs.append(x)
    ys.append(y)
    
    
#ーーーーーーグラフーーーーーーー
fig, ax = plt.subplots()

lines=[]
for i in range(len(th1s)):
    line, = ax.plot([], [], 'o-', lw=2, label=th1s[i])
    lines.append(line,)

time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, " ", transform=ax.transAxes)


def init():
    for line in lines:
        line.set_data([], [])
    time_text.set_text("")
    return line, time_text

def animate(i):
    for j in range(len(lines)):
        line = lines[j]
        thisx = [0, xs[j][i]]
        thisy = [0, ys[j][i]]
        line.set_data(thisx, thisy)

    time_text.set_text(time_template % (i*dt))
    return line, time_text


ani = FuncAnimation(fig, animate, frames=np.arange(0, len(t))
                    , interval=25, blit=True,init_func=init)

ax.set_xlim(-L-0.5,L+0.5)
ax.set_ylim(-L-0.5,L+0.5)
ax.legend(title="init $\\theta$")
ax.set_aspect("equal")
ax.grid()

plt.show()

3ani.save("XXX.gif", writer="pillow", fps=20)

 

youtu.be

 

 

次は二重振り子によるカオスをアニメーションで描写したいと思います!

 

それでは🌏