【Science】 Matplotlibで振り子ー二重振り子 カオスを知る
こんにちは!
今日は二重振り子でカオスを見てみたいと思います!
単振り子に基礎的なコードを載せています。
二重振り子
ドキュメント
モジュール
まずは必要なモジュールです。
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)
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)
カオス
最後にカオスを描写してみます。
先ほどよりも振り子数が一気に増えます。
%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)
それでは🌏