概要

同符号の電荷を持つ2物体がクーロン斥力を及ぼしあいながら運動する系を対象とし、1次元モデル[1]を2次元へ拡張したシミュレーションを構築した。運動方程式を数値的に解き、全運動量および力学的エネルギーの保存則が満たされることを確認するとともに、数値積分における時間刻み幅の選定が結果の精度に与える影響について検討した。

1. 序論

2物体が及ぼしあう斥力によって近接後に反発する現象は、完全弾性衝突の一種とみなすことができる。文献[1]では、同符号の電荷を持つ2物体を1次元上に配置し、クーロンの法則に基づく力を用いてこの過程を数値的に再現している。本稿では同モデルを2次元平面上の運動に拡張し、実装および数値計算上の妥当性を検証する。

2. モデル

質量 $m_A,\ m_B$ を持つ2物体 A, B が、2次元平面上の位置 $(x_A,y_A),\ (x_B,y_B)$ にあるとする。両者の距離を

r=\sqrt{(x_A-x_B)^2+(y_A-y_B)^2}

とし、電荷 $Q_A=Q_B=1$、クーロン定数 $k=1$ とすると、両物体間に働く力の大きさは

f=\frac{k Q_A Q_B}{r^2}=\frac{1}{r^2}

で与えられる。力は相手から遠ざかる向きに働くため、各物体の加速度は

\begin{aligned}
a_{x,A} &= \frac{f}{m_A}\cdot\frac{x_A-x_B}{r}, &
a_{y,A} &= \frac{f}{m_A}\cdot\frac{y_A-y_B}{r} \\[4pt]
a_{x,B} &= \frac{f}{m_B}\cdot\frac{x_B-x_A}{r}, &
a_{y,B} &= \frac{f}{m_B}\cdot\frac{y_B-y_A}{r}
\end{aligned}

と表される。両式はニュートンの第三法則 $\vec F_A=-\vec F_B$ を満たしており、外力が存在しない孤立系であることから、全運動量

\vec P = m_A\vec v_A + m_B\vec v_B

および力学的エネルギー(運動エネルギーと位置エネルギーの和)は時間に対して保存される。

3. 数値計算手法

運動方程式の時間積分には、速度を先に更新した後に位置を更新するセミインプリシット・オイラー法(シンプレクティック・オイラー法)を用いた。時刻 $t$ における各量から時刻 $t+dt$ における量を次式で計算する。

\begin{aligned}
v_{x,A}(t+dt) &= v_{x,A}(t) + a_{x,A}(t)\,dt, &
x_A(t+dt) &= x_A(t) + v_{x,A}(t+dt)\,dt
\end{aligned}

($y$ 成分および B についても同様)。

初期条件を表1に示す。

表1. 初期条件

変数	値
$m_A$	1
$m_B$	2
$(x_A,y_A)$	$(-1,\ 1)$
$(x_B,y_B)$	$(-1,\ -1)$
$(v_{x,A},v_{y,A})$	$(10,\ -10)$
$(v_{x,B},v_{y,B})$	$(10,\ 10)$
計算区間	$0 \le t < 0.4$

実装したPythonコードを以下に示す。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import matplotlib.animation as animation

fig = plt.figure()
ims = []

m_A = 1
m_B = 2

x_A, y_A = -1, 1
x_B, y_B = -1, -1

vx_A, vx_B = 10, 10
vy_A, vy_B = -10, 10

t = 0
dt = 1e-5
frame_interval = 100
step = 0

x_A_ary, x_B_ary = [], []
y_A_ary, y_B_ary = [], []
t_ary = []
Py_A_ary, Py_B_ary = [], []
K_A_ary, K_B_ary = [], []
f_ary = []

while t < 0.4:
    r = ((x_A - x_B) ** 2 + (y_A - y_B) ** 2) ** 0.5
    f = 1 / ((x_B - x_A) ** 2 + (y_B - y_A) ** 2)

    ax_A = f * (x_A - x_B) / (m_A * r)
    ay_A = f * (y_A - y_B) / (m_A * r)
    ax_B = f * (x_B - x_A) / (m_B * r)
    ay_B = f * (y_B - y_A) / (m_B * r)

    P_A = m_A * vy_A
    P_B = m_B * vy_B

    K_A = 0.5 * m_A * (vx_A ** 2 + vy_A ** 2)
    K_B = 0.5 * m_B * (vx_B ** 2 + vy_B ** 2)

    vx_A += ax_A * dt
    vy_A += ay_A * dt
    vx_B += ax_B * dt
    vy_B += ay_B * dt

    x_A += vx_A * dt
    x_B += vx_B * dt
    y_A += vy_A * dt
    y_B += vy_B * dt

    t += dt
    step += 1

    x_A_ary.append(x_A); x_B_ary.append(x_B)
    y_A_ary.append(y_A); y_B_ary.append(y_B)
    Py_A_ary.append(P_A); Py_B_ary.append(P_B)
    K_A_ary.append(K_A); K_B_ary.append(K_B)
    f_ary.append(f)
    t_ary.append(t)

    if step % frame_interval == 0:
        plt.xlim(-3, 3)
        plt.ylim(-10, 10)
        im1 = plt.plot(x_A, y_A, "o", color="red")
        im2 = plt.plot(x_B, y_B, "o", color="blue")
        ims.append(im1 + im2)

ani = animation.ArtistAnimation(fig, ims, interval=100)
ani.save("move_AB_2D.gif")
plt.show()

4. 時間刻み幅の検討

本モデルでは力が $f=1/r^2$ で距離に対して発散するため、最接近距離が小さいパラメータ条件では、最接近時刻付近において力が短時間で急激に変化する。この急激な変化を固定刻みのオイラー法で追跡する場合、時間刻み幅 $dt$ が粗いと数値誤差が顕著になり、保存されるべき全エネルギーが変動する。

表1の初期条件のもとで最接近距離は約 $r_\min = 7.5\times10^{-3}$ であり、$dt$ を変化させたときの全エネルギー $E=K_A+K_B+1/r$ の変動範囲を表2に示す。

表2. 時間刻み幅と全エネルギーの変動幅

$dt$	最接近距離 $r_\min$	全エネルギーの範囲
$1\times10^{-3}$	0.0050	268 〜 1157
$1\times10^{-4}$	0.0075	296 〜 306
$1\times10^{-5}$	0.0075	300.0 〜 301.0
$1\times10^{-6}$	0.0075	300.45 〜 300.55

(初期全エネルギーは約300.5である。)

この結果から、$dt=1\times10^{-3}$ では全エネルギーの変動が初期値の3倍以上に達し、保存則の検証には不適切であることが分かる。一方、$dt=1\times10^{-5}$ 以下では変動幅が初期値の1%未満に収まり、十分な精度が得られることを確認した。本稿では以降 $dt=1\times10^{-5}$ を採用する。

なお、$dt$ を縮小すると計算ステップ数が増加するため(本条件では40,000ステップ)、アニメーション用のフレーム取得を毎ステップ行うとフレーム数が過大となる。そこで物理量の時間積分は全ステップで行う一方、アニメーション用のフレームのみ frame_interval=100 ステップに1回の頻度で間引いて取得し、計算精度と描画コストを両立させた。

5. 結果

図1に2物体の軌跡を示す。物体Aは物体Bに向かって直進し、最接近付近で急激に進行方向を変えて離れていく、弾性衝突に特徴的な軌道が得られた。

図2に運動量の時間変化を示す。物体A・Bのy方向運動量は最接近時に急激に入れ替わるが、和は計算区間を通して一定に保たれている。

図3に運動エネルギーの時間変化を示す。最接近時に運動エネルギーの和が一時的に減少し(位置エネルギーへの変換に対応)、離脱後は初期値近傍に回復する。

図4に2物体間に働く力の時間変化を示す。最接近の瞬間にのみ鋭いピークを持ち、それ以外の区間ではほぼ0に近い値を取る。

図	ファイル
図1 軌跡	`AB軌跡.png`
図2 運動量	`運動量_2D.png`
図3 運動エネルギー	`運動エネルギー_2D.png`
図4 相互作用力	`物体に働く力_2D.png`
図5 アニメーション	`move_AB_2D.gif`

6. 結論

クーロン斥力を及ぼしあう2物体系を2次元に拡張してシミュレーションし、全運動量および力学的エネルギーが保存されることを数値的に確認した。また、力が距離の2乗に反比例して発散する系では、最接近距離に応じて時間刻み幅を十分に小さく取らなければ保存則が満たされないことを定量的に示した。

7. 今後の課題

距離に応じて時間刻み幅を動的に変化させる適応刻み法の導入
リープフロッグ法やRunge-Kutta-Nyström法など、より高精度な数値積分法との比較
3体以上の多体系への拡張

参考文献

[1] arairuca, "2つの正電荷が近づいて衝突するモデルをPythonで作ってみた", Qiita, https://qiita.com/arairuca/items/216bd1597297e08a45c5

付録: Pythonコード全文

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
import matplotlib.animation as animation

fig = plt.figure()

ims = []

m_A=1
m_B=2

x_A=-1
y_A=1

x_B=-1
y_B=-1

vx_A=10
vx_B=10

vy_A=-10
vy_B=10

t=0

dt=1e-5
frame_interval=100
step=0

x_A_ary=[]
x_B_ary=[]

y_A_ary=[]
y_B_ary=[]

vx_A_ary=[]
vx_B_ary=[]


vy_A_ary=[]
vy_B_ary=[]

ax_A_ary=[]
ax_B_ary=[]

ay_A_ary=[]
ay_B_ary=[]

t_ary=[]

Py_A_ary=[]
Py_B_ary=[]

K_A_ary=[]
K_B_ary=[]

f_ary=[]

while t<0.4:
    r=((x_A-x_B)**2+(y_A-y_B)**2)**0.5
    f=1/((x_B-x_A)**2+(y_B-y_A)**2)
    ax_A=f*(x_A-x_B)/(m_A*r)
    ay_A=f*(y_A-y_B)/(m_A*r)
    ax_B=f*(x_B-x_A)/(m_B*r)
    ay_B=f*(y_B-y_A)/(m_B*r)

    P_A=m_A*vy_A
    P_B=m_B*vy_B
    
    K_A=0.5*m_A*(vx_A**2+vy_A**2)
    K_B=0.5*m_B*(vx_B**2+vy_B**2)
    
    vx_A=vx_A+ax_A*dt
    vy_A=vy_A+ay_A*dt
    
    vx_B=vx_B+ax_B*dt
    vy_B=vy_B+ay_B*dt
    
    x_A=x_A+vx_A*dt
    x_B=x_B+vx_B*dt
    
    y_A=y_A+vy_A*dt
    y_B=y_B+vy_B*dt
    
    t=t+dt
    step=step+1
    ax_A_ary.append(ax_A)
    ax_B_ary.append(ax_B)
    
    vx_A_ary.append(vx_A)
    vx_B_ary.append(vx_B)

    x_A_ary.append(x_A)
    x_B_ary.append(x_B)
    
    ay_A_ary.append(ay_A)
    ay_B_ary.append(ay_B)
    
    vy_A_ary.append(vy_A)
    vy_B_ary.append(vy_B)

    y_A_ary.append(y_A)
    y_B_ary.append(y_B)
    
    
    Py_A_ary.append(P_A)
    Py_B_ary.append(P_B)
    
    K_A_ary.append(K_A)
    K_B_ary.append(K_B)
    
    f_ary.append(f)
    
    t_ary.append(t)
    
    
    if step%frame_interval==0:
        plt.xlim(-3,3)
        plt.ylim(-10,10)
        im1=plt.plot(x_A,y_A,"o",color="red")
        im2=plt.plot(x_B,y_B,"o",color="blue")
        ims.append(im1+im2)  # グラフを配列 ims に追加


# 10枚のプロットを 100ms ごとに表示
ani = animation.ArtistAnimation(fig, ims, interval=100)
 # 保存
ani.save('move_AB_2D.gif')
plt.show()

plt.plot(t_ary,Py_A_ary,color="red",label="A")
plt.plot(t_ary,Py_B_ary,color="blue",label="B")
plt.xlabel("時刻")
plt.ylabel("運動量")
plt.legend()
plt.savefig("運動量_2D.png")
plt.show()

plt.plot(t_ary,K_A_ary,color="red",label="A")
plt.plot(t_ary,K_B_ary,color="blue",label="B")
plt.xlabel("時刻")
plt.ylabel("運動エネルギー")
plt.legend()
plt.savefig("運動エネルギー_2D.png")
plt.show()


plt.plot(t_ary,f_ary,color="blue")

plt.xlabel("時刻")
plt.ylabel("2物体に働く力")
plt.legend()
plt.savefig("物体に働く力_2D.png")
plt.show()


plt.plot(x_A_ary,y_A_ary,label="A")
plt.plot(x_B_ary,y_B_ary,label="B")
plt.savefig("AB軌跡.png")
plt.show()