バネで固定された粒子のシミュレーション

確率過程のシミュレーションを考える前に,今回はひとまず通常の確率的ではない(=決定論的)運動のシミュレーションを扱う.例題として原点にバネで固定された粒子のシミュレーションをオイラー法を用いて行う.粒子はバネから受ける力の他に,速度に比例した抵抗力を受けるものとする.

運動方程式

(1)   dR(t)/dt = V(t)
(2)   m(dV(t)/dt) = - zV(t) - kR(t)
  • R(t): 粒子の位置座標
  • V(t): 粒子の速度
  • t: 時間
  • m: 粒子の質量
  • k: バネ定数
  • z: 粒子の摩擦係数

オイラー法

時刻tを小さな間隔Dtで区切り,t=t_i=Dt*i (i=0,1,2,...)と離散化して考え,R(t_i)をR_iと表示する.

オイラー法では,以下の様に微分を前進差分で置き換える.

dR(t)/dt = (R_i+1 - R_i)/Dt

dV(t)/dt = (V_i+1 - V_i)/Dt

よって,運動方程式は以下の様に差分形式に変換できる(差分形式はコンピュータが得意な形)

R_i+1 = R_i + V_i * Dt

V_i+1 = V_i * (1 - z/m * Dt) - k/m * Dt * R_i

使用するライブラリのインポート

% matplotlib nbagg
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

変数の設定

dim  = 2                            # 系の次元
nums = 1000                         # シミュレーションのステップ数
R    = np.zeros(dim)                # 粒子の位置座標(瞬間値)
V    = np.zeros(dim)                # 粒子の速度(瞬間値)
Rs   = np.zeros([dim,nums])         # 粒子の位置座標(各ステップの値)
Vs   = np.zeros([dim,nums])         # 粒子の速度(各ステップの値)
Et   = np.zeros(nums)               # 系の全エネルギー(各ステップの値)
time = np.zeros(nums)               # 時刻(各ステップの値)

オイラー法でシミュレーション&動画表示

# System parameters
m    = 1.0                          # 粒子の質量
k    = 1.0                          # バネ定数
zeta = 0.0                          # 粒子の摩擦係数
# Initial condition
R[0] = 1.                           # 粒子位置の初期値(x座標)
R[1] = 0.                           # 粒子位置の初期値(y座標)
V[0] = 0.                           # 粒子速度の初期値(x座標)
V[1] = 0.                           # 粒子速度の初期値(y座標)
# set dt
dt   = 0.05*np.sqrt(k/m)             # dtの設定(バネの振動数の定数倍とする)
# set box size
box  = 5                            # 動画作成領域のサイズ
# First set up the figure, the axis, and the plot element we want to animate
fig, ax = plt.subplots(figsize=(7.5,7.5))
ax = plt.axes(xlim=(-box/2,box/2),ylim=(-box/2,box/2))
particles, = ax.plot([], [], 'ro', ms=10)
line, = ax.plot([], [], lw=1)
title = ax.text(0.5, 1.05, '', transform = ax.transAxes, va='center')
xdata, ydata = [], []

def init():
    particles.set_data([], [])
    line.set_data([], [])
    title.set_text('')
    return particles,line,title 

def animate(i):
    global R,V,F,Rs,Vs,time,Et
    V = V*(1-zeta/m*dt)-k/m*dt*R   # オイラー法で粒子速度の積分
#    V = V*(1-zeta/m*dt)-k/m*dt*np.linalg.norm(R)**2*R   # オイラー法で粒子速度の積分
    R = R + V*dt                    # オイラー法で粒子位置の積分
    xdata.append(R[0])
    ydata.append(R[1])
    particles.set_data(R[0], R[1])
    line.set_data(xdata, ydata)
    title.set_text("t = "+str(i*dt)+"  Et = "+str(Et[i]))
    Rs[0:dim,i]=R                   # 粒子位置の瞬間値を保存
    Vs[0:dim,i]=V                   # 粒子速度の瞬間値を保存
    time[i]=i*dt                    # 時刻の瞬間値を保存
    Et[i]=0.5*m*np.linalg.norm(V)**2+0.5*k*np.linalg.norm(R)**2   # 系の全エネルギーの瞬間値を保存
#    Et[i]=0.5*m*np.linalg.norm(V)**2+0.25*k*np.linalg.norm(R)**4   # 系の全エネルギーの瞬間値を保存
    return particles,line,title

anim = animation.FuncAnimation(fig, animate, init_func=init,
                               frames=nums, interval=5, blit=True, repeat=False)

保存された種々のデータ(A(t), B(t), ...)を時刻(t)の関数としてプロット

fig, ax = plt.subplots(figsize=(7.5,7.5))
ax.set_xlabel("$t$", fontsize=20)
ax.plot(time,Et)                # 各時刻の系の全エネルギーをプロット(摩擦がなければ一定,あれば減少するはず)
ax.plot(time,Rs[0])             # 各時刻の粒子位置(x座標)をプロット(動画で見たとおりに振動するはず)
ax.plot(time,Rs[1])             # 各時刻の粒子位置(y座標)をプロット(動画で見たとおりに振動するはず)
ax.legend(['$E_t(t)$','$R_x(t)$','$R_y(t)$'], fontsize=14)
plt.show()

A(t) vs B(t) のようなパラメトリックプロットも可能

fig, ax = plt.subplots(figsize=(7.5,7.5))
ax.plot(Rs[0,0:nums],Rs[1,0:nums])  # これで粒子の軌跡が平面上にパラメトリックプロットされるはず
plt.show()

宿題

  • シミュレーションの設定(粒子の初期条件,バネ定数,粒子質量,摩擦係数,バネの線形/非線形などなど)を変更して何度か実行し,目的のシミュレーションが実行できていることを確認する.その1つについて系の全エネルギーを時刻の関数としてプロットしたグラフをメールで提出すること.シミュレーションの設定条件もメールに明記すること.