古典物理学と言えばニュートンの運動方程式で、何かの軌道をプロットするのが定番です。三体以上の問題は特殊なケースしか閉じた解が出てこないので、数値解析が行われています。銀河にある星の数をシミュレーションすると、スーパーコンピューターな数値演算能力が必要な世界になるわけですが、二個や三個であればPC上のRで計算できるので、試しにやってみましょう。
1. 計算する式
まずは物理の教科書に書いてありそうなニュートンの運動方程式の確認から*1。
が万有引力定数、が質量、が距離、は質点の数、は質点、が時点を表します。
力の成分をとに分解します。成分は省略。
両辺をで割って、、と置いて、常微分方程式に変形しましょう。
これをルンゲ=クッタ法などで計算します。二体問題でも解析解を出すのには色々と苦労するので、手抜きな感じがたまりません。
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. 計算結果
パラメーターに大きく依存するわけですが、上のコードだと下のようなグラフが描けます。最近は波ばかりプロットしていたので、グルグルしていて爽快ですね。
なおオイラー法などで計算すると、誤差が大きくなって全く違う図になります。