古典力学 – シミュ君日記

2021年9月19日

【第８話】サイクロイド曲線に経路が拘束された粒子の運動【ルンゲクッタで行こう！⑧】

############################################################################
#　[008]サイクロイド曲線に経路が拘束された粒子の運動
############################################################################
import os
import numpy as np
import random 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
random.seed(0)

#図全体
fig = plt.figure(figsize=(15, 8))
#全体設定
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']
#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
#時間区間
t_min = 0
t_max = 20.0
#時間区間数
NT =   20000
skip =    20

#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#重力加速度
g = np.array([0.0, 0.0, -10.0])

#################################################################
## 経路の設定
#################################################################
class Path:
	#コンストラクタ
	def __init__(self):
		#始点と終点のx座標
		self.x_start = -1
		self.x_end = 1
		#サイクロイド曲線パラメータr
		self.r = ( self.x_end - self.x_start ) / ( 2 * np.pi )
		#始点と終点のz座標
		self.path_z0 = 2 * self.r ; #最下点が0となるように
		#サイクロイド曲線始点の位置ベクトル
		self.start = np.array([self.x_start, 0.0, self.path_z0])

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

	#接線ベクトルを指定する媒介変数関数
	def tangent( self, theta ):
		A = 1.0 / np.sqrt( 2.0 )
		x = A * np.sqrt( 1.0 - np.cos(theta))
		y = 0
		z = - A * np.sqrt( 1.0 + np.cos(theta))
		if( np.sin(theta) < 0 ): z = - z
		return np.array([x, y, z])

	#曲率ベクトルを指定する媒介変数関数
	def curvature( self, theta ):
		A = - 1.0 / ( 4.0 * self.r )
		x = - A * np.sqrt( ( 1.0 + np.cos(theta) ) / ( 1.0 - np.cos(theta) ) )
		y = 0
		z = - A
		if( np.abs(x) == np.inf ):
			x = 0
			z = 0
		if( np.sin(theta) < 0 ): x = - x

		return np.array([x, y, z])
	#媒介変数の取得
	def getTheta( self, r, v ):
		#球体の相対位置ベクトル
		bar_r = r -self.start
		if( np.abs( v[0] ) > 0.1 ):
			_r = np.abs(bar_r[2]) / 2.0 * (  1.0 + (v[2]/v[0])**2 )
		else :
			_r = self.r

		if( bar_r[0] < _r * np.pi ) :
			theta = np.arccos( 1.0 + bar_r[2]/_r )
		else:
			theta = 2.0 * np.pi - np.arccos( 1.0 + bar_r[2]/_r )


		return theta

#経路の生成
path = Path()

#補正パラメータ
compensationK = 100.0;     #補正ばね弾性係数
compensationGamma = 1.0; #補正粘性抵抗係数
compensationFactor = 1;  #補正因子

#粒子に加わる力
def F (self, t, r_t, v_t ):
	#外場（重力）の計算
	fe = m * g

	#経路自体の運動
	v0 = np.array([0.0, 0.0, 0.0])
	a0 = np.array([0.0, 0.0, 0.0])

	#相対速度
	bar_v = v_t - v0

	#媒介変数の取得
	theta = path.getTheta ( r_t, bar_v )

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

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

	#加速度を計算
	a = a0 + d2l_dt2 * tangent + dl_dt * dl_dt * curvature
	a_length = np.sqrt(np.sum( a**2 ))

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

	return m * a

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		#現時刻の位置と速度
		self.r = r0.copy()
		self.v = v0.copy()
		#１ステップ前の時刻の位置と速度
		self._r = r0.copy()
		self._v = v0.copy()
		#差分
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		#時間間隔
		self.dt = dt
		#内部
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F(self, t, r_t, v_t ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始
#################################################################
#初期条件
particles =[]
for i in range(20):
	#球体の初期位置
	theta0 = np.pi
	#初速度
	x0 = path.r * ( theta0 - np.sin(theta0) ) + path.start[0]
	y0 = 0
	z0 = - path.r * ( 1.0 - np.cos(theta0)) + path.start[2]
	#初期条件
	r0 = np.array([ x0, y0, z0 ])
	v0 = np.array([ 0.170 * (i+1), 0.0, 0.0 ])

	particles.append(  RK4( m, dt, r0, v0 ) )

####################################
## 描画用関数
####################################
#アニメーション作成用
ims=[]

#円座標点配列の生成
thetas = np.linspace(0, 2.0*np.pi, 200)
p = path.position( thetas )
cirsx = p[0]
cirsy = p[2]

def plot( t ):
	#各コマを描画
	img = plt.plot(cirsx, cirsy, color = "black", linestyle='dotted', linewidth = 1 )
	i = 0
	for particle in particles:
		img += plt.plot([particle.r[0]], [particle.r[2]], color = colors[i], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
		i += 1
		if (i>9): i = 0

	time = plt.text(0.8, 0.717, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )
	E = 1.0/2.0 * m * np.sum( particles[0].v**2 ) - m * np.dot( g, particles[0].r )
	print( "t = {:.4f}".format(t)," E =",  E )

####################################
## 各時刻における運動方程式の計算
####################################
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t)

	#ルンゲ・クッタ法による時間発展
	for particle in particles:
		particle.timeEvolution( t )
	#位置と速度を更新
	for particle in particles:
		particle.r += particle.dr
		particle.v += particle.dv

#################################################################
## 分子運動アニメーション
#################################################################
plt.title( u"サイクロイド曲線経路に拘束された球体の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.1, 1.1])
plt.ylim([-0.1, 0.7])
plt.subplots_adjust(left = 0.12, right = 0.97, bottom = 0.10, top = 0.90)

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

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

2021年9月17日2021年9月19日

【第７話】経路が拘束された粒子の運動【ルンゲクッタで行こう！⑦】

############################################################################
#　[007]経路に拘束された剛体球の運動
############################################################################
import os
import numpy as np
import random 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
random.seed(0)

#図全体
fig = plt.figure(figsize=(8, 8))
#全体設定
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']
#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
#時間区間
t_min = 0
t_max = 10.0
#時間区間数
NT =   1000
skip =   10
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#重力加速度
g = np.array([0.0, 0.0, -10.0])

#################################################################
## 経路の設定
#################################################################
class Path:
	#コンストラクタ
	def __init__(self):
		#円の半径
		self.R = 1
		#円の中心位置ベクトル
		self.center = np.array([0.0, 0.0, self.R])

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

	#接線ベクトルを指定する媒介変数関数
	def tangent( self, theta ):
		x = -np.sin(theta)
		y = 0
		z = np.cos(theta)
		return np.array([x, y, z])
	#曲率ベクトルを指定する媒介変数関数
	def curvature( self, theta ):
		x = -np.cos(theta) / self.R
		y = 0
		z = -np.sin(theta) / self.R
		return np.array([x, y, z])
	#媒介変数の取得
	def getTheta( self, r ):
		#相対位置ベクトル
		bar_r = r - self.center
		bar_r_length = np.sqrt(np.sum( bar_r**2 ))
		sinTheta = bar_r[2] / bar_r_length
		if( sinTheta > 0 ):
			theta = np.arccos( bar_r[0] / bar_r_length )
		else:
			theta = 2.0 * np.pi - np.arccos( bar_r[0] / bar_r_length )

		return theta

#経路の生成
path = Path()

#補正パラメータ
compensationK = 1.0;     #補正ばね弾性係数
compensationGamma = 1.0; #補正粘性抵抗係数
compensationFactor = 0;  #補正因子

#粒子に加わる力
def F (self, t, r_t, v_t ):
	#外場（重力）の計算
	fe = m * g

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

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

	#経路自体の運動
	v0 = np.array([0.0, 0.0, 0.0])
	a0 = np.array([0.0, 0.0, 0.0])

	#相対速度
	bar_v = v_t - v0

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

	#加速度を計算
	a = a0 + d2l_dt2 * tangent + dl_dt * dl_dt * curvature
	a_length = np.sqrt(np.sum( a**2 ))

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

	return m * a

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		#現時刻の位置と速度
		self.r = r0.copy()
		self.v = v0.copy()
		#１ステップ前の時刻の位置と速度
		self._r = r0.copy()
		self._v = v0.copy()
		#差分
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		#時間間隔
		self.dt = dt
		#内部
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F(self, t, r_t, v_t ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始
#################################################################
#初期条件
particles =[]
for i in range(20):
	particles.append(  RK4( m, dt, np.array([ 0.0, 0.0, 2.0 ]), np.array([ 0.001 * (i+1), 0.0, 0.0 ] ) ) )

####################################
## 描画用関数
####################################
#アニメーション作成用
ims=[]

#円座標点配列の生成
thetas = np.linspace(0, 2*np.pi, 200)
cirsx = np.cos( thetas )
cirsy = np.sin( thetas ) + 1.0

def plot( t ):
	#各コマを描画
	img = plt.plot(cirsx, cirsy, color = "black", linestyle='dotted', linewidth = 1 )
	i = 0
	for particle in particles:
		img += plt.plot([particle.r[0]], [particle.r[2]], color = colors[i], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
		i += 1
		if (i>9): i = 0

	time = plt.text(0.8, 2.18, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )
	E = 1.0/2.0 * m * np.sum( particles[0].v**2 ) - m * np.dot( g, particles[0].r )
	print( "t = {:.4f}".format(t)," E =",  E )

####################################
## 各時刻における運動方程式の計算
####################################
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t)

	#ルンゲ・クッタ法による時間発展
	for particle in particles:
		particle.timeEvolution( t )
	#位置と速度を更新
	for particle in particles:
		particle.r += particle.dr
		particle.v += particle.dv

#################################################################
## 分子運動アニメーション
#################################################################
plt.title( u"円経路に拘束された球体の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.1, 1.1])
plt.ylim([-0.1, 2.1])
plt.subplots_adjust(left = 0.12, right = 0.97, bottom = 0.10, top = 0.90)

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)

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

2021年9月6日2021年9月6日

【第６話】剛体球同士の衝突力計算アルゴリズム【ルンゲクッタで行こう！⑥】

############################################################################
#　箱の中のN個の自由粒子（古典力学）
############################################################################
import os
import numpy as np
import random 
import matplotlib.pyplot as plt
import matplotlib.animation as animation
random.seed(0)

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
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']
#################################################################
## 物理系の設定
#################################################################
#質量
m_min = 1.0
m_max = 10.0
#粒子の半径
radius_min = 0.002
radius_max = 0.02
#粒子数
NM = 100
#箱の長さ
L = 1.0
#粒子に加わる力
fex = np.array([0.0, 0.0, 0.0])
#初期位置と速度
v_max = 2.0

#粒子に加わる力
def F (self, t, r_t, v_t ):
	return self.fex + self.fc - self.m * self.beta * v_t

#時間区間
t_min = 0
t_max = 5.0
#時間区間数
NT =    500
skip =    1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#
lw = 0.01
######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		#現時刻の位置と速度
		self.r = r0.copy()
		self.v = v0.copy()
		#１ステップ前の時刻の位置と速度
		self._r = r0.copy()
		self._v = v0.copy()
		#差分
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		#時間間隔
		self.dt = dt
		#外力
		self.fex = np.array([0.0, 0.0, 0.0])
		#衝突力
		self.fc = np.array([0.0, 0.0, 0.0])
		#粘性係数
		self.beta = 0
		#半径
		self.radius = 0
		#内部
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F(self, t, r_t, v_t ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 時間発展の計算
#################################################################

#壁との衝突判定（戻値：法線ベクトル）
def CollisionWallDetection( particle ):
	#衝突リスト
	collisionlist = []
	
	if(particle.r[0] + particle.radius > L):
		collisionlist.append({ "type": 1, "normal" : np.array([-1.0, 0.0, 0.0]) })
	if(particle.r[0] - particle.radius < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([1.0, 0.0, 0.0]) })
	if(particle.r[1] + particle.radius > L):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, -1.0, 0.0]) })
	if(particle.r[1] - particle.radius < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, 1.0, 0.0]) })
	if(particle.r[2] + particle.radius > L):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, 0.0, -1.0]) })
	if(particle.r[2] - particle.radius < 0):
		collisionlist.append({ "type": 1, "normal" : np.array([0.0, 0.0, 1.0]) })

	return collisionlist

#粒子の運動エネルギー
def T():
	_T = 0
	for nm in range(NM):
		_T += 1.0/2.0 *  particles[nm].m * np.sum( particles[nm].v**2 )
	return  _T
#粒子の位置エネルギー
def U():
	_R = 0
	for nm in range(NM):
		_R += - np.dot(particles[nm].fex, particles[nm].r)
	return _R

#系の時間発展の計算
def timeEvolution( t, timeScale ):

	#複数衝突
	cylinderCollisionList = []

	for nm in range(NM):
		#時間スケールの指定
		particles[nm].dt = dt / timeScale

		#ルンゲ・クッタ法による時間発展
		particles[nm].fc = np.array([0.0, 0.0, 0.0])
		particles[nm].timeEvolution( t )
		#位置と速度を更新
		particles[nm].r += particles[nm].dr
		particles[nm].v += particles[nm].dv
		#壁との衝突リスト
		CollisionWallList = CollisionWallDetection( particles[nm] )
		#衝突あり
		if( len(CollisionWallList) > 0 ):

			#２つ以上の壁との衝突がある場合　→　高精度計算モードへ
			if( len(CollisionWallList) > 1 ):
				return False

			#以下、衝突が１箇所
			normal = CollisionWallList[ 0 ]["normal"]

			#いったん巻き戻しておく
			particles[nm].r -= particles[nm].dr
			particles[nm].v -= particles[nm].dv

			#衝突力
			if( CollisionWallList[ 0 ]["type"] == 1 ): #壁との衝突
				#print(nm)
				#相対速度
				barV = - particles[nm].v
				particles[nm].fc = (2.0 * particles[nm].m * np.dot(barV, normal) * normal ) / particles[nm].dt - np.dot(particles[nm].fex, normal) * normal
				#位置と速度を更新
				particles[nm].timeEvolution( t )
				particles[nm].r += particles[nm].dr
				particles[nm].v += particles[nm].dv


	#粒子同士の衝突判定
	for nm1 in range(NM):
		for nm2 in range(NM):
			if( nm1 >= nm2 ): continue
			#２点間ベクトル
			r12 = particles[nm2].r - particles[nm1].r
			r12_length = np.sqrt(np.sum( r12**2 ))

			if( r12_length < particles[nm1].radius + particles[nm2].radius ):
				#print( t, nm1, nm2 ) 
				#すでに衝突力が与えられている場合　→　高精度計算モードへ
				if( (particles[nm1].fc[0] != 0 or particles[nm1].fc[1] != 0 or particles[nm1].fc[2] != 0) 
				 or ( particles[nm2].fc[0] != 0 or particles[nm2].fc[1] != 0 or particles[nm2].fc[2] != 0 )):
					#時間発展：失敗
					return False

				n12 = r12 / r12_length
				v12 = particles[nm2].v - particles[nm1].v
				m1 = particles[nm1].m
				m2 = particles[nm2].m
				barf =  particles[nm2].fex / m2 - particles[nm1].fex / m1 
				f21 = 2.0 * m1 * m2 / ( m1 + m2 ) * np.dot(v12, n12) * n12 / particles[nm1].dt
				f21 += m1 * m2 / ( m1 + m2 ) * np.dot(barf, n12) * n12
				f12 = -f21

				#衝突力の設定
				particles[nm1].fc = f21
				particles[nm2].fc = f12
				#いったん巻き戻しておく
				particles[nm1].r -= particles[nm1].dr
				particles[nm1].v -= particles[nm1].dv
				particles[nm2].r -= particles[nm2].dr
				particles[nm2].v -= particles[nm2].dv
				#時間発展を計算
				particles[nm1].timeEvolution( t )
				particles[nm2].timeEvolution( t )
				#位置と速度を更新
				particles[nm1].r += particles[nm1].dr
				particles[nm1].v += particles[nm1].dv
				particles[nm2].r += particles[nm2].dr
				particles[nm2].v += particles[nm2].dv

				#２点間ベクトル
				r12 = particles[nm2].r - particles[nm1].r
				r12_length = np.sqrt(np.sum( r12**2 ))
				#計算後の２粒子が重なっている場合は失敗
				if( r12_length < particles[nm1].radius + particles[nm2].radius ):
					#時間発展：失敗
					return False

	#計算後、壁とあるいは粒子同士が重なっている場合は失敗
	for nm1 in range(NM):
		#壁との衝突リスト
		CollisionWallList = CollisionWallDetection( particles[nm1] )
		#衝突あり
		if( len(CollisionWallList) > 0 ):
			#時間発展：失敗
			return False
		for nm2 in range(NM):
			if( nm1 >= nm2 ): continue
			#２点間ベクトル
			r12 = particles[nm2].r - particles[nm1].r
			r12_length = np.sqrt(np.sum( r12**2 ))
			if( r12_length < particles[nm1].radius + particles[nm2].radius ):
				#時間発展：失敗
				return False


	#時間発展：正常
	return True


#高精度計算モード
def highPrecisionMode(t, level):
	print ( "レベル" , level, "t=", round(t, level+3 ))

	timeScale = 10**level
	for nm in range(NM):
		particles[nm].r = particles[nm]._r.copy()
		particles[nm].v = particles[nm]._v.copy()

	for _n in range(0, 10):
		t_level = t + dt / timeScale * _n

		result_level = timeEvolution( t_level, timeScale )

		if( result_level == False ): 
			#もっと高精度が必要
			level += 1
			highPrecisionMode(t_level, level)
			level -= 1
			
		else:
			#時間発展前の位置と速度を保持
			for nm in range(NM):
				particles[nm]._r = particles[nm].r.copy()
				particles[nm]._v = particles[nm].v.copy()		


#################################################################
## 計算開始
#################################################################
#粒子
particles = []

for nm in range(NM):

	radius = radius_min + ( radius_max - radius_min) * random.random()
	m = m_min + ( m_max - m_min) * random.random()

	theta_v = 2.0 * np.pi * random.random()
	phi_v = 0 * 2.0 * np.pi * random.random()
	_v_max = v_max * random.random()
	v0 = np.array([
		_v_max * np.sin( theta_v ) * np.cos( phi_v ), 
		_v_max * np.sin( theta_v ) * np.sin( phi_v ), 
		_v_max * np.cos( theta_v )
	])

	for i in range(1000):
		r0 = np.array([
			radius + (L - 2.0 * radius) * random.random(), 
			L / 2.0, 
			radius + (L - 2.0 * radius) * random.random()
		])
		flag = True
		for p in particles:
			#２点間ベクトル
			r12 = r0 - p.r
			r12_length =  np.sqrt(np.sum( r12**2 ))
			if( r12_length < radius + p.radius ):
				flag = False
				break
		if(flag == False): continue

		particles.append( RK4( m, dt, r0, v0 ))
		particles[ len( particles ) - 1 ].radius = radius
		break


#粒子外力の設定
for nm in range(NM):
	particles[nm].fex = fex

####################################
## 描画用関数
####################################
#アニメーション作成用
ims=[]
#ピストン位置の配列
zs = []

def plot( t ):
	#各コマを描画
	for nm in range(NM):
		markersize = particles[nm].radius * 1000
		colorRG = 1.0 - particles[nm].m / m_max
		if( nm == 0 ): img  = plt.plot([particles[nm].r[0]], [particles[nm].r[2]], color = (colorRG, colorRG, 1.0), marker = 'o', markersize = markersize, linestyle='solid', linewidth = 0 )
		else: img += plt.plot([particles[nm].r[0]], [particles[nm].r[2]], color = (colorRG, colorRG, 1.0), marker = 'o', markersize = markersize, linestyle='solid', linewidth = 0 )
	time = plt.text(0.9, 1.12, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )
	E = T() + U()
	print( "t = {:.4f}".format(t)," E =",  E )


####################################
## 各時刻における運動方程式の計算
####################################
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t)

	#時間発展前の位置と速度を保持
	for nm in range(NM):
		particles[nm]._r = particles[nm].r.copy()
		particles[nm]._v = particles[nm].v.copy()		

	level = 0

	#時間発展の計算
	result = timeEvolution( t, 1 )
	
	#複数衝突を検知した場合
	if( result == False ): 
		level += 1
		highPrecisionMode(t, level )



#################################################################
## 分子運動アニメーション
#################################################################
plt.title( u"箱の中の粒子の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-0.1, 1.1])
plt.ylim([-0.1, 1.1])


#箱の概形
plt.vlines([0-lw], -2*lw, L + 2*lw, "black", linestyles='solid', linewidth = 10 )
plt.vlines([L+lw], -2*lw, L + 2*lw, "black", linestyles='solid', linewidth = 10 )
plt.hlines([-lw], 0, L, "black", linestyles='solid', linewidth = 10 )
plt.hlines([L+lw], 0, L, "black", linestyles='solid', linewidth = 10 )

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output100.mp4", writer="ffmpeg", dpi=300)

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

2021年8月30日2021年9月14日

【第５話】シリンダー内の粒子の運動計算アルゴリズム【ルンゲクッタで行こう！⑤】

############################################################################
#　シリンダーの中の１粒子運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
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']

#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
M = 10.0
radius = 0.05
#箱の長さ
L = 2.0
#重力加速度
g = np.array([0.0, 0.0, -10.0])

def F ( t, m, r_t, v_t, fc ):
	return m * g + fc

#初期位置と速度
r0 = np.array([0.0 , 0.0, 0.0])
v0 = np.array([0.5, 0.0, 2.5])
R0 = np.array([0.0 , 0.0, L/2])
V0 = np.array([0.0, 0.0, 0.0])
#時間区間
t_min = 0
t_max = 10.0
#時間区間数
NT =  2000
skip =  1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

#衝突判定（戻値：法線ベクトル）
def CollisionDetection( r, R ):
	if(r[0] + radius > L/2):
		return 1, np.array([-1.0, 0.0, 0.0])
	if(r[0] - radius < -L/2):
		return 1, np.array([1.0, 0.0, 0.0])
	if(r[1] + radius > L/2):
		return 1, np.array([0.0, -1.0, 0.0])
	if(r[1] - radius < -L/2):
		return 1, np.array([0.0, 1.0, 0.0])
	if(r[2] + radius> R[2]):
		return 2, np.array([0.0, 0.0, -1.0])
	if(r[2] - radius < -L/2):
		return 1, np.array([0.0, 0.0, 1.0])
	#衝突なし
	return False, np.array([0.0, 0.0, 0.0])
	
######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, m, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.m = m
		self.r = r0.copy()
		self.v = v0.copy()
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		self.dt = dt
		self.fc = np.array([0.0, 0.0, 0.0])
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F( t, self.m, r_t, v_t, self.fc ) / self.m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 


def T(v, V):
	return 1.0/2.0 * m * np.sum( v**2 ) + 1.0/2.0 * M * np.sum( V**2 )
def U(r, R):
	return  - m * np.dot(g, r) - M * np.dot(g, R)
#################################################################
## 計算開始
#インスタンス
rk4_m = RK4(m, dt, r0, v0)
rk4_M = RK4(M, dt, R0, V0)

#アニメーション作成用
ims=[]
def plot( t, r, v, R, V ):
	E = T(v, V) + U(r, R)
	print( "{:.2f}".format(t), E )
	#各コマを描画
	img  = plt.plot([r[0]], [r[2]], colors[0], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img += plt.plot([-L/2, L/2], [R[2],R[2]], colors[1], linestyle='solid', linewidth = 10 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4_m.r, rk4_m.v, rk4_M.r, rk4_M.v )

	#ルンゲ・クッタ法による時間発展
	rk4_m.fc = np.array([0.0, 0.0, 0.0])
	rk4_M.fc = np.array([0.0, 0.0, 0.0])
	rk4_m.timeEvolution( t )
	rk4_M.timeEvolution( t )
	#位置と速度を更新
	rk4_m.r += rk4_m.dr
	rk4_m.v += rk4_m.dv
	rk4_M.r += rk4_M.dr
	rk4_M.v += rk4_M.dv

	CollisionFlag , n = CollisionDetection( rk4_m.r, rk4_M.r )
	if( CollisionFlag ):
		#いったん巻き戻して
		rk4_m.r -= rk4_m.dr
		rk4_m.v -= rk4_m.dv
		rk4_M.r -= rk4_M.dr
		rk4_M.v -= rk4_M.dv
		#衝突力
		if( CollisionFlag == 1 ):
			barV = - rk4_m.v
			rk4_m.fc = (2.0 * m * np.dot(barV, n) * n ) / dt - m * np.dot(g, n) * n
			rk4_M.fc = np.array([0.0, 0.0, 0.0])
		elif( CollisionFlag == 2 ):
			barV = rk4_M.v - rk4_m.v
			rk4_m.fc = 2.0 * m * M / ( m + M ) * np.dot(barV, n) * n / dt
			rk4_M.fc = - rk4_m.fc

		rk4_m.timeEvolution( t )
		rk4_M.timeEvolution( t )
		#位置と速度を更新
		rk4_m.r += rk4_m.dr
		rk4_m.v += rk4_m.dv
		rk4_M.r += rk4_M.dr
		rk4_M.v += rk4_M.dv

plt.title( u"シリンダーの中の粒子の運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#箱の概形
plt.vlines([-L/2], -L/2, L/2, "black", linestyles='solid', linewidth = 10 )
plt.vlines([L/2], -L/2, L/2, "black", linestyles='solid', linewidth = 10 )
plt.hlines([-L/2], -L/2, L/2, "black", linestyles='solid', linewidth = 10 )

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

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

2021年8月27日2021年8月30日

【第４話】強制振り子の計算アルゴリズム【ルンゲクッタで行こう！④】

############################################################################
#　強制振り子運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
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']

#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0

#支点の初期位置
box_r0 = np.array([0.0, 0.0, 0.0])

#初期位置と速度
r0 = np.array([0.0, 0.0, -1.0])
v0 = np.array([0.0, 0.0, 0.0])

#位置ベクトル
L0 = r0 - box_r0
#ひもの長さ
L0_abs = np.sqrt( np.sum( L0**2 ) )

#強制振動の大きさと角振動数
box_l = 0.2
box_omega = 0.5 * (2.0 * np.pi) 

#支点の運動関数
def boxFunction( t ):
	box_rx = box_r0[0] + box_l * np.sin( box_omega * t)
	box_vx = box_l * box_omega * np.cos( box_omega * t)
	box_ax = - box_l * box_omega**2 * np.sin( box_omega * t)
	box_r = np.array([box_rx, 0.0, 0.0])
	box_v = np.array([box_vx, 0.0, 0.0])
	box_a = np.array([box_ax, 0.0, 0.0])
	return box_r, box_v, box_a


#重力加速度
g = np.array([0.0, 0.0, -10.0])
#補正パラメータ
compensationK = 0.0
compensationGamma = 0.0

#張力ベクトル
def S ( t, r_t, v_t):
	#支点の位置ベクトル
	box_r, box_v, box_a = boxFunction(t)	
	#おもり位置ベクトル
	L_t = r_t - box_r
	#おもり速度ベクトル
	v_bar = v_t - box_v
	#ひもの長さ
	L_abs = np.sqrt( np.sum(L_t**2) )

	return -m / L_abs ** 2 * ( np.sum(v_bar**2) - np.dot(box_a, L_t) + np.dot(g, L_t) ) * L_t

def F ( t, r_t, v_t ):
	#支点の位置ベクトル
	box_r, box_v, box_a = boxFunction(t)

	#おもり位置ベクトル
	L_t = r_t - box_r
	#おもり速度ベクトル
	v_bar = v_t - box_v
	#ひもの長さ
	L_abs = np.sqrt( np.sum( L_t**2 ) )
	#おもり速度の大きさ
	v_bar_abs = np.sqrt( np.sum( v_bar**2 ) )
	#単位方向ベクトル
	L_hat = L_t.copy() / L_abs
	#張力ベクトル
	S_ = S( t, r_t, v_t ) 
	#張力の大きさ
	S_coefficient = np.sqrt( np.sum( S_**2 ) )

	#補正力ベクトル
	compensationK = S_coefficient 
	compensationGamma = v_bar_abs * 10
	fc1 = - compensationK * ( L_abs - L0_abs ) * L_hat
	fc2 = - compensationGamma * np.dot(v_bar, L_hat) * L_hat

	#合力ベクトル
	return  S_ + m * g + fc1 + fc2


#時間区間
t_min = 0
t_max = 40.0
#時間区間数
NT = 4000
skip =  2
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.r = r0.copy()
		self.v = v0.copy()
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		self.dt = dt
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F( t, r_t, v_t ) / m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)

		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始

#インスタンス
rk4 = RK4(dt, r0, v0)

#アニメーション作成用
ims=[]
def plot( t, r ):
	box_r, box_v, box_a = boxFunction(t)
	print( "{:.2f}".format(t),  np.sqrt((r[0]-box_r[0])**2 + (r[2]-box_r[2])**2) - 1  )

	#各コマを描画
	img  = plt.plot([box_r[0]], [box_r[2]], colors[0], marker = 's', markersize = 10, linestyle='solid', linewidth = 0 )
	img  += plt.plot([r[0]], [r[2]], colors[1], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img  += plt.plot([box_r[0], r[0]], [box_r[2], r[2]], colors[1], linestyle='solid', linewidth = 1 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4.r )
	#ルンゲ・クッタ法による時間発展
	rk4.timeEvolution( t )
	#位置と速度を更新
	rk4.r += rk4.dr
	rk4.v += rk4.dv

plt.title( u"強制単振り子運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#アニメーションの生成
ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
ani.save("output.mp4", writer="ffmpeg", dpi=300)

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

2021年8月24日2021年9月3日

【第３話】拘束条件の改良！【ルンゲクッタで行こう！③】

20210.08.27 アルゴリズムを改良しました。破綻しません。

############################################################################
#　単振り子運動シミュレーション（補正力あり）
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
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']

#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
#ひもの長さ
r_abs0 = 1.0
#重力加速度
g = np.array([0.0, 0.0, -10.0])
#補正パラメータ
compensationK = 100.0
compensationGamma = 0.0

#張力ベクトル
def S ( t, r_t, v_t ):
	r_abs = np.sqrt( np.sum(r_t**2) )
	return -m / r_abs ** 2 * ( np.sum(v_t**2) + np.dot(g, r_t) ) * r_t

def F ( t, r_t, v_t ):
	r = np.sqrt( np.sum( r_t**2 ) )
	v = np.sqrt( np.sum( v_t**2 ) )
	r_hat = r_t.copy() / r
	v_hat = v_t.copy() / v
	fc  = - compensationK * ( r - r_abs0 ) * r_hat
	fc += - compensationGamma * np.dot(v, r_hat) * r_hat
	return S( t, r_t, v_t ) + m * g + fc


#初期位置と速度
r0 = np.array([0.0 , 0.0, -r_abs0])
v0 = np.array([6.0, 0.0, 0.0])
#時間区間
t_min = 0
t_max = 10.0
#時間区間数
NT = 1000
skip =  1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.r = r0.copy()
		self.v = v0.copy()
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		self.dt = dt
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F( t, r_t, v_t ) / m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)	
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始

#インスタンス
rk4 = RK4(dt, r0, v0)

#アニメーション作成用
ims=[]
def plot( t, r ):
	print( "{:.2f}".format(t),  np.sqrt(r[0]**2 + r[2]**2) - 1  )
	#各コマを描画
	img  = plt.plot([0], [0], colors[0], marker = 'o', markersize = 10, linestyle='solid', linewidth = 0 )
	img  += plt.plot([r[0]], [r[2]], colors[1], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img  += plt.plot([0, r[0]], [0, r[2]], colors[1], linestyle='solid', linewidth = 1 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4.r )
	#ルンゲ・クッタ法による時間発展
	rk4.timeEvolution( t )
	#位置と速度を更新
	rk4.r += rk4.dr
	rk4.v += rk4.dv



plt.title( u"単振り子運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#アニメーションの生成
#ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)

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

2021年8月22日2021年9月14日

【第２話】単振り子の作り方【ルンゲクッタで行こう！②】

20210.08.27 アルゴリズムを改良しました。破綻しません。

############################################################################
#　単振り子運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

#図全体
fig = plt.figure(figsize=(10, 10))
#全体設定
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']

#################################################################
## 物理系の設定
#################################################################
#質量
m = 1.0
#ひもの長さ
r = 1.0
#重力加速度
g = np.array([0.0, 0.0, -10.0])

#張力ベクトル
def S ( t, r_t, v_t ):
	vv = np.sum(v_t**2)
	bb = np.dot(g, r_t)
	r_abs = np.sqrt( np.sum(r_t**2) )
	S_abs =  -m * (vv + bb) / (r_abs ** 2)
	return S_abs * r_t 

def F ( t, r_t, v_t ):
	return S( t, r_t, v_t ) + m * g

#初期位置と速度
r0 = np.array([0.0 , 0.0, -r])
v0 = np.array([6.0, 0.0, 0.0])
#時間区間
t_min = 0
t_max = 10.0
#時間区間数
NT = 1000
skip = 1
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
	#コンストラクタ
	def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
		self.r = r0.copy()
		self.v = v0.copy()
		self.dr = np.array([0.0, 0.0, 0.0])
		self.dv = np.array([0.0, 0.0, 0.0])
		self.dt = dt
		self.__v1 = np.array([0.0, 0.0, 0.0])
		self.__v2 = np.array([0.0, 0.0, 0.0])
		self.__v3 = np.array([0.0, 0.0, 0.0])
		self.__v4 = np.array([0.0, 0.0, 0.0])
		self.__a1 = np.array([0.0, 0.0, 0.0])
		self.__a2 = np.array([0.0, 0.0, 0.0])
		self.__a3 = np.array([0.0, 0.0, 0.0])
		self.__a4 = np.array([0.0, 0.0, 0.0])

	#速度を与えるメソッド
	def V(self, t, r, v):
		return v.copy()
	##########################################################
	#加速度を与えるメソッド
	def A(self, t, r_t, v_t):
		return F( t, r_t, v_t ) / m
	##########################################################

	#時間発展を計算するメソッド
	def timeEvolution(self, t):
		#１段目
		self.__v1 = self.V(
			t, 
			self.r, 
			self.v
		)
		self.__a1 = self.A(
			t, 
			self.r, 
			self.v
		)
		#２段目
		self.__v2 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0
		)
		self.__a2 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v1 * self.dt / 2.0, 
			self.v + self.__a1 * self.dt / 2.0			
		)
		#３段目
		self.__v3 = self.V(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0					
		)
		self.__a3 = self.A(
			t + self.dt / 2.0, 
			self.r + self.__v2 * self.dt / 2.0, 
			self.v + self.__a2 * self.dt / 2.0			
		)
		#４段目
		self.__v4 = self.V(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt					
		)
		self.__a4 = self.A(
			t + self.dt, 
			self.r + self.__v3 * self.dt, 
			self.v + self.__a3 * self.dt			
		)
		#差分の計算
		self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
		self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 計算開始

#インスタンス
rk4 = RK4(dt, r0, v0)

#アニメーション作成用
ims=[]
def plot( t, r ):
	print( "{:.2f}".format(t),  np.sqrt(r[0]**2 + r[2]**2) - 1  )
	#各コマを描画
	img  = plt.plot([0], [0], colors[0], marker = 'o', markersize = 10, linestyle='solid', linewidth = 0 )
	img  += plt.plot([r[0]], [r[2]], colors[1], marker = 'o', markersize = 20, linestyle='solid', linewidth = 0 )
	img  += plt.plot([0, r[0]], [0, r[2]], colors[1], linestyle='solid', linewidth = 1 )
	time = plt.text(0.8, 1.25, "t = " +  str("{:.2f}".format(t)) + r"$[{\rm s}]$" , ha='left', va='center', fontsize=18)
	#テキストをグラフに追加
	img.append( time )
	ims.append( img )

#各時刻における運動方程式の計算
for tn in range(len(ts)):
	t = ts[tn]
	if( tn%skip == 0 ): 
		plot(t, rk4.r )
	#ルンゲ・クッタ法による時間発展
	rk4.timeEvolution( t )
	#位置と速度を更新
	rk4.r += rk4.dr
	rk4.v += rk4.dv


plt.title( u"単振り子運動", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(r"$x\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(r"$z\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)

#描画範囲を設定
plt.xlim([-1.2, 1.2])
plt.ylim([-1.2, 1.2])

#アニメーションの生成
#ani = animation.ArtistAnimation(fig, ims, interval=10)
#ani.save("output.html", writer=animation.HTMLWriter())
#ani.save("output.gif", writer="imagemagick")
#ani.save("output.mp4", writer="ffmpeg", dpi=300)

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

2021年8月17日2021年8月22日

【第１話】ルンゲクッタで行こう！①【Pythonコピペで古典力学完全攻略マニュアル】

############################################################################
#　単振動運動シミュレーション
############################################################################
import numpy as np
import matplotlib.pyplot as plt

#図全体
fig = plt.figure(figsize=(12, 8))
#全体設定
plt.rcParams['font.family'] = 'Times New Roman' #フォント
plt.rcParams['font.size'] = 24 #フォントサイズ
plt.rcParams["mathtext.fontset"] = 'cm' #数式用フォント
#カラーリストの取得
colors = plt.rcParams['axes.prop_cycle'].by_key()['color']

#################################################################
## 物理定数
#################################################################
#質量
m = 1.0
#ばね定数
k = 1.0

######################################
#　４次のルンゲ・クッタ　クラス
######################################
class RK4:
    #コンストラクタ
    def __init__(self, dt = 0.01, r0 = np.array([0.0, 0.0, 0.0]), v0 = np.array([0.0, 0.0, 0.0]) ):
        self.r = r0.copy()
        self.v = v0.copy()
        self.dr = np.array([0.0, 0.0, 0.0])
        self.dv = np.array([0.0, 0.0, 0.0])
        self.dt = dt
        self.__v1 = np.array([0.0, 0.0, 0.0])
        self.__v2 = np.array([0.0, 0.0, 0.0])
        self.__v3 = np.array([0.0, 0.0, 0.0])
        self.__v4 = np.array([0.0, 0.0, 0.0])
        self.__a1 = np.array([0.0, 0.0, 0.0])
        self.__a2 = np.array([0.0, 0.0, 0.0])
        self.__a3 = np.array([0.0, 0.0, 0.0])
        self.__a4 = np.array([0.0, 0.0, 0.0])

    #速度を与えるメソッド
    def V(self, t, r, v):
        return v.copy()
    ##########################################################
    #加速度を与えるメソッド
    def A(self, t, r, v):
        return -r.copy() * k / m
    ##########################################################

    #時間発展を計算するメソッド
    def timeEvolution(self, t):
        #１段目
        self.__v1 = self.V(
            t, 
            self.r, 
            self.v
        )
        self.__a1 = self.A(
            t, 
            self.r, 
            self.v
        );
        #２段目
        self.__v2 = self.V(
            t + self.dt / 2.0, 
            self.r + self.__v1 * self.dt / 2.0, 
            self.v + self.__a1 * self.dt / 2.0
        )
        self.__a2 = self.A(
            t + self.dt / 2.0, 
            self.r + self.__v1 * self.dt / 2.0, 
            self.v + self.__a1 * self.dt / 2.0            
        )
        #３段目
        self.__v3 = self.V(
            t + self.dt / 2.0, 
            self.r + self.__v2 * self.dt / 2.0, 
            self.v + self.__a2 * self.dt / 2.0                    
        )
        self.__a3 = self.A(
            t + self.dt / 2.0, 
            self.r + self.__v2 * self.dt / 2.0, 
            self.v + self.__a2 * self.dt / 2.0            
        )
        #４段目
        self.__v4 = self.V(
            t + self.dt, 
            self.r + self.__v3 * self.dt, 
            self.v + self.__a3 * self.dt                    
        );
        self.__a4 = self.A(
            t + self.dt, 
            self.r + self.__v3 * self.dt, 
            self.v + self.__a3 * self.dt            
        );        
        #差分の計算
        self.dr = self.dt * ( self.__v1 + 2.0 * self.__v2 + 2.0 * self.__v3 + self.__v4 ) / 6.0 
        self.dv = self.dt * ( self.__a1 + 2.0 * self.__a2 + 2.0 * self.__a3 + self.__a4 ) / 6.0 

#################################################################
## 物理系の設定
#################################################################

#初期位置と速度
r0 = np.array([0.0, 0.0, 10.0])
v0 = np.array([0.0, 0.0, 0.0])
#時間区間
t_min = 0
t_max = 20.0

#時間区間数
NT = 1000
skip = 10
#時間間隔
dt = (t_max - t_min) / NT
#座標点配列の生成
ts = np.linspace(t_min, t_max, NT + 1)
#インスタンス
rk4 = RK4(dt, r0, v0)
#グラフ描画用位置データ
tgs = []
rgs = []
ags = []

#解析解：初期位置のみ指定
def AnalyticalSolution( t ):
    omega = np.sqrt( k / m )
    return r0[2] *np.cos(omega * t)

#各時刻における運動方程式の計算
for tn in range(len(ts)):
    t = ts[tn]
    if( tn == 0 ):
        tgs.append( t ) 
        rgs.append( rk4.r[2] )
        ags.append( AnalyticalSolution( t ) )
        continue
    #ルンゲ・クッタ法による時間発展
    rk4.timeEvolution( t )
    #位置と速度を更新
    rk4.r += rk4.dr
    rk4.v += rk4.dv
    if( tn%skip == 0  ): 
        tgs.append( t ) 
        rgs.append( rk4.r[2] )
        ags.append( AnalyticalSolution( t ) )

plt.title( u"単振動運動（数値解と解析解）", fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.xlabel(u"時刻" + r"$\,[{\rm s}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
plt.ylabel(u"位置" + r"$\,[{\rm m}]$" , fontsize=20, fontname="Yu Gothic", fontweight=1000)
#罫線の描画
plt.grid(which = "major", axis = "x", alpha = 0.8, linestyle = "-", linewidth = 1)
plt.grid(which = "major", axis = "y", alpha = 0.8, linestyle = "-", linewidth = 1)
#スケールの指定
#ax = plt.gca()
#ax.set_yscale('log')
#ax.set_xscale('log')
#描画範囲を設定
plt.xlim([0, 20])
plt.ylim([-10.5, 10.5])
#プロット
plt.plot(tgs, ags, linestyle='solid', linewidth = 1 )
plt.plot(tgs, rgs, marker = 'o', markersize = 10, linestyle='solid', linewidth = 0 )
#グラフの表示
plt.show()