Rで多体問題を数値解析 - 餡子付゛録゛

古典物理学と言えばニュートンの運動方程式で、何かの軌道をプロットするのが定番です。三体以上の問題は特殊なケースしか閉じた解が出てこないので、数値解析が行われています。銀河にある星の数をシミュレーションすると、スーパーコンピューターな数値演算能力が必要な世界になるわけですが、二個や三個であればPC上のRで計算できるので、試しにやってみましょう。

1. 計算する式

まずは物理の教科書に書いてありそうなニュートンの運動方程式の確認から*1。

$m_i \frac{d^2r_i}{dt^2}=-\sum^n_{j=1,j \neq i}G\frac{m_i m_j}{r^2_{ij}}$

$G$ が万有引力定数、 $m$ が質量、 $r$ が距離、 $n$ は質点の数、 $i$ は質点、 $t$ が時点を表します。
力の成分を $x$ と $y$ に分解します。 $z$ 成分は省略。

$m_i \frac{d^2x_i}{dt^2}=-\sum^n_{j=1,j \neq i}G\frac{m_i m_j}{r^2_{ij}} \cdot \frac{x_j-x_i}{r_{ij}}$

$m_i \frac{d^2y_i}{dt^2}=-\sum^n_{j=1,j \neq i}G\frac{m_i m_j}{r^2_{ij}} \cdot \frac{y_j-y_i}{r_{ij}}$

$r_{ij} = \sqrt{(x_j-x_i)^2 + (y_j-y_i)^2}$

両辺を $m_i$ で割って、 $x2=\frac{d^2x_i}{dt^2}$ 、 $y2=\frac{d^2y_i}{dt^2}$ と置いて、常微分方程式に変形しましょう。

$x_1=\frac{dx_i}{dt}$
$y_1=\frac{dy_i}{dt}$
$x_2=-\sum^n_{j=1,j \neq i}G\frac{m_j}{r^2_{ij}} \cdot \frac{x_j-x_i}{r_{ij}}$
$y_2=-\sum^n_{j=1,j \neq i}G\frac{m_i m_j}{r^2_{ij}} \cdot \frac{y_j-y_i}{r_{ij}}$

これをルンゲ＝クッタ法などで計算します。二体問題でも解析解を出すのには色々と苦労するので、手抜きな感じがたまりません。

2. ソースコード

N=3でも4でもいいのですが、N=2のソースコードを提示します。まずは演算部分。

# deSolveパッケージを使います
library("deSolve")

# 常微分方程式化したモデル
model <- function(t, X, params){
  with(as.list(params),{
    # 引数を分解
    x <- X[1:n]
    y <- X[(n+1):(2*n)]
    dx1 <- X[(2*n+1):(3*n)]
    dy1 <- X[(3*n+1):(4*n)]
    # 戻り値の領域を確保
    dx2 <- numeric(n)
    dy2 <- numeric(n)
    for(i in 1:n){
      r3 <- sqrt((x[-i]-x[i])^2+(y[-i]-y[i])^2)^3
      dx2[i] <- sum(-G*m[-i]*(x[-i]-x[i])/r3)
      dy2[i] <- sum(-G*m[-i]*(y[-i]-y[i])/r3)
    }
    list(c(dx1, dy1, dx2, dy2))
  })
}

# 計算のループ条件。0.1刻みで250まで。
times <- seq(0, 250, 0.1)
# 質点の数
n <- 2
# 質点の重さ
m <- c(500, 1)
# 初期座標と初期速度を以下の並び順で指定
# x[1], x[2], y[1], y[2], dx[1], dx[2], dy[1], dy[2]
X <- c(100, 400, 100, 400, 1, -1, 1, 0)
# 常微分方程式として計算する
out1 <- ode(X, times, model, parms = list(n=n, G=-9.80665, m=m), method="rk4")

演算結果は数値の羅列で分かりづらいので、グラフにしましょう。

# プロット範囲を指定
xlim <- c(0, 800)
ylim <- c(0, 800)

# N時点までのプロット
plot_as_of_N <- function(N){
  plot(out1[,2][1:N], out1[,2+n][1:N], xlim=xlim, ylim=ylim, type="l", xlab="x", ylab="y", lty=2)
  for(i in 2:n){
    par(new=T)
    plot(out1[,i+1][1:N], out1[,i+1+n][1:N], xlim=xlim, ylim=ylim, type="l", lty=i+1, xlab="", ylab="", main=paste("t=", round(N)))
  }
  for(i in 1:n){
    points(out1[N,i+1], out1[N,i+1+n], pch=i, cex=2)
  }
  legend("topleft", lty=seq(2, length(m)+1), legend=sprintf("m=%d", m), bty="n")
}

# 同時に2×2個をプロット
par(mfrow = c(2,2), oma = c(0, 0, 0, 0))
for(i in seq(1, length(out1[,1]), (length(out1[,1])-1)/3)){
  plot_as_of_N(i)
}