２つの物体の運動（２体問題）のシミュレーション
Feynmanら（1963、p. 9-6）の方法による
質量
と
の２つの物体が引力
で引き合っているときの運動のシミュレーションプログラムを、岡本（1997）のpp. 139-142の解説に従って作成してみた。
次式
で表される運動方程式を、Feynmanら（1963）のp. 9-6（第９章の6ページ目）に説明されている方法に基づいて差分で表し、リスト１に示すプログラムを作成した。リスト１のプログラムをNumpyのndarray型を用いたものに変更したものは、リスト２に示した。
プログラムファイルは圧縮ファイルtwomovingobjsfile.zipとしたので、ダウンロード展開すれば自由に利用できる。
リスト１のプログラムでは、差分による更新は以下のようにコーディングされている。ただし、簡単のために
としている。
x0_t = x0 + d_t * vx0h
y0_t = y0 + d_t * vy0h
x1_t = x1 + d_t * vx1h
y1_t = y1 + d_t * vy1h
r3 = ((x0_t - x1_t)**2 +
(y0_t - y1_t)**2) ** (3/2)
vx0h_t = vx0h + d_t * (-8) *
(x0_t - x1_t) / r3
vy0h_t = vy0h + d_t * (-8) *
(y0_t - y1_t) / r3
vx1h_t = vx1h + d_t * (x0_t
- x1_t) / r3
vy1h_t = vy1h + d_t * (y0_t
- y1_t) / r3
更新後の２物体の位置は、次のコードにより画面に表示される。
screen.fill((0, 0, 0))
pg.draw.circle(screen, (255,
200,255), (x_scale(x0), y_scale(y0)), 5)
pg.draw.circle(screen, (0, 255,
255), (x_scale(x1), y_scale(y1)), 20)
pg.display.flip()
リスト１のプログラムの実行例を、図１に示す。
図１ CentOS Stream 9上で実行
実行環境は、Python3.12/Anaconda/CentOS Stream 9である。プログラムを実行すると、画面内で２つの物体（青と黄）の運動が見られるが、表示画面の右上角のX印をクリックするとプログラムの実行終了となる。
Feynman, R. P., Leighton, R. B., & Sands, M. (1963). The Feynman Lectures on Physics, Vol. I. Addison-Wesley Publishing Company.
岡本安晴（1997）Delphiプログラミング入門――Object Pascalの基礎から始める――．CQ出版社
#
#
Simulation of
# Two Moving
Bodies
#
import numpy as np
import pygame as pg
import sys
pg.init()
screen = pg.display.set_mode(size=(800,700))
screen.fill((0, 0, 0))
pg.display.flip()
def x_scale(x):
return 400 + x * 100
def y_scale(y):
return 350 + y * 100
#
# The
initial values
#
x0 = 3.0
y0 = 0.0
x1 = 0.0
y1 = 0.0
vx0h = 0.0
vy0h = -0.8
vx1h = 0.0
vy1h = 0.1
d_t = 0.001
while True:
#
# Calculation of values
for update
#
x0_t = x0 + d_t * vx0h
y0_t = y0 + d_t * vy0h
x1_t = x1 + d_t * vx1h
y1_t = y1 + d_t * vy1h
r3 = ((x0_t - x1_t)**2 + (y0_t
- y1_t)**2) ** (3/2)
vx0h_t = vx0h + d_t * (-8) *
(x0_t - x1_t) / r3
vy0h_t = vy0h + d_t * (-8) *
(y0_t - y1_t) / r3
vx1h_t = vx1h + d_t * (x0_t
- x1_t) / r3
vy1h_t = vy1h + d_t * (y0_t
- y1_t) / r3
#
#
Update the values
# cf. Feynman et al. (1963),
p. 9-6.
#
x0 = x0_t
y0 = y0_t
x1 = x1_t
y1 = y1_t
vx0h = vx0h_t
vy0h = vy0h_t
vx1h = vx1h_t
vy1h = vy1h_t
#
#
Display the updated values
#
screen.fill((0, 0, 0))
pg.draw.circle(screen, (255,
200,255), (x_scale(x0), y_scale(y0)), 5)
pg.draw.circle(screen, (0,
255, 255), (x_scale(x1), y_scale(y1)), 20)
pg.display.flip()
for e in pg.event.get():
if
e.type == pg.QUIT:
# Is the windows's
close button clicked?
pg.quit()
sys.exit()
pg.time.wait(1)
#
#
Simulation of
# Two Moving
Bodies
#
import numpy as np
import pygame as pg
import sys
pg.init()
screen =
pg.display.set_mode(size=(800,700))
screen.fill((0, 0, 0))
pg.display.flip()
def x_scale(x):
return 400 + x * 100
def y_scale(y):
return 350 + y * 100
#
# The
initial values
#
x0 = np.array([3.0, 0.0])
x1 = np.array([0.0,0.0])
v0h = np.array([0.0, -0.8])
v1h = np.array([0.0, 0.1])
d_t = 0.001
while True:
#
# Calculation of values
for update
#
x0_t = x0 + d_t * v0h
x1_t = x1 + d_t * v1h
r2 = np.sum((x0_t -
x1_t)**2)
v0h_t = v0h + d_t * (-8) *
(x0_t - x1_t) / (r2**(3/2))
v1h_t = v1h + d_t * (x0_t -
x1_t) / (r2**(3/2))
#
#
Update the values
# cf. Feynman et al. (1963),
p. 9-6.
#
x0 = x0_t
x1 = x1_t
v0h = v0h_t
v1h = v1h_t
#
#
Display the updated values
#
screen.fill((0, 0, 0))
pg.draw.circle(screen, (255,
200,255), (x_scale(x0[0]), y_scale(x0[1])), 5)
pg.draw.circle(screen, (0,
255, 255), (x_scale(x1[0]), y_scale(x1[1])), 20)
pg.display.flip()
for e in pg.event.get():
if
e.type == pg.QUIT:
# Is the windows's
close button clicked?
pg.quit()
sys.exit()
pg.time.wait(1)