@syotalosposted at 2021-11-13

N重振り子のカオスを判定

Q&A

解決したいこと

N重振り子のカオスを判定するために、最後の質点のみ初期値を10^-10だけずらした二つの系の角度の差を縦軸に、時間を横軸にプロットしました。角度差はlog10（）を用いているため、その初期値は-10になるはずです。
しかし、微妙に違った値をとってしまいます。これは私のコードに問題があるのでしょうか。それともこの結果が正しい結果なのでしょうか。詳しい方、ご教授願います。。

該当するソースコード

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([60. , 60,60])
x1 = np.array([np.radians(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([60. ,60, 60+10**-10])
x2 = np.array([np.radians(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(-15,5)

plt.plot(tpoints,np.log10(theta_tips2-theta_tips1))
plt.show()

@bizzpaper posted at 2021-11-13

弧度法に変換しているからではないでしょうか？

import numpy as np

a = 60.
b = 60. + 1e-10

print(b - a)
# 1.000017846308765e-10

x = np.radians(a)
y = np.radians(b)

print(y - x)
# 1.7454926393156711e-12

また、コーディングが少し目に余ったのでアドバイスをしておきます。

グローバル変数が多すぎる

グロバール変数は便利ですが、むやみやたらに使用するとバグの温床になります。
変数のスコープは使用する範囲をでないようできる限り、狭くするのが望ましいです（「変数スコープ」で検索）。

配列演算

numpyにはブロードキャストという機能があります。
これは、同じ大きさの配列であれば要素ごとに演算してくれる大変便利な機能です。
この機能を使えば、以下の二つのコードの結果は同じになります。

import numpy as np

a = np.full(1, 10, float)
b = np.full(2, 10, float)
c = np.zeros_like(a)

for i in range(n):
	c[i] = a[i] * b[i]

import numpy as np

a = np.full(1, 10, float)
b = np.full(2, 10, float)

c = a * b

さらに、これには処理の高速化という効用もあります。
pythonは数値計算によく使われるCなど他の言語と比べ、処理速度が格段に遅いです。
そのため、こういった機能は活用していくのが望ましいです。

関数化

シート１とシート２にはほとんど同じ処理を行っています。
それをコピーしてべた書きで対応しいるとシートがさらに増えたとき詰みます。
したがって、処理を俯瞰し、切り分けできるところは関数化していくことが望ましいです（「プログラム再利用性」で検索）。

最初は時間がかかるかもしれませんが、処理の切り出しを練習しておいて損はないです。

また、本気でやりたいなら、関数の上位のクラスを勉強することを強くお勧めします。
私はクラスを勉強してから、格段にプログラムを書く力が向上したと感じています。

ファイル分割

結果の解析ではトライ＆エラーが比較的多く行われることが想定されます。
そのため、解析する度に特に変更のないメイン計算を行ってしまうため、待ち時間が非常に多く発生します。
以下の例のように目的に合わせてファイルを分割することで、この待ち時間を軽減できます。

目的	説明
メイン計算	メイン計算を行う。
解析1	解析を行う。
プロット1	解析１の結果のグラフを作成する。
解析2	解析を行う。処理時間が短ければ、一緒にプロットしてもよい。
プロット2	解析２の結果のグラフを作成する。

それぞれの結果はCSVなどで一時的に保存し、読み込んで使用するとよいでしょう。

お節介だったかもしれませんが、@syotalos さんは根気強くこのテーマに取り組んでおられるようだったので、お力になってあげたいと思いアドバイスさせていただきました！

0Like

Are you sure you want to delete the question?