import numpy as np
import matplotlib.pyplot as plt
import math
h = 0.01 # 時間刻み幅
xe = 0.3 # バネの自然長[m]
k = 50 # バネ定数[N/m]
m = 0.5 # おもりの質量[kg]
v0 = 0 # おもりの初速度[m/s](静かに手を離すので必ず0)
def f_v(t, x, v):
""" dv/dt """
return - k / m * (x - xe)
def f_x(t, x, v):
""" dx/dt """
return v
def calc_numerical_solution(t, x, v):
""" 4次のルンゲ・クッタ法により次ステップの値を計算する """
k1 = f_v(t, x, v) * h
l1 = f_x(t, x, v) * h
k2 = f_v(t + h / 2, x + l1 / 2, v + k1 / 2) * h
l2 = f_x(t + h / 2, x + l1 / 2, v + k1 / 2) * h
k3 = f_v(t + h / 2, x + l2 / 2, v + k2 / 2) * h
l3 = f_x(t + h / 2, x + l2 / 2, v + k2 / 2) * h
k4 = f_v(t + h, x + l3, v + k3) * h
l4 = f_x(t + h, x + l3, v + k3) * h
v_next = v + (k1 + 2 * k2 + 2 * k3 + k4) / 6
x_next = x + (l1 + 2 * l2 + 2 * l3 + l4) / 6
return x_next, v_next
def calc_theoretical_solution(t):
""" 理論解を計算する """
omega = math.sqrt(k / m)
x = (x0 - xe) * math.cos(omega * t) + xe
v = -(x0 - xe) * omega * math.sin(omega * t)
return x, v
def simple_harmonic_motion(x0, t_list):
""" 単振動の計算 """
xn = x0
vn = v0
# t=0の理論解を計算する
xt, vt = calc_theoretical_solution(0)
# 初期条件を解のリストに追加
xn_list = [xn] # xの数値解リスト
vn_list = [vn] # vの数値解リスト
xt_list = [xt] # xの理論解リスト
vt_list = [vt] # vの理論解リスト
# t>0について解の計算を行う
for t in t_list[1:]:
# 数値解を計算
xn, vn = calc_numerical_solution(t, xn, vn)
xn_list.append(xn)
vn_list.append(vn)
# 理論解を計算
xt, vt = calc_theoretical_solution(t)
xt_list.append(xt)
vt_list.append(vt)
return xn_list, vn_list, xt_list, vt_list
def plot_graph(xn_list, vn_list, xt_list, vt_list):
fig = plt.figure(figsize=(8, 4))
fig.suptitle('Horizontal simple harmonic motion')
# 位置のプロット
plt.subplot(2, 1, 1)
plt.plot(t_list, xn_list, label='numerical',
color='r', linewidth=1) # 数値解
plt.plot(t_list, xt_list, label='theoretical', color='g', linestyle='None',
marker='o', markersize=1) # 理論解
plt.hlines(xe, t_list[0], t_list[-1], linestyles='dotted',
colors='k', linewidth=1) # 自然長の線
plt.ylabel("x [m]")
# 凡例
plt.legend(['numerical', 'theoretical'],
loc='upper left', bbox_to_anchor=(1, 1))
# 速度のプロット
plt.subplot(2, 1, 2)
plt.plot(t_list, vn_list, color='r', linewidth=1) # 数値解
plt.plot(t_list, vt_list, color='g', linestyle='None',
marker='o', markersize=1) # 理論解
plt.hlines(0, t_list[0], t_list[-1], linestyles='dotted',
colors='k', linewidth=1) # 速度0の線
plt.xlabel("time [s]")
plt.ylabel("velocity [m/s]")
fig.tight_layout()
plt.show()
if __name__ == '__main__':
x0 = float(input('初期位置[x]:'))
t_end = float(input('計算終了時刻[s]:'))
# 計算時間ステップのリストを生成する(最後の時刻t_endも含めたいので+hしている)
t_list = np.arange(0, t_end + h, h)
# 単振動の計算
xn_list, vn_list, xt_list, vt_list = simple_harmonic_motion(x0, t_list)
plot_graph(xn_list, vn_list, xt_list, vt_list)