リアプノフ指数の値が負から正へ
解決したいこと
4重振り子において、4つ目の質点の初期角度を10^-10だけずらした2つの系の角度の差をとり、リアプノフ指数でカオスを判定しました。そうすると、最終的には、正に収束していくのですが、最初の方はわずかですが負の値をとりました。これは、結果として正しい結果なのでしょうか、それとも私のコードがどこかで間違えてしまっているのでしょうか？
また、これが正しい結果だとした場合、リアプノフ指数が負から正になるというのは、どういったことを意味しているのでしょうか？
該当するソースコード
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import matplotlib.animation as animation
import pandas as pd
df_param = pd.read_excel('python N重振り子.xlsx',header=1,sheet_name=None, engine='openpyxl')
df_param['Sheet1'].index = ['質量m(kg)', '長さL(m)','初速度0(m/s)']
df_param['Sheet1'] = df_param['Sheet1'].drop('個数',axis=1)
m1 = np.array([df_param['Sheet1'].iat[0, i] for i in range(len(df_param['Sheet1'].columns))]).astype(np.float64).round(2)
L1 = np.array([df_param['Sheet1'].iat[1, i] for i in range(len(df_param['Sheet1'].columns))]).astype(np.float64).round(2)
theta1 = np.array([90,90,90,90])
x1 = np.array([(i) for i in theta1]).astype(np.float64)
v1 = np.array([df_param['Sheet1'].iat[2, i] for i in range(len(df_param['Sheet1'].columns))]).astype(np.float64).round(2)
df_param['Sheet2'].index = ['質量m(kg)', '長さL(m)','初速度0(m/s)']
df_param['Sheet2'] = df_param['Sheet2'].drop('個数',axis=1)
m2 = np.array([df_param['Sheet2'].iat[0, i] for i in range(len(df_param['Sheet2'].columns))]).astype(np.float64).round(2)
L2 = np.array([df_param['Sheet2'].iat[1, i] for i in range(len(df_param['Sheet2'].columns))]).astype(np.float64).round(2)
theta2 = np.array([90,90,90,90+10**-10])
x2 = np.array([(i) for i in theta2]).astype(np.float64)
v2 = np.array([df_param['Sheet2'].iat[2, i] for i in range(len(df_param['Sheet2'].columns))]).astype(np.float64).round(2)
ms=[m1,m2]
Ls=[L1,L2]
thetas=[theta1,theta2]
xs=[x1,x2]
vs=[v1,v2]
#初期条件のコピー
xs0 = xs.copy()
print('m=',ms,type(ms),len(ms))
print('L=',Ls,type(Ls),len(Ls))
print('theta=',thetas,type(thetas),len(thetas))
print('x=',xs,type(xs),len(xs))
print('v=',vs,type(vs),len(vs))
print('列数=',len(df_param['Sheet1'].columns))
#初期位置の確認
n = len(df_param['Sheet1'].columns)
x_ini_cor1 = np.zeros(n,dtype=np.float64)
y_ini_cor1 = np.zeros(n,dtype=np.float64)
def ini_cor_func1(j):
if j == 0:
x_ini_cor1[j] = L1[j]*np.sin(x1[j])
y_ini_cor1[j] = -L1[j]*np.cos(x1[j])
else:
x_ini_cor1[j] = L1[j]*np.sin(x1[j]) + x_ini_cor1[j-1]
y_ini_cor1[j] = -L1[j]*np.cos(x1[j])+ y_ini_cor1[j-1]
return x_ini_cor1[j], y_ini_cor1[j]
for j in range(n):
x_ini_cor1[j] , y_ini_cor1[j] = ini_cor_func1(j)
x_ini_cor2 = np.zeros(n,dtype=np.float64)
y_ini_cor2 = np.zeros(n,dtype=np.float64)
def ini_cor_func2(j):
if j == 0:
x_ini_cor2[j] = L2[j]*np.sin(x2[j])
y_ini_cor2[j] = -L2[j]*np.cos(x2[j])
else:
x_ini_cor2[j] = L2[j]*np.sin(x2[j]) + x_ini_cor2[j-1]
y_ini_cor2[j] = -L2[j]*np.cos(x2[j])+ y_ini_cor2[j-1]
return x_ini_cor2[j], y_ini_cor2[j]
for j in range(n):
x_ini_cor2[j] , y_ini_cor2[j] = ini_cor_func2(j)
xplot1_ = np.insert(x_ini_cor1,0,0)
yplot1_ = np.insert(y_ini_cor1,0,0)
xplot2_ = np.insert(x_ini_cor2,0,0)
yplot2_ = np.insert(y_ini_cor2,0,0)
init = 0
end = 100
dt = 0.05
h = dt
loop = int(end/h)
n = len(df_param['Sheet2'].columns)
g = 9.8
# initial state
t = init
tpoints = np.arange(init, end , h)
xpoints1 = []
vpoints1 = []
# A = np.zeros((n,n),dtype=np.float64)
# B = np.zeros((n,n),dtype=np.float64)
E1 = -np.ones_like(x1)
def N_func1(t, x1, v1):
A1 = np.zeros((n,n),dtype=np.float64)
B1= np.zeros((n,n),dtype=np.float64)
for i in range(n):
for j in range(n):
for k in range(max(i,j),n):
A1[i][j] += m1[k]
B1[i][j] += m1[k]
if i == j:
A1[i][j] *= L1[j]
B1[i][j] *= g * np.sin(x1[i])
else:
A1[i][j] *= L1[j]*np.cos(x1[i]-x1[j])
B1[i][j] *= L1[j]*v1[j]**2*np.sin(x1[i]-x1[j])
#逆行列の計算
inv_A1 = np.linalg.inv(A1)
#inv_A*Bを計算
inv_A1_B1 = np.dot(inv_A1, B1)
F1 = np.dot(inv_A1_B1, E1)
return F1
xpoints1 = []
vpoints1 = []
#配列要素数の定義
j11 = np.zeros_like(v1)
k11 = np.zeros_like(x1)
j12 = np.zeros_like(v1)
k12 = np.zeros_like(x1)
j13 = np.zeros_like(v1)
k13 = np.zeros_like(x1)
j14 = np.zeros_like(v1)
k14 = np.zeros_like(x1)
def RK1(t,x1,v1):
vt1 = v1.copy()
xt1 = x1.copy()
xpoints1.append(xt1)
vpoints1.append(vt1)
j11 = N_func1(t, x1, v1) * h
k11 = v1 * h
j12 = N_func1(t + h / 2, x1 + k11 / 2, v1 + j11 / 2) * h
k12 = (v1 + j11/ 2)* h
j13 = N_func1(t + h / 2, x1 + k12 / 2, v1 + j12 / 2) * h
k13 = (v1 + j12/ 2)* h
j14 = N_func1(t + h, x1 + k13, v1 + j13)*h
k14 = (v1 + j13)* h
v1 += (j11 + 2*j12 + 2*j13 + j14)/6
x1 += (k11 + 2*k12 + 2*k13 + k14)/6
return x1,v1,xpoints1 ,vpoints1
E2 = -np.ones_like(x2)
def N_func2(t, x2, v2):
A2 = np.zeros((n,n),dtype=np.float64)
B2 = np.zeros((n,n),dtype=np.float64)
for i in range(n):
for j in range(n):
for k in range(max(i,j),n):
A2[i][j] += m2[k]
B2[i][j] += m2[k]
if i == j:
A2[i][j] *= L2[j]
B2[i][j] *= g * np.sin(x2[i])
else:
A2[i][j] *= L2[j]*np.cos(x2[i]-x2[j])
B2[i][j] *= L2[j]*v2[j]**2*np.sin(x2[i]-x2[j])
#逆行列の計算
inv_A2 = np.linalg.inv(A2)
#inv_A*Bを計算
inv_A2_B2 = np.dot(inv_A2, B2)
F2 = np.dot(inv_A2_B2, E2)
return F2
xpoints2 = []
vpoints2 = []
#配列要素数の定義
j21 = np.zeros_like(v2)
k21 = np.zeros_like(x2)
j22 = np.zeros_like(v2)
k22 = np.zeros_like(x2)
j23 = np.zeros_like(v2)
k23 = np.zeros_like(x2)
j24 = np.zeros_like(v2)
k24 = np.zeros_like(x2)
def RK2(t,x2,v2):
vt2 = v2.copy()
xt2 = x2.copy()
xpoints2.append(xt2)
vpoints2.append(vt2)
j21 = N_func2(t, x2, v2) * h
k21 = v2 * h
j22 = N_func2(t + h / 2, x2 + k21 / 2, v2 + j21 / 2) * h
k22 = (v2 + j21/ 2)* h
j23 = N_func2(t + h / 2, x2 + k22 / 2, v2 + j22 / 2) * h
k23 = (v2 + j22/ 2)* h
j24 = N_func2(t + h, x2 + k23, v2 + j23)*h
k24 = (v2 + j23)* h
v2 += (j21 + 2*j22 + 2*j23 + j24)/6
x2 += (k21 + 2*k22 + 2*k23 + k24)/6
return x2,v2,xpoints2 ,vpoints2
# from ipykernel import kernelapp as app
for t in range(len(tpoints)):
# for t in range(2):
x1, v1, xpoints1, vpoints1= RK1(t,x1,v1)
for t in range(len(tpoints)):
# for t in range(2):
x2, v2, xpoints2, vpoints2= RK2(t,x2,v2)
xpoints1 = np.array(xpoints1)
vpoints1 = np.array(vpoints1)
xpoints2 = np.array(xpoints2)
vpoints2 = np.array(vpoints2)
theta_tips1 = []
for x in xpoints1:
theta_tips1.append(x[-1]) # 配列の末尾へは-1でアクセスできます
theta_tips2 = []
for x in xpoints2:
theta_tips2.append(x[-1]) # 配列の末尾へは-1でアクセスできます
theta_tips1=np.array(theta_tips1)
theta_tips2=np.array(theta_tips2)
fig = plt.figure()
ax = fig.add_subplot(111)
#ax.set_xlim(0,100)
#ax.set_ylim()
#plt.plot(tpoints,np.log10(theta_tips2-theta_tips1))
plt.show()
def lyapunov_exponent():
lyapunov_exponent = []
d_0 = abs((theta_tips2[0]-theta_tips1[0])) #初期値段階での差異の絶対値
for i in range(int(len(theta_tips2))):
d_i = abs(theta_tips2[i] - theta_tips1[i]) # i期での差異の絶対値
d_rate = (d_i/d_0) # i期と初期値のずれ幅の比率
ln_d_rate = np.log10(d_rate)
lmd = (1/(i+1)) *ln_d_rate
lyapunov_exponent.append(lmd)
return lyapunov_exponent
print(lyapunov_exponent())
lyapunov_exponent = lyapunov_exponent()
ax.set_title('Lyapunov Exponent')
plt.plot(lyapunov_exponent, label="Lyapunov Exponent")
plt.show()
0
