RuntaScience diary

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

Welcome to my blog

About

Matplotlibで振り子ー二重振り子 カオスを知る

こんにちは!
今日は二重振り子でカオスを見てみたいと思います!

 

単振り子に基礎的なコードを載せています。

 

 

runtascience.hatenablog.com

 

 

 

 

 二重振り子

ドキュメント

matplotlib.org

 

モジュール

まずは必要なモジュールです。

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import matplotlib.animation as animation

 

軌跡

 まずは二重振り子の軌跡を描いてみます。

 

%matplotlib nbagg
#物理量設定
G = 9.8 #重力加速度(m/s^2)
L1 = 1.0  #振り子(上)の長さ(m)
L2 = 1.0  #振り子(下)の長さ(m)
M1 = 1.0  #振り子(上)の重さ(kg)
M2 = 1.0  #振り子(下)の重さ(kg)


def derivs(state, t):

    dydx = np.zeros_like(state)
    dydx[0] = state[1]

    delta = state[2] - state[0]
    den1 = (M1+M2) * L1 - M2 * L1 * np.cos(delta) * np.cos(delta)
    dydx[1] = ((M2 * L1 * state[1] * state[1] * np.sin(delta) * np.cos(delta)
                + M2 * G * np.sin(state[2]) * np.cos(delta)
                + M2 * L2 * state[3] * state[3] * np.sin(delta)
                - (M1+M2) * G * np.sin(state[0]))
               / den1)

    dydx[2] = state[3]

    den2 = (L2/L1) * den1
    dydx[3] = ((- M2 * L2 * state[3] * state[3] * np.sin(delta) * np.cos(delta)
                + (M1+M2) * G * np.sin(state[0]) * np.cos(delta)
                - (M1+M2) * L1 * state[1] * state[1] * np.sin(delta)
                - (M1+M2) * G * np.sin(state[2]))
               / den2)

    return dydx

#初期値
dt = 0.05 #時間間隔
t = np.arange(0, 20, dt) #時間間隔分
th1 = 120.0 #振り子(上)初期角度
w1 = 0.0 #振り子(上)角速度
th2 = -10.0 #振り子(下)初期角度
w2 = 0.0 #振り子(下)角速度

state1 = np.radians([th1, w1, th2, w2]) # 初期状態

#微分方程式を解く
y01 = integrate.odeint(derivs, state1, t)


x1 = L1 * np.sin(y01[:, 0])
y1 = -L1 * np.cos(y01[:, 0])
x2 = L2 * np.sin(y01[:, 2]) + x1
y2 = -L2 * np.cos(y01[:, 2]) + y1

#--------------------グラフ-----------------------
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect("equal")
ax.grid()

line, = ax.plot([], [], "o-", lw=2)
locus, = ax.plot([], [], "r-", linewidth=2) #軌跡用
time_template = "time = %.1fs"
time_text = ax.text(0.05, 0.9, "", transform=ax.transAxes)

xlocus, ylocus = [], []
def init():
    line.set_data([], [])
    locus.set_data([], [])
    time_text.set_text("")
    return line,locus, time_text


def animate(i):
    thisx1 = [0, x1[i], x2[i]]
    thisy1 = [0, y1[i], y2[i]]
    xlocus.append(x2[i])
    ylocus.append(y2[i])

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



ani = animation.FuncAnimation(fig, animate, range(1, len(y01)),
                              interval=dt*1000, blit=True, init_func=init)
plt.show()
# ani.save("XXX.gif", writer="pillow", fps=20)

 

youtu.be

 

 

 

2個

 次にいきなり複数の振り子にいくのはなんなので、2つの振り子を描写してみます。

 

%matplotlib nbagg
#物理量設定
G = 9.8 #重力加速度(m/s^2)
L1 = 1.0  #振り子(上)の長さ(m)
L2 = 1.0  #振り子(下)の長さ(m)
M1 = 1.0  #振り子(上)の重さ(kg)
M2 = 1.0  #振り子(下)の重さ(kg)


def derivs(state, t):

    dydx = np.zeros_like(state)
    dydx[0] = state[1]

    delta = state[2] - state[0]
    den1 = (M1+M2) * L1 - M2 * L1 * np.cos(delta) * np.cos(delta)
    dydx[1] = ((M2 * L1 * state[1] * state[1] * np.sin(delta) * np.cos(delta)
                + M2 * G * np.sin(state[2]) * np.cos(delta)
                + M2 * L2 * state[3] * state[3] * np.sin(delta)
                - (M1+M2) * G * np.sin(state[0]))
               / den1)

    dydx[2] = state[3]

    den2 = (L2/L1) * den1
    dydx[3] = ((- M2 * L2 * state[3] * state[3] * np.sin(delta) * np.cos(delta)
                + (M1+M2) * G * np.sin(state[0]) * np.cos(delta)
                - (M1+M2) * L1 * state[1] * state[1] * np.sin(delta)
                - (M1+M2) * G * np.sin(state[2]))
               / den2)

    return dydx

dt = 0.05
t = np.arange(0, 20, dt)

th1 = 120.0
th12 = 120.1
w1 = 0.0
th2 = -10.0
w2 = 0.0

#初期状態
state1 = np.radians([th1, w1, th2, w2])
state2 = np.radians([th12, w1, th2, w2])

#微分方程式を解く
y01 = integrate.odeint(derivs, state1, t)
y02 = integrate.odeint(derivs, state2, t)

x1 = L1 * np.sin(y01[:, 0])
y1 = -L1 * np.cos(y01[:, 0])
x2 = L2 * np.sin(y01[:, 2]) + x1
y2 = -L2 * np.cos(y01[:, 2]) + y1

x3 = L1 * np.sin(y02[:, 0])
y3 = -L1 * np.cos(y02[:, 0])
x4 = L2 * np.sin(y02[:, 2]) + x3
y4 = -L2 * np.cos(y02[:, 2]) + y3

fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect("equal")
ax.grid()

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


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


def animate(i):
    thisx1 = [0, x1[i], x2[i]]
    thisy1 = [0, y1[i], y2[i]]
    thisx2 = [0, x3[i], x4[i]]
    thisy2 = [0, y3[i], y4[i]]

    line.set_data(thisx1, thisy1)
    line2.set_data(thisx2, thisy2)
    time_text.set_text(time_template % (i*dt))
    return line,line2, time_text



ani = animation.FuncAnimation(fig, animate, range(1, len(y01)),
                              interval=dt*1000, blit=True, init_func=init)
plt.show()
#ani.save("XXX.gif", writer="pillow", fps=20)

 

youtu.be

 

 

カオス

 最後にカオスを描写してみます。

先ほどよりも振り子数が一気に増えます。

 

%matplotlib nbagg
#物理量設定
G = 9.8 #重力加速度(m/s^2)
L1 = 1.0  #振り子(上)の長さ(m)
L2 = 1.0  #振り子(下)の長さ(m)
M1 = 1.0  #振り子(上)の重さ(kg)
M2 = 1.0  #振り子(下)の重さ(kg)


def derivs(state, t):

    dydx = np.zeros_like(state)
    dydx[0] = state[1]

    delta = state[2] - state[0]
    den1 = (M1+M2) * L1 - M2 * L1 * np.cos(delta) * np.cos(delta)
    dydx[1] = ((M2 * L1 * state[1] * state[1] * np.sin(delta) * np.cos(delta)
                + M2 * G * np.sin(state[2]) * np.cos(delta)
                + M2 * L2 * state[3] * state[3] * np.sin(delta)
                - (M1+M2) * G * np.sin(state[0]))
               / den1)

    dydx[2] = state[3]

    den2 = (L2/L1) * den1
    dydx[3] = ((- M2 * L2 * state[3] * state[3] * np.sin(delta) * np.cos(delta)
                + (M1+M2) * G * np.sin(state[0]) * np.cos(delta)
                - (M1+M2) * L1 * state[1] * state[1] * np.sin(delta)
                - (M1+M2) * G * np.sin(state[2]))
               / den2)

    return dydx

dt = 0.05
t = np.arange(0, 20, dt)


th1s = np.arange(120.0,120.01,0.0001)
w1 = 0.0
th2 = -10.0
w2 = 0.0


#初期状態
states=[]
for th1 in th1s:
    state = np.radians([th1, w1, th2, w2])
    states.append(state)
    
#微分方程式を解く
ys=[]
for state in states:
    y = integrate.odeint(derivs, state, t)
    ys.append(y)
    
x1s,x2s,y1s,y2s=[],[],[],[]
for y in ys:
    x1 = L1 * np.sin(y[:, 0])
    y1 = -L1 * np.cos(y[:, 0])
    x2 = L2 * np.sin(y[:, 2]) + x1
    y2 = -L2 * np.cos(y[:, 2]) + y1
    x1s.append(x1)
    y1s.append(y1)
    x2s.append(x2)
    y2s.append(y2)

fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(-2, 2), ylim=(-2, 2))
ax.set_aspect("equal")
ax.grid()

lines=[]
for i in range(len(th1s)):
    line, = ax.plot([], [], "o-", lw=2)
    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 lines, time_text


def animate(i):
    for j in range(len(lines)):
        line = lines[j]
        thisx1 = [0, x1s[j][i], x2s[j][i]]
        thisy1 = [0, y1s[j][i], y2s[j][i]]

        line.set_data(thisx1, thisy1)
    time_text.set_text(time_template % (i*dt))
    return lines, time_text


ani = animation.FuncAnimation(fig, animate, range(1, len(y)),
                              interval=dt*1000, blit=True, init_func=init)
plt.show()

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

 

 

youtu.be

 

 

 

それでは🌏