【Python】保存量補正力の導入

この前に紹介した伸び縮みしないひもによる張力による単振り子のように、ひもの長さのように本来絶対に保存してほしい量がある数値計算の場合、計算誤差の蓄積で破綻してしまうことがよくあるね。陽的解法である４次のルンゲクッタなどのアルゴリズムの場合、誤差が単調増加することが知られているね。今回利用しているPythonのscipy.integrate.solve_ivpは１ステップあたりの相対誤差と絶対誤差を設定できるので問題ないように思えるけれども、単調増加することは避けられないので、抜本的な対応を必要だね。そこで、ひもの長さを補正するための力（ひもの長さ補正力）を次のとおり定義して実装してみるね。

$$
\boldsymbol{f}_{\rm c} = -k_{\rm c}\, (|\boldsymbol{r}|-r_0) \hat{\boldsymbol{r}} - \beta_{\rm c}(\bar{\boldsymbol{v}}\cdot \hat{\boldsymbol{r}})\hat{\boldsymbol{r}}
$$

第１項目はもとの長さに戻すための補正ばね弾性力、第２項目はひもの長さ方向の速度変化を抑えるための補正粘性抵抗力だよ。$${k_c}$$と$${\beta_c}$$はそれぞれ補正ばね弾性係数と補正粘性抵抗係数だけれども、これは単なる定数ではなくて、ひもに加わる張力の大きさ $${S}$$ に比例した値とすべきだね（そうじゃないと張力が大きくなるほど補正力が相対的に小さくなってしまうからね）。

次のグラフは計算時間に対するひもの長さの誤差$$|\Delta L|$$だよ（相対誤差：$${10^{-5}}$$、絶対誤差：$${10^{-5}}$$）。補正力ゼロ（$${k_c = 0,\ \beta_c=0}$$）、補正粘性力のみ（$${k_c = 0,\ \beta_c=0.5}$$）、補正弾性力のみ（$${k_c = 0.3,\ \beta_c=0}$$）、補正弾性力＋補正粘性力（$${k_c = 0.3,\ \beta_c=0.5}$$）の４パターンを示しているよ。補正力なしの場合、100秒後のひもの長さの誤差がすでに２０％以上生じてしまっている状態で、補正粘性力を加えると大分抑え込めてるいるけれども、単調増加は抑えられないね。補正ばね弾性力を加えると誤差が一定値に抑えることができていることがわかるね。

ちなみにscipy.integrate.solve_ivpの相対誤差と絶対誤差を$${10^{-10}}$$と設定したときの計算結果は次のとおりだよ。計算時間は大幅に増えるけれども誤差の大きさは抑えられてるね。だけれども補正力がない場合の単調増加となる傾向は同じだね。やはり補正力は大切だね。

プログラムソース（Python）：単振り子

#################################################################
TITLE = "張力＋重力による単振子運動"
#################################################################
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   = 100.0       #計算終了時刻
DT = 0.01               #時間間隔
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])
#質量
m = 1.0
#粘性係数
gamma = 0.0
#ひもの長さ
L0 = 1
#初期速度
vec_v0 = np.array([0, 0, 0])
#初期角度
theta = np.pi / 2.0
###############################
# 関数定義
###############################
#張力＋重力の微分方程式
def equations1(t, X , vec_g, m, compensationK, compensationGamma):
	vec_r =  np.array([ X[0], X[1], X[2] ])
	vec_v =  np.array([ X[3], X[4], X[5] ])
	#紐の長さ
	L = np.sqrt( np.dot(vec_r, vec_r) )
	#速度の大きさ
	v = np.sqrt( np.dot(vec_v, vec_v) )
	#張力
	vec_S = - m / L**2 * ( v**2 + np.dot(vec_g, vec_r) ) * vec_r
	#位置ベクトルと速度ベクトルの単位ベクトル
	hat_r = vec_r / L
	if( v > 0 ): hat_v = vec_v / v
	else: hat_v = vec_v * 0
	#補正力
	f_c = - compensationK * ( L - L0 ) * hat_r - compensationGamma * np.dot(hat_v, hat_r) * hat_r
	#質点に加わる力
	F = vec_S + m * vec_g + f_c * np.sqrt( np.dot(vec_S, vec_S) )
	#ニュートンの運動方程式
	vec_a = F / m
	return [vec_v[0], vec_v[1], vec_v[2], vec_a[0], vec_a[1], vec_a[2]]

# theta 依存性を可視化
def calculation( L0, compensationK, compensationGamma ):
	#初期位置
	vec_r0 = np.array([ L0 * np.sin( theta ), 0, L0 * np.cos( theta )])
	solution = integrate.solve_ivp(
		equations1,                     #常微分方程式
		method = 'RK45',                #計算方法
		args = (vec_g, m, compensationK, compensationGamma),              #引数
		t_span=( times[0], times[-1]),  #積分区間
		rtol = 1.e-5,                   #相対誤差
		atol = 1.e-5,                   #絶対誤差
		y0 = np.concatenate([vec_r0, vec_v0]), #初期値
		t_eval = times                  #結果を取得する時刻
	)
	#計算結果
	return solution.y[ 0 ], solution.y[ 1 ], solution.y[ 2 ]

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

	x1, y1, z1 = calculation( L0, 0.0, 0.0 )
	x2, y2, z2 = calculation( L0, 0.0, 0.5 )
	x3, y3, z3 = calculation( L0, 0.3, 0.0 )
	x4, y4, z4 = calculation( L0, 0.3, 0.5 )

	#描画用リスト
	y1l = []
	y2l = []
	y3l = []
	y4l = []
	for i in range(len(times)):
		L1 = np.sqrt( x1[i]**2 + y1[i]**2 + z1[i]**2 )
		L2 = np.sqrt( x2[i]**2 + y2[i]**2 + z2[i]**2 )
		L3 = np.sqrt( x3[i]**2 + y3[i]**2 + z3[i]**2 )
		L4 = np.sqrt( x4[i]**2 + y4[i]**2 + z4[i]**2 )
		y1l.append( abs(L1 - L0) )
		y2l.append( abs(L2 - L0) )
		y3l.append( abs(L3 - L0) )
		y4l.append( abs(L4 - L0) )

	#################################################################
	# グラフの設定
	#################################################################
	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, 0])
	plt.yscale("log")  

	plt.plot(times, y1l, linestyle='solid', linewidth = 3)
	plt.plot(times, y2l, linestyle='solid', linewidth = 3)
	plt.plot(times, y3l, linestyle='solid', linewidth = 3)
	plt.plot(times, y4l, linestyle='solid', linewidth = 3)
	#グラフの表示
	plt.show()