【Python】任意多重結合における長さ補正力の導入
前にも述べたけれども、ルンゲ・クッタ法では誤差が単調増加してしまうので、どんなに誤差設定（相対誤差・絶対誤差）を小さくしても長い経過時間後には計算が破綻してしまうね。そこで、以前にも導入したひもの長さを初期値に戻す補正力を導入するよ。下図は$${i}$$番目と$${j}$$番目のおもりをつなぐひもの補正力を表す模式図だよ。
ひもの長さの初期値を$${ L^{(0)}_{ij}}$$、相対速度ベクトルを$${\bar{\boldsymbol{v}}_{ij}}$$、相対位置ベクトルの単位ベクトルを$${\hat{\boldsymbol{n}}_{ij}}$$としたときの、補正力は次のように与えられるね。
$$
\boldsymbol{f}^{( c)}_{ij} = \left[ k_{\rm c}\, (L_{ij}-L^{(0)}_{ij}) + \beta_{\rm c}(\bar{\boldsymbol{v}}_{ij}\cdot \hat{\boldsymbol{n}}_{ij} ) \right]\hat{\boldsymbol{n}}_{ij}
$$
$${ k_c, \beta_c}$$は、それぞれ補正ばね弾性係数、補正粘性係数だね。このひもの長さ補正力を用いて、２重振子の各おもりに加わる力は次のようになるね。
$$
\boldsymbol{F}_1 = \boldsymbol{f}_{01} + \boldsymbol{f}_{21}+ \boldsymbol{F}^{\,({\rm ex})}_1 + \alpha_1\left( \boldsymbol{f}^{(c)}_{01}+ \boldsymbol{f}^{(c)}_{21} \right)
$$
$$
\boldsymbol{F}_2 = \boldsymbol{f}_{12}+ \boldsymbol{F}^{\,({\rm ex})}_2 + \alpha_2 \boldsymbol{f}^{(c)}_{12}
$$
ただし、おもりに加わっている力が大きいほど補正力が大きくなってほしいので、補正力の係数を次のとおりに導入するね。
$$
\alpha_1=| \boldsymbol{f}_{01}| +| \boldsymbol{f}_{12}|
$$
$$
\alpha_2=| \boldsymbol{f}_{12}|
$$
計算結果
先の２重振子運動をルンゲクッタ法（相対誤差・絶対誤差：$${10^{-4}}$$）と設定して、補正ばね弾性係数（$${k_c}$$）、補正粘性係数（$${\beta_c}$$）を変えてひもの長さの誤差（$${|\Delta L|}$$）の時間変化（100秒間）を計算するよ。あえて計算誤差が大きくなるように２つのおもりの質量を$${m_1=10, m_2=100}$$としているよ。
上図は、$${ k_c =0}$$に固定して、$${\beta_c = 0\sim0.4}$$まで変化させたときの結果だよ。予想どおり、$${\beta_c }$$を大きくするほど誤差が小さくなっているね。$${\beta_c =0}$$のときは誤差が２程度と、既に計算が破綻している状況だけれども、$${\beta_c =0.1}$$で誤差は1/10を達成しているね。$${\beta_c =0.4}$$で誤差は1/100近くまで減らせているね。ひもの長さ方向の速度変化を減衰させるだけですごい効果だね。
補正力が大きくなればなるほど、ルンゲクッタ法の計算時間が増大するので、$${\beta_c }$$はできるだけ小さな値のほうが良いので、とりあえず$${\beta_c =0.2}$$で固定しておいて、今度は$${k_c =0\sim0.4}$$まで変化させて計算した結果が上図だよ。$${k_c }$$が大きくなるほど、ひもの長さの誤差が減っているね。$${k_c =0.03}$$で誤差1/100を十分に達成しているね。これでルンゲクッタ法を用いても長時間計算が破綻せずに済むね。
【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 = 100.0 #計算終了時刻
DT = 0.01 #時間間隔
###############################
# 計算グローバル変数
###############################
#時刻
times = np.arange(TIME_START, TIME_END, DT)
#重力加速度
vec_g = np.array([0, 0, -10])
#質量
m1 = 1.0
m2 = 1.0
#粘性係数
gamma = 0.0
#ひもの長さ
L01_0 = 1.0
L12_0 = 1.0
#初速度
v0 = 0.01
#支点の振幅と角振動数
l = 0.0
omega_box = 0 * np.pi
###############################
# 関数定義
###############################
#支点の位置と速度
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 equations(t, X, compensationK, compensationGamma):
vec_r1 = np.array([ X[0], X[1], X[2] ])
vec_v1 = np.array([ X[3], X[4], X[5] ])
vec_r2 = np.array([ X[6], X[7], X[8] ])
vec_v2 = np.array([ X[9], X[10], X[11] ])
#支点の位置ベクトル・速度ベクトル・加速度ベクトル
vec_r0, vec_v0, vec_a0 = box( t )
#相対ベクトル
vec_r01 = vec_r0 - vec_r1
vec_r12 = vec_r1 - vec_r2
#ひもの長さ
L01 = np.sqrt( np.dot(vec_r01, vec_r01) )
L12 = np.sqrt( np.dot(vec_r12, vec_r12) )
#相対単位ベクトル
vec_n01 = vec_r01 / L01
vec_n12 = vec_r12 / L12
#相対速度ベクトル
vec_v01 = vec_v0 - vec_v1
vec_v12 = vec_v1 - vec_v2
#相対速度の大きさ
v01 = np.sqrt( np.dot(vec_v01, vec_v01) )
v12 = np.sqrt( np.dot(vec_v12, vec_v12) )
#外力による加速度ベクトル
vec_a1 = vec_g
vec_a2 = vec_g
#相対加速度ベクトル
vec_a01 = vec_a0 - vec_a1
vec_a12 = vec_a1 - vec_a2
#A係数
A01 = -1.0 / L01 * ( v01**2 + np.dot(vec_a01, vec_r01) )
A12 = -1.0 / L12 * ( v12**2 + np.dot(vec_a12, vec_r12) )
#張力係数
f01 = - ( m1 * (m1 + m2) / ( m1 + m2 * ( 1.0 - np.dot( vec_n01, vec_n12)**2 ) ) ) * A01 - ( m1 * m2 * np.dot( vec_n01, vec_n12)/ ( m1 + m2 * ( 1.0 - np.dot( vec_n01, vec_n12)**2 ) ) ) * A12
f12 = - ( m1 * m2 * np.dot( vec_n01, vec_n12) / ( m1 + m2 * ( 1.0 - np.dot( vec_n01, vec_n12)**2 ) ) ) * A01 - ( m1 * m2 / ( m1 + m2 * ( 1.0 - np.dot( vec_n01, vec_n12)**2 ) ) ) * A12
#位置ベクトルと速度ベクトルの単位ベクトル
hat_n01 = vec_r01 / L01
hat_n12 = vec_r12 / L12
#補正力
vec_f01_c = (compensationK * ( L01 - L01_0 ) + compensationGamma * np.dot(vec_v01, hat_n01) ) * hat_n01
vec_f12_c = (compensationK * ( L12 - L12_0 ) + compensationGamma * np.dot(vec_v12, hat_n12) ) * hat_n12
#質点に加わる力
vec_F1 = f01 * vec_n01 - f12 * vec_n12 + m1 * vec_g + (vec_f01_c - vec_f12_c) * ( np.abs(f01) + np.abs(f12) )
vec_F2 = f12 * vec_n12 + m2 * vec_g + vec_f12_c * ( np.abs(f12) )
#ニュートンの運動方程式
vec_a1 = vec_F1 / m1
vec_a2 = vec_F2 / m2
return [
vec_v1[0], vec_v1[1], vec_v1[2], vec_a1[0], vec_a1[1], vec_a1[2],
vec_v2[0], vec_v2[1], vec_v2[2], vec_a2[0], vec_a2[1], vec_a2[2]
]
# theta 依存性を可視化
def calculation(vec_r1, vec_v1, vec_r2, vec_v2, compensationK, compensationGamma):
solution = integrate.solve_ivp(
equations, #常微分方程式
method = 'RK45', #計算方法
args = (compensationK, compensationGamma), #引数
t_span=( times[0], times[-1]), #積分区間
rtol = 1.e-4, #相対誤差
atol = 1.e-4, #絶対誤差
y0 = np.concatenate([vec_r1, vec_v1, vec_r2, vec_v2]), #初期値
t_eval = times #結果を取得する時刻
)
#計算結果
vec_r1_t = np.array( [solution.y[0], solution.y[1], solution.y[2]] )
vec_v1_t = np.array( [solution.y[3], solution.y[4], solution.y[5]] )
vec_r2_t = np.array( [solution.y[6], solution.y[7], solution.y[8]] )
vec_v2_t = np.array( [solution.y[9], solution.y[10], solution.y[11]] )
return vec_r1_t, vec_v1_t, vec_r2_t, vec_v2_t
##########################################################
if __name__ == "__main__":
def calError(compensationK, compensationGamma ):
#初期位置と速度
vec_r1 = np.array([ 0, 0, L01_0])
vec_v1 = np.array([ 0, 0, 0])
vec_r2 = np.array([ 0, 0, L01_0 + L12_0])
vec_v2 = np.array([ v0, 0, 0])
vec_r1_t, vec_v1_t, vec_r2_t, vec_v2_t = calculation( vec_r1, vec_v1, vec_r2, vec_v2, compensationK, compensationGamma)
result = []
ts = []
for i in range(len(times)):
t = i * DT
vec_r_box, vec_v_box, vec_a_box = box( t )
L01 = np.sqrt( (vec_r1_t[ 0 ][ i ] - vec_r_box[ 0 ])**2 + (vec_r1_t[ 1 ][ i ] - vec_r_box[ 1 ])**2 + (vec_r1_t[ 2 ][ i ] - vec_r_box[ 2 ])**2 )
L12 = np.sqrt( (vec_r2_t[ 0 ][ i ] - vec_r1_t[ 0 ][ i ])**2 + (vec_r2_t[ 1 ][ i ] - vec_r1_t[ 1 ][ i ])**2 + (vec_r2_t[ 2 ][ i ] - vec_r1_t[ 2 ][ i ])**2 )
DeltaL = np.abs( L01 - L01_0 ) + np.abs( L12 - L12_0 )
if(i%50==49):
print( L01, L12 )
ts.append( t )
result.append( DeltaL )
return result, ts
# d0, ts = calError(0.0, 0.0)
# d1, ts = calError(0.0, 0.1)
# d2, ts = calError(0.0, 0.2)
# d3, ts = calError(0.0, 0.3)
# d4, ts = calError(0.0, 0.4)
d0, ts = calError(0.00, 0.2)
d1, ts = calError(0.01, 0.2)
d2, ts = calError(0.02, 0.2)
d3, ts = calError(0.03, 0.2)
d4, ts = calError(0.04, 0.2)
#################################################################
# グラフの設定
#################################################################
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, d0, linestyle='solid', linewidth = 4)
plt.plot(ts, d1, linestyle='solid', linewidth = 4)
plt.plot(ts, d2, linestyle='solid', linewidth = 4)
plt.plot(ts, d3, linestyle='solid', linewidth = 4)
plt.plot(ts, d4, linestyle='solid', linewidth = 4)
#グラフの表示
plt.show()