【Python】任意多重結合振子によるカオス運動

上記の２重振子運動における張力の計算アルゴリズムは、手動で導出しました。今回は、以前に解説した任意多重結合による運動の計算アルゴリズムに従って、張力に関する連立方程式もPythonで計算できるように改造するね。

計算結果

次のアニメーションは固定された支点からつながった三角形型のおもりと、そこからぶら下がる２個のおもり、合計５個のおもりと６本の接続の系だね。接続数が６本だから、張力に関して６本の連立方程式となるのだけれども、 np.linalg.solve関数を用いて数値的に計算するよ。それなりに計算できているね。

Pythonプログラムソース

Pythonプログラムソースを用意しました。よかったら試してみてください。

#################################################################
TITLE = "任意の結合振子によるカオス運動"
#################################################################
import os
import numpy as np
import scipy.integrate as integrate
import matplotlib.pyplot as plt
import matplotlib.animation as animation
import colorsys
INF = float('inf')
###############################
# 実験室グローバル変数
###############################
TIME_START = 0.0        #計算開始時刻
TIME_END   = 1.0       #計算終了時刻
DT = 0.01               #時間間隔
PLOT_FLAG = True        #グラフ描画
ANIMATION_FLAG = True   #アニメーション作成
ANIMATION_INTERVAL = 3 #アニメーション作成時の描画インターバル
OUTPUT_FLAG = False
OUTPUT_PASS = "./data_20241004/"
if( OUTPUT_FLAG and os.path.isdir(OUTPUT_PASS) == False) : os.mkdir(OUTPUT_PASS)
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])

#支点の振幅と角振動数
l = 0.0
omega_box = 0 * np.pi

#質量
ml = [ INF, 1.0, 1.0, 1.0, 1.0, 1.0  ]
NumberOfFix = 1

#結合リスト
connections = []
connections.append( [ 0, 1 ] )
connections.append( [ 1, 2 ] )
connections.append( [ 2, 3 ] )
connections.append( [ 1, 3 ] )
connections.append( [ 2, 4 ] )
connections.append( [ 3, 5 ] )

#初期位置と速度
vec_r = [
	np.array([ 0, 0, 1.0]),
	np.array([ -1.0/2.0, 0, 1.0 + np.sqrt(3.0)/2.0 ]),
	np.array([  1.0/2.0, 0, 1.0 + np.sqrt(3.0)/2.0 ]),
	np.array([ -1.0/2.0, 0, 1.0 + np.sqrt(3.0)/2.0 + 1.0]),
	np.array([  1.0/2.0, 0, 1.0 + np.sqrt(3.0)/2.0 + 1.0])
] 
vec_v = [ 
	np.array([ 0, 0, 0]), 
	np.array([ 0, 0, 0]), 
	np.array([ 0, 0, 0]), 
	np.array([ 0, 0, 0]), 
	np.array([ 0.01, 0, 0])
] 

#ひもの長さ（初期値）
L0 = []
for n in range(len(connections)):
	i = connections[ n ][ 0 ]
	j = connections[ n ][ 1 ]
	if( i == 0 ):
		ri = np.array([ 0, 0, 0 ])
	else:
		ri = vec_r[ i - NumberOfFix ]
	rj = vec_r[ j - NumberOfFix ]
	rij = ri - rj
	L0.append( np.sqrt( np.dot( rij, rij ) ) )

###############################
# 関数定義
###############################
#相対誤差
rtol = 1e-6
#絶対誤差                   
atol = 1e-6                   
#補正ばね弾性係数
compensationK = 0.06
#補正粘性抵抗係数
compensationGamma = 0.5

#支点の位置と速度
def box(t):
	vec_r_box = np.array([ l * np.sin( omega_box * t ) , 0, 0])
	vec_v_box = np.array([ l * omega_box * np.cos( omega_box * t ), 0, 0])
	vec_a_box = np.array([ -l * omega_box**2 * np.sin( omega_box * t ), 0, 0])
	return vec_r_box, vec_v_box, vec_a_box

#接続番号を取得
def getConectionNumber( i, j ):
	if( i == j ): return -1
	for m in range( len(connections) ):
		connection = connections[ m ]
		if( i in connection and j in connection ): return m
	return -1

#張力＋重力の微分方程式
def equations(t, X, compensationK, compensationGamma):
	vec_r = [0] * len(ml)
	vec_v = [0] * len(ml)
	vec_a = [0] * len(ml)
	j = 0
	for i in range(len(ml)):
		if(ml[ i ] == INF):
			#支点の位置ベクトル・速度ベクトル・加速度ベクトル
			vec_r[ i ], vec_v[ i ], vec_a[ i ] = box( t )
		else:
			vec_a[ i ]	= vec_g
			vec_r[ i ] = np.array([ X[ 6 * j + 0 ], X[ 6 * j + 1 ], X[ 6 * j + 2 ] ])
			vec_v[ i ] = np.array([ X[ 6 * j + 3 ], X[ 6 * j + 4 ], X[ 6 * j + 5 ] ])
			j += 1

	N = len(connections)
	vec_rij = [ 0 ] * N
	vec_vij = [ 0 ] * N
	vec_aij = [ 0 ] * N
	vec_nij = [ 0 ] * N
	Lij = [ 0 ] * N
	vij = [ 0 ] * N
	for n in range(N):
		i = connections[ n ][ 0 ]
		j = connections[ n ][ 1 ]
		#相対位置ベクトル
		vec_rij[ n ] = vec_r[ i ] - vec_r[ j ]
		#ひもの長さ
		Lij[ n ] =  np.sqrt( np.dot(vec_rij[n], vec_rij[n]) ) 
		#相対単位ベクトル
		vec_nij[ n ] = vec_rij[ n ] / Lij[ n ]
		#相対速度ベクトル
		vec_vij[ n ] = vec_v[ i ] - vec_v[ j ]
		#相対速度の大きさ
		vij[ n ] = np.sqrt( np.dot(vec_vij[n], vec_vij[n]))
		#相対加速度ベクトル
		vec_aij[ n ] = vec_a[ i ] - vec_a[ j ]

	#a係数（右辺）
	A = np.array( [ 0.0 ] * N )
	for n in range( N ):
		A[ n ] = -1.0 / Lij[ n ] * ( vij[ n ]**2 + np.dot(vec_aij[ n ], vec_rij[ n ]) ) 
	#f成分（左辺）
	F = np.zeros(( N, N ))
	for n in range( N ):
		i = connections[ n ][ 0 ]
		j = connections[ n ][ 1 ]
		for k in range( len( ml ) ):
			nki = getConectionNumber( k, i )
			nkj = getConectionNumber( k, j )
			if( nki > -1):
				vec_n = vec_nij[ nki ]
				if( i < k ) : vec_n = -vec_n
				F[ n ][ nki ] += np.dot( vec_n, vec_nij[ n ]) / ml[ i ]
			if( nkj > -1):
				vec_n = vec_nij[ nkj ]
				if( j < k ) : vec_n = -vec_n
				F[ n ][ nkj ] -= np.dot( vec_n, vec_nij[ n ]) / ml[ j ]

	#print( np.linalg.det(F) )
	# 求解
	vec_f = np.linalg.solve(F, A)
	
	#質点に加わる力
	vec_F = [0.0] * len(ml)
	for i in range(len(ml)):
		if(ml[ i ] == INF): continue
		vec_F[ i ] = ml[ i ] * vec_g
		for j in range( len( ml )):
			nn = getConectionNumber( i, j )
			if( nn > -1):
				f = vec_f[ nn ] * vec_nij[ nn ]
				if( i < j ): f = -f
				vec_F[ i ] += f

	#補正力
	vec_fc = [0.0] * N
	for n in range(N):
		i = connections[ n ][ 0 ]
		j = connections[ n ][ 1 ]
		vec_fc[ n ] = (compensationK * ( Lij[ n ] - L0[ n ] ) + compensationGamma * np.dot(vec_vij[ n ], vec_nij[ n ]) ) * vec_nij[ n ]
		if( i > j ): vec_fc[ n ] = -vec_fc[ n ] 
		
	#補正因子
	r = [ 0.0 ] * len(ml)
	for i in range(len(ml)):
		if(ml[ i ] == INF): continue
		for j in range(len(ml)):
			nn = getConectionNumber( i, j )
			if( nn > -1 ):
				r[ i ] += np.sqrt( np.dot( vec_f[ nn ],  vec_f[ nn ]) )  

	#ニュートンの運動方程式
	vec_a = [0.0] * len(ml)
	n = 0
	for i in range(len(ml)):
		if(ml[ i ] == INF): continue
		vec_a[ i ] = (vec_F[ i ] + vec_fc[ n ] * r[ i ] )/ ml[ i ]
		n += 1

	H = []
	for i in range(len(ml)):
		if(ml[ i ] == INF): continue
		H.append( vec_v[ i ][ 0 ])
		H.append( vec_v[ i ][ 1 ])
		H.append( vec_v[ i ][ 2 ])
		H.append( vec_a[ i ][ 0 ])
		H.append( vec_a[ i ][ 1 ])
		H.append( vec_a[ i ][ 2 ])
	return H

def calculation(vec_r, vec_v, compensationK, compensationGamma):
	init = []
	for i in range(len(vec_r)):
		init.append( vec_r[i] )
		init.append( vec_v[i] )

	solution = integrate.solve_ivp(
		equations,                     #常微分方程式
		method = 'RK45',                #計算方法
		args = (compensationK, compensationGamma),  #引数
		t_span=( times[0], times[-1]),  #積分区間
		rtol = rtol,                   #相対誤差
		atol = atol,                   #絶対誤差
		y0 = np.concatenate(init),      #初期値
		t_eval = times                  #結果を取得する時刻
	)
	#計算結果
	vec_r_t = []
	vec_v_t = []
	for i in range(len(vec_r)):
		vec_r_t.append( np.array( [solution.y[ 6 * i + 0 ], solution.y[ 6 * i + 1 ], solution.y[ 6 * i + 2 ]] ) )
		vec_v_t.append( np.array( [solution.y[ 6 * i + 3 ], solution.y[ 6 * i + 4 ], solution.y[ 6 * i + 5 ]] ) )

	return vec_r_t, vec_v_t

##########################################################
if __name__ == "__main__":

	vec_r_t, vec_v_t = calculation( vec_r, vec_v, compensationK, compensationGamma)

	DeltaLs = []
	ts = []
	for tn in range(len(times)):
		t = tn * DT
		vec_r_box, vec_v_box, vec_a_box = box( t )

		Lij = [0.0] * len( connections )
		for cn in range(len( connections )) :
			i = connections[ cn ][ 0 ]
			j = connections[ cn ][ 1 ]
			if(ml[ i ] != INF):
				ri = np.array( [ vec_r_t[ i - NumberOfFix ][ 0 ][ tn ], vec_r_t[ i - NumberOfFix ][ 1 ][ tn ], vec_r_t[ i - NumberOfFix ][ 2 ][ tn ] ])
			else:
				ri = vec_r_box
			if(ml[ j ] != INF):
				rj = np.array( [ vec_r_t[ j - NumberOfFix  ][ 0 ][ tn ], vec_r_t[ j - NumberOfFix ][ 1 ][ tn ], vec_r_t[ j - NumberOfFix ][ 2 ][ tn ] ])
			else:
				rj = vec_r_box			
			
			rij = ri - rj
			Lij[ cn ] =  np.sqrt( np.dot( rij, rij ) )

		DeltaL = 0
		for cn in range(len( connections )) :
			DeltaL += np.abs( Lij[ cn ] - L0[ cn ] ) / len( connections )
		
		if(tn%10==9): 
			print( DeltaL )
			ts.append( t )
			DeltaLs.append( DeltaL )

	if( OUTPUT_FLAG ):
		for tn in range(len(times)):
			t = tn * DT
			vec_r_box, vec_v_box, vec_a_box = box( t )
			path = OUTPUT_PASS + str( tn ) + '.txt'
			f = open(path, 'w')
			for j in range( len(ml) ):
				if( ml[j] == INF ):
					for xyz in range(3):
						f.write( str( round( vec_r_box[ xyz ], 2)) + "\t")
					for xyz in range(3):
						f.write( str( round( vec_v_box[ xyz ], 2)) + "\t")
					f.write("\n")
				else:
					for xyz in range(3):
						f.write( str( round( vec_r_t[ j - NumberOfFix][ 0 ][ tn ], 2)) + "\t")
					for xyz in range(3):
						f.write( str( round( vec_v_t[ j - NumberOfFix][ 0 ][ tn ], 2)) + "\t")
					f.write("\n")
			f.close()

	if( PLOT_FLAG ):
		#################################################################
		# 計算誤差の計算
		#################################################################
		fig = plt.figure( figsize=(12, 8) ) #画像サイズ
		plt.subplots_adjust(left = 0.085, right = 0.97, bottom = 0.10, top = 0.95) #グラフサイズ
		plt.rcParams['font.family'] = 'Times New Roman' #フォント
		plt.rcParams['font.size'] = 16 #フォントサイズ
		plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
		#カラーリストの取得
		colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
		#グラフタイトルと軸ラベルの設定 
		plt.title( TITLE, fontname="Yu Gothic", fontsize=20, fontweight=500 )
		plt.xlabel( r"$t \,[{\rm s}]$", fontname="Yu Gothic", fontsize=20, fontweight=500 )
		plt.ylabel( r"$|\Delta L| \,[{\rm m}]$" , fontname="Yu Gothic", fontsize=20, fontweight=500 )
		#罫線の描画
		plt.grid( which = "major", axis = "x", alpha = 0.7, linewidth = 1 )
		plt.grid( which = "major", axis = "y", alpha = 0.7, linewidth = 1 )

		#描画範囲を設定
		plt.xlim([0, TIME_END])
		#plt.ylim([0, 2.0])
		plt.yscale("log")  
		#グラフの描画
		plt.plot(ts, DeltaLs, linestyle='solid', linewidth = 4)

		#グラフの表示
		plt.show()


	if( ANIMATION_FLAG ):

		#################################################################
		# アニメーションの設定
		#################################################################
		fig = plt.figure( figsize=(8, 6) ) #画像サイズ
		plt.subplots_adjust(left = 0.095, right = 0.97, bottom = 0.10, top = 0.92) #グラフサイズ
		plt.rcParams['font.family'] = 'Times New Roman' #フォント
		plt.rcParams['font.size'] = 16 #フォントサイズ
		plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
		#カラーリストの取得
		colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

		#アニメーション作成用
		imgs = []

		#アニメーション分割数
		AN = len( times )

		#アニメーションの各グラフを生成
		for an in range( AN ):
			t = DT * an
			vec_r_box, vec_v_box, vec_a_box = box( t )
			if( an % ANIMATION_INTERVAL != 0): continue

			c = colorsys.hls_to_rgb( 0.2 , 0.5 , 1)

			img  = plt.plot([ vec_r_box[0]], [vec_r_box[2]], marker = 's', markersize = 5, color = colors[ 0 ])
			img += plt.plot([ vec_r_box[ 0 ], vec_r_t[ 0 ][ 0 ][ an ] ], [vec_r_box[ 2 ],  vec_r_t[ 0 ][ 2 ][ an ]], linestyle='solid', linewidth = 2, color = colors[ 0 ])

			NumberOfFix = 1
			for cn in range(len( connections ) - 1):
				i = connections[ cn + 1 ][ 0 ]
				j = connections[ cn + 1 ][ 1 ]
				img += plt.plot([ vec_r_t[ i - NumberOfFix ][ 0 ][ an ], vec_r_t[  j - NumberOfFix ][ 0 ][ an ] ], [vec_r_t[ i - NumberOfFix ][ 2 ][ an ], vec_r_t[ j - NumberOfFix ][ 2 ][ an ]], linestyle='solid', linewidth = 2, color = colors[ 0 ])

			NumberOfFix = 1
			for i in range(len(ml)):
				if( i == 0): continue 
				img += plt.plot([ vec_r_t[ i - NumberOfFix ][ 0 ][ an ] ], [vec_r_t[ i - NumberOfFix ][ 2 ][ an ]],  marker = 'o', markersize = 20, color = colors[i])

			#時刻をグラフに追加
			time = plt.text(13.5, 0.35, "t = " +  str("{:.2f}".format(t)), fontsize=18)
			img.append( time )
			
			imgs.append( img )

		#グラフタイトルと軸ラベルの設定 
		plt.title( TITLE, fontname="Yu Gothic", fontsize=20, fontweight=500 )
		plt.xlabel( r"$x \,[{\rm m}]$", fontname="Yu Gothic", fontsize=20, fontweight=500 )
		plt.ylabel( r"$z \,[{\rm m}]$" , fontname="Yu Gothic", fontsize=20, fontweight=500 )
		#罫線の描画
		plt.grid( which = "major", axis = "x", alpha = 0.7, linewidth = 1 )
		plt.grid( which = "major", axis = "y", alpha = 0.7, linewidth = 1 )

		#描画範囲を設定
		plt.xlim([-4.2, 4.2])
		plt.ylim([-3.2, 3.2])
		#アニメーションの生成
		ani = animation.ArtistAnimation(fig, imgs, interval=50)
		#アニメーションの保存
		ani.save("output.html", writer=animation.HTMLWriter())
		#ani.save("output.gif", writer="imagemagick")
		#ani.save("output.mp4", writer="ffmpeg", dpi=360)
		#グラフの表示
		plt.show()