見出し画像

【Python】楕円経路に拘束された質点の運動シミュレーション

シミュ君

導出した経路に拘束された質点運動をシミュレーションするための計算アルゴリズムを用いた２つ目の例として楕円経路を取り上げるよ。計算に必要な経路に関する各種ベクトル（経路ベクトル・接線ベクトル・曲率ベクトル）の導出から始めるよ。

楕円経路の経路ベクトル・接線ベクトル・曲率ベクトル

楕円とは２次元平面上において、ある２点からの距離の和が等距離となる点の集合で定義される図形だね。２点をつなぐ直線をx軸、その直線の垂直方向をy軸と定義し、２点の中点を原点と定義するね。この２点は焦点と呼ばれ焦点間距離を$${R}$$、２点からの距離の和を$${L}$$（$${=L_1+L_2}$$）としたとき、上記の条件を満たす点（x, y）の集合を表す図形が楕円ってことだよ。

楕円の定義

楕円を表す方程式は

$$
a^2 = \frac{L^2}{4} \ , \ , b^2 = \frac{L^2-R^2}{4}
$$

として次のように表されるね。

中心が原点の楕円の式

$$
\frac{x^2}{a^2} + \frac{y^2}{b^2} = 1
$$

経路は媒介変数で表す必要があるね。

楕円経路の媒介変数表示

$$
\left\{\begin{matrix}x(\theta )= a\cos \theta & \cr y(\theta )=b\sin \theta& \end{matrix}\right.
$$

楕円の媒介変数表示の模式図

円経路のときと同様に、$${\theta=0}$$を始点、$${\theta=2\pi}$$を終点とした経路を経路ベクトル$${\boldsymbol{r}_{\rm path} (\theta)}$$と定義し、接線ベクトルと曲率ベクトルの導出を行うね。これらのベクトル量の計算に必要な$${\theta}$$の$${l}$$の１回微分は

$$
\frac{d l}{d\theta } = \sqrt{a^2\sin^2\theta+b^2\cos^2\theta} \ \to \ \frac{d \theta}{dl } = \frac{1}{\sqrt{a^2\sin^2\theta+b^2\cos^2\theta}}
$$

２回微分は

$$
\frac{d^2 \theta}{dl^2 } = \frac{d\theta}{dl}\,\frac{d }{d\theta}\left(\frac{1}{\sqrt{a^2\sin^2\theta+b^2\cos^2\theta}}\right)=\frac{(b^2-a^2)\sin\theta\cos\theta}{\left(a^2\sin^2\theta+b^2\cos^2\theta\right)^2}
$$

と与えられるね。楕円経路の接線ベクトルと曲率ベクトルは次のとおり得られるね

楕円経路の接線ベクトル

$$
\boldsymbol{t}_{\rm path} (\theta)= \frac{1}{\sqrt{a^2\sin^2\theta+b^2\cos^2\theta}} \,(-a\sin \theta ,b\cos\theta)
$$

楕円経路の曲率ベクトル

$$
\boldsymbol{\chi} (\theta) = \frac{-ab}{\left(a^2\sin^2\theta+b^2\cos^2\theta\right)^2}\,(b\cos\theta, a\sin\theta )
$$

ちなみに曲率ベクトルの大きさは円経路とは異なり$${\theta}$$に依存するね。

$$
\chi(\theta)=|\boldsymbol{\chi} (\theta)|= \frac{ab}{\left(a^2\sin^2\theta+b^2\cos^2\theta\right)^{3/2}}
$$

$${a>b}$$の場合、曲率の最小値と最大値はそれぞれ$${\chi(\pi/2)=b/a}$$と$${\chi(\pi/2)=a/b}$$だね。さっそく計算結果を示すよ。

計算結果

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
###############################
# 実験室グローバル変数
###############################
TIME_START = 0.0        #計算開始時刻
TIME_END   = 10.0       #計算終了時刻
DT = 0.01               #時間間隔
PLOT_FLAG = False        #グラフ描画
ANIMATION_FLAG = True  #アニメーション作成
ANIMATION_INTERVAL = 5 #アニメーション作成時の描画インターバル
OUTPUT_FLAG = False
OUTPUT_PASS = "./data_20241112/"
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])
#質量
m = 1.0
#初速度の大きさ
#v0s = [ 0.2 ]
v0s = np.arange(0.02, 0.05+0.0001, 0.0001)
###############################
# 関数定義
###############################
#相対誤差
rtol = 1e-9
#絶対誤差                   
atol = 1e-9                   
#補正ばね弾性係数
compensationK = 0.004
#補正粘性抵抗係数
compensationGamma = 0.4

# 経路の設定
class Path:
	#コンストラクタ
	def __init__(self):
		#楕円の長径と短径
		self.a = 0.7
		self.b = 1.4
		#楕円の中心位置ベクトル
		self.center = np.array([0.0, 0.0, self.b])

	#経路の位置ベクトルを指定する媒介変数関数
	def position( self, theta ):
		x = self.a * np.cos( theta ) + self.center[0]
		y = 0
		z = self.b * np.sin( theta ) + self.center[2]
		return np.array([x, y, z])

	#接線ベクトルを指定する媒介変数関数
	def tangent( self, theta ):
		K = 1.0 / np.sqrt( (self.a * np.sin(theta))**2 + (self.b * np.cos(theta))**2 )
		x = -K * self.a * np.sin(theta)
		y = 0
		z =  K * self.b * np.cos(theta)
		return np.array([x, y, z])
	#曲率ベクトルを指定する媒介変数関数
	def curvature( self, theta ):
		K = -self.a * self.b * 1.0 / ( (self.a * np.sin(theta))**2 + (self.b * np.cos(theta))**2 )**2
		x = K * self.b * np.cos(theta)
		y = 0
		z = K * self.a * np.sin(theta)
		return np.array([x, y, z])
	#媒介変数の取得
	def getTheta( self, r ):
		#相対位置ベクトル
		bar_r = r - self.center

		#楕円方程式の右辺
		c = np.sqrt( (bar_r[0] / self.a)**2 + (bar_r[2] / self.b )**2 )
		#楕円方程式のパラメータの変換
		_a = self.a * c
		_b = self.a * c
		sinTheta = bar_r[2] / _b

		if (sinTheta > 0): theta = np.arccos(bar_r[0] / _a)
		else: theta = 2 * np.pi - np.arccos(bar_r[0] / _b)

		return theta


#経路拘束力の微分方程式
def equations(t, X, path, compensationK, compensationGamma):
	vec_r =  np.array([ X[0], X[1], X[2] ])
	vec_v =  np.array([ X[3], X[4], X[5] ])

	#媒介変数の取得
	theta = path.getTheta ( vec_r )

	#媒介変数に対する位置ベクトル、接線ベクトル、曲率ベクトルの計算
	position = path.position( theta )
	tangent = path.tangent( theta )
	curvature = path.curvature( theta )

	#経路中心の位置ベクトル・速度ベクトル・加速度ベクトル
	vec_r_box = np.copy( path.center )
	vec_v_box = np.array([0.0, 0.0, 0.0])
	vec_a_box = np.array([0.0, 0.0, 0.0])

	#相対速度
	bar_v = vec_v - vec_v_box

	#微係数dl/dtとd^2l/dt^2を計算
	dl_dt = np.dot( bar_v, tangent )
	d2l_dt2 = np.dot( vec_g, tangent ) - np.dot( vec_a_box, tangent) 

	#加速度を計算
	vec_a = vec_a_box + d2l_dt2 * tangent + dl_dt * dl_dt * curvature

	#ズレ位置ベクトル
	bar_r = vec_r - position
	#補正ばね弾性力
	fk = -compensationK * bar_r
	#法線ベクトル
	n1 = curvature / np.sqrt(np.sum( curvature**2 ))
	n2 = np.cross(n1, tangent )
	#補正粘性抵抗力
	fgamma1 = -compensationGamma * np.dot( bar_v, n1 ) * n1 
	fgamma2 = -compensationGamma * np.dot( bar_v, n2 ) * n2
	#経路補正力を加える
	vec_a += (fk + fgamma1 + fgamma2)/m * np.abs(m*vec_a)

	return [vec_v[0], vec_v[1], vec_v[2], vec_a[0], vec_a[1], vec_a[2]]

# 初速度依存性を可視化
def calculation(vec_r0, vec_v0, path, compensationK, compensationGamma):

	solution = integrate.solve_ivp(
		equations,                     #常微分方程式
		method = 'RK45',                #計算方法
		args = (path, compensationK, compensationGamma),              #引数
		t_span=( times[0], times[-1]),  #積分区間
		rtol = rtol,                   #相対誤差
		atol = atol,                   #絶対誤差
		y0 = np.concatenate([vec_r0, vec_v0]), #初期値
		t_eval = times                  #結果を取得する時刻
	)
	#計算結果
	return solution.y[ 0 ], solution.y[ 1 ], solution.y[ 2 ], solution.y[ 3 ], solution.y[ 4 ], solution.y[ 5 ]

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

	#経路の生成
	path = Path()
	#描画用リスト
	xls = []
	yls = []
	zls = []
	vxls = []
	vyls = []
	vzls = []
	#初速度 による違い 
	for i in range(len(v0s)):
		print( i )

		#初期位置と速度
		vec_r0 = path.position( np.pi/2 )
		vec_v0 = np.array([ v0s[i], 0, 0])

		x, y, z, vx, vy, vz = calculation( vec_r0, vec_v0, path, compensationK, compensationGamma )
		xls.append( x )
		yls.append( y )
		zls.append( z )
		vxls.append( vx )
		vyls.append( vy )
		vzls.append( vz )

		vec_r = np.array([ x[-1], y[-1], z[-1] ])
		vec_v = np.array([ vx[-1],vy[-1],vz[-1] ])

		E_start = 1.0/2.0 * m * np.sum( vec_v0**2 ) - m * np.dot( vec_g, vec_r0 )
		E_end =   1.0/2.0 * m * np.sum( vec_v**2 )  - m * np.dot( vec_g, vec_r )
		DeltaE = np.abs(E_end - E_start)  
		print( E_start, E_end, " |ΔE/E(0)| =", DeltaE / E_start )

	if( OUTPUT_FLAG ):
		for i in range(len(times)):
			f_path = OUTPUT_PASS + str( i ) + '.txt'
			f = open(f_path, 'w')

			for j in range( len(xls) ):
				f.write( str(round( xls[ j ][ i ], 4)) + "\t" +  str(round( yls[ j ][ i ], 4)) + "\t" + str(round( zls[ j ][ i ], 4)) + "\t" +  str(round( vxls[ j ][ i ], 4)) + "\t" + str(round( vyls[ j ][ i ], 4)) + "\t" +  str(round( vzls[ j ][ i ], 4)) + "\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"$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([0, TIME_END])
		plt.ylim([0, 2])

		#グラフの描画
		for i in range(len(v0s)):
			c = colorsys.hls_to_rgb( i/len(v0s) , 0.5 , 1)
			plt.plot(times, zls[i], linestyle='solid', linewidth = 4)#, color = c)

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

	if( ANIMATION_FLAG ):

		#楕円座標点配列の生成
		thetas = np.linspace(0, 2*np.pi, 200)
		cirsx = path.a * np.cos( thetas )
		cirsy = path.b * np.sin( thetas ) + path.center[ 2 ]

		#################################################################
		# アニメーションの設定
		#################################################################
		fig = plt.figure(figsize=(8*4/5, 10*4/5)) #画像サイズ
		plt.subplots_adjust(left = 0.07, right = 0.98, bottom = 0.07, 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']

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

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

		#軌跡リスト
		r_trajectorys = [ 0.0 ] * len( v0s )

		#アニメーションの各グラフを生成
		for an in range( AN ):
			t = DT * an

			if( an % ANIMATION_INTERVAL != 0): continue

			#軌跡の描画
			img = plt.plot(cirsx, cirsy, color = "black", linestyle='dotted', linewidth = 1 )

			for i in range( len( v0s ) ):
				c = colorsys.hls_to_rgb( i/len(v0s) , 0.5 , 1)
				if( i == 0 ): 
					img += plt.plot([xls[i][an]], [zls[i][an]], color = c, marker = 'o', markersize = 20)
				else:
					img += plt.plot([xls[i][an]], [zls[i][an]], color = c, marker = 'o', markersize = 20)

			
			#時刻をグラフに追加
			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([-1.2, 1.2])
		plt.ylim([-0.1, 2.9])
		#アニメーションの生成
		ani = animation.ArtistAnimation(fig, imgs, interval=50)
		#アニメーションの保存
		#ani.save("output.html", writer=animation.HTMLWriter())
		ani.save("output2.gif", writer="imagemagick")
		ani.save("output.mp4", writer="ffmpeg", dpi=360)
		#グラフの表示
		#plt.show()

いいなと思ったら応援しよう！