読者です 読者をやめる 読者になる 読者になる

餡子付゛録゛

ソフトウェア開発ツールの便利な使い方を紹介。

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が時点を表します。
力の成分をxyに分解します。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)
}

z軸を加えたり、質点の数を増やすのは難しく無いと思います。例えば以下のように三つのパラメーターをいじれば、三体問題化できます。

n <- 3
m <- c(500, 2, 1)
X <- c(100, 400, 800, 100, 400, 0, 1, -1, 0, 1, 1, 0)

衝突処理が無いので質点が近づきすぎると、どこかに飛んで行くのでパラメーターの設定には注意してください。

3. 計算結果

パラメーターに大きく依存するわけですが、上のコードだと下のようなグラフが描けます。最近は波ばかりプロットしていたので、グルグルしていて爽快ですね。

なおオイラー法などで計算すると、誤差が大きくなって全く違う図になります。

*1:詳しい説明は富久(2010)の第5章を参照