🌔

Julia言語の常微分方程式ソルバーを使って円制限三体問題の軌道計算をしてみる〜並列計算/モンテカルロ・シミュレーション編〜

2023/02/03に公開

Julia

1.はじめに

本記事では，宇宙軌道力学分野で使われる円制限三体問題を題材にして，常微分方程式 (Ordinary Differential Equation; ODE) ライブラリDifferentialEquations.jlの使い方を説明するシリーズ記事の第３弾である．第１弾ではDifferentialEquations.jlの入門と境界値問題ライブラリBoundaryValueDiffEq.jlの使い方を紹介した．

第２弾では，ODEの計算途中に"ある条件"を満たしたタイミングで"ある処理"を実行するようなイベント処理（Event Handler）機能の使い方を紹介した．

第３弾では，DifferentialEquations.jlの拡張機能である「並列計算によるモンテカルロ・シミュレーション」を紹介する．「モンテカルロ・シミュレーションとは何ぞや？」という説明も後ほど行うが，この機能も軌道計算をする上で大活躍する機能である．

本記事で紹介する拡張機能を用いると，いよいよ下記のような美しい図が描けるようになる．

地球・月系の三体問題で近地点通過時の位相 $\theta$ と軌道長半径 $a$ をプロットした図

今回も第２弾の記事で扱った平面円制限三体問題の運動方程式を用いる．

using LinearAlgebra

# Definition of Equation of Motion of Planar CRTBP
function eom_pcrtbp!(dx, x, μ, t)
    # Preparation
    r1 = x[1:2] + [μ, 0.0]
    r2 = x[1:2] - [1.0 - μ, 0.0]
    a_cori = [2.0x[4], -2.0x[3]]
    a_cf = [x[1], x[2]]
    # dx/dt
    dx[1:2] = x[3:4]
    dx[3:4] = -(1.0 - μ) * r1 / (norm(r1)^3) - μ * r2 / (norm(r2)^3) + a_cori + a_cf
end

2. 並列計算によるモンテカルロ・シミュレーション

モンテカルロ・シミュレーション (Monte Carlo Simulations) とは，ランダムに初期条件やパラメータを変えながら複数回シミュレーションをして，その統計的結果を得ることを指す．DifferentialEquations.jlでは，並列計算によるモンテカルロ・シミュレーションを簡単に実装するアンサンブル・シミュレーション (Ensemble Simulations)機能が用意されている．この機能を用いると，モンテカルロ・シミュレーション以外にも，グリッドサーチ (Grid Search) のように事前に用意したデータセットを用いて条件を変えながら，複数回シミュレーションを行うことも可能である．

2.1. アンサンブル・シミュレーションの使い方

アンサンブル・シミュレーションは，

ODEProblemを定義して解きたいODEの問題設定を行う
試行ごとに実施したい処理（例えば，ODEの初期値の変更等）をprob_func()という関数の形で定義する
1.と2.を引数として，EnsembleProblemを定義する
solve()を実行して，アンサンブルシミュレーションを解く（このとき，試行回数をtrajectoriesという形で定義する）

という流れで実行できる．Julia言語を用いた具体的な実装方法を説明しよう．

まずODEProblemを定義する．ここで，パラメータ・初期条件・伝播時間等は，後ほどprob_func()の中で変更可能なので，ひとまず仮置きする値を入れる．

using DifferentialEquations

# ODEの条件・問題設定（後でアンサンブル・シミュレーションを行う際に値を変更できるので，ひとまず仮置き）
μ = 0.01215058426994 # 地球・月系を仮定
x0 = [0.8369304414031459, -4.265566441114701, -4.076205405804159, -0.43142925920758124]
tspan = (0.0, 3.0pi)
prob = ODEProblem(eom_pcrtbp!, x0, tspan, μ)

次に試行ごとに実施したい処理をprob_func()という関数系で定義する．例えば，位置ベクトルの $x$ 成分をランダムに変更しながら，モンテカルロ・シミュレーションを実施したい場合は，以下のようなprob_func()を定義すれば良い．初期条件をランダムに変更するのではなく，予め定義した初期条件を使うことも可能である（参考URL）．この機能を使うことで，いわゆるグリッドサーチを走らせることができる．

function prob_func(prob, i, repeat)
    x0 = [rand(), -4.265566441114701, -4.076205405804159, -0.43142925920758124]
    remake(prob, u0=x0)
end

ここで，remake()という関数は，与えられたODEProblemの問題設定に修正を加える関数であり，上記の例では初期条件u0を新たなx0で置き換えることを意味している．

あとはEnsembleProblemを定義して，解くだけである．試行回数を10回としてモンテカルロ・シミュレーションを解く場合は，以下のように書けば良い．

ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
sim = solve(ensemble_prob,Vern7(),EnsembleDistributed(),trajectories=10)

solve()で問題を解くときにtrajectories=10を与えることで，試行回数10回のモンテカルロ・シミュレーションが可能である．

このくらいの機能であれば，forループを用いて自作で実装しても簡単な気がするが，EnsembleProblemにはシミュレーション実行結果に対して統計的な処理を行う機能（EnsembleSummary）もあり，こういう機能を使い出すと便利さを実感しそうである．

2.2. 並列計算方法

ここまでで説明した方法は並列計算を用いないシンプルなアンサンブル・シミュレーションだった．並列計算を用いたい場合は，Julia言語の標準ライブラリであるDistributedモジュールを用いることで簡単に実装できる．Julia言語のマルチプロセス処理と分散計算のマニュアルに詳しい説明が書かれているが，ここでは詳細な説明には踏み込まない．並列計算と言っても，複数のコンピュータを用いる方法，１つのコンピュータ内の複数コアを用いる方法，GPUを用いる方法等があるが，ここでは1つのコンピュータ内の複数コアを用いる方法を紹介する．

複数コアを用いたマルチプロセスを立ち上げるためには，以下のプログラムを走らせる（参考ページ）

using Distributed
addprocs(exeflags="--project=$(Base.active_project())")

ここで，addprocs()の引数にexeflags="--project=$(Base.active_project())"を与えることで，現在のJuliaプロジェクト環境を複数プロセスに渡すことができる．

マルチプロセスを立ち上げた後は，並列計算で利用するモジュールや関数・グローバル変数を定義する際に，@everywhereマクロを宣言の頭に付けることで，全てのプロセスにこれらの情報を渡すことができ，並列計算が実行できる．

# 例①：全てのプロセスでモジュールを利用
@everywhere using DifferentialEquations
@everywhere using LinearAlgebra

# 例②：全てのプロセスで関数を定義
@everywhere function prob_func(prob, i, repeat)
    x0 = [rand(), -4.265566441114701, -4.076205405804159, -0.43142925920758124]
    remake(prob, u0=x0)
end

EnsembleProblemに対する具体的な並列計算の実装結果を第4章に示す．

3. 円制限三体問題の軌道力学に関する補足情報

最初に紹介した近地点通過時の位相 $\theta$ と軌道長半径 $a$ をプロットした図が描きたい人向けに少し情報を補足する．（軌道力学に興味はなく）並列計算・モンテカルロシミュレーションにのみ興味がある読者は読み飛ばしていただきたい．

3.1. ヤコビ定数

力学の問題において，運動量や力学的エネルギーという保存量があるように，円制限三体問題において，宇宙機（第三天体）の運動を規定するヤコビ定数 (Jacobi Constant) という保存量が存在する．言い換えると，円制限三体問題で伝播された軌道に沿ってヤコビ定数を評価すると，同じ値（＝値が保存されている）になる．導出は割愛するが，ヤコビ定数の定義は以下のようになる．

\begin{align} C = (x^2 + y^2) + \frac{2(1-\mu)}{r_1} + \frac{2\mu}{r_2} - (v_x^2+v_y^2+v_z^2) \end{align}

教材によって，ヤコビ定数の定義には，若干の差がある（上式を２倍しているものや上式の正負を反転しているもの）ため，論文を読む際にはどのように定義されているのかを確認しながら，読む必要がある．

3.2. 近地点通過時の位相 $\theta$ の計算

CRTBPで伝播された軌道において，近地点通過時の状態量（位置・速度）を取得するためには，イベント処理機能を用いる．

イベント機能を用いて近地点通過時の状態量（位置・速度）を取得すると，（第一天体から見た）回転座標系での位相角 $\theta$ は以下の式で計算できる．

\begin{align} \theta = \tan^{-1}\left(\frac{y}{x+\mu}\right) \end{align}

$\tan^{-1}$ を計算する際には，象限の扱いに気を付ける必要がある．C言語等では $[-\pi, \pi]$ の値域を持つatan2()という関数が用意されているが，Julia言語では，atan()関数に与える引数の数で自動判別される（Julia言語のatan関数）．xとyに近地点通過時の位置がArray{Float64}形式で格納されている場合， $\theta$ を計算するコードは以下のようになる．

θ = atan.(y, x.+μ)

ここで，.（ドット）が２箇所に打たれているが，Julia言語ではドットを付けることによって，配列の成分ごとの計算が可能である．例えば，x.+μとすることで，配列xのそれぞれの成分にμを加える演算が可能である．更に，atan.()とすることで，関数の引数が配列であっても，それぞれの成分に対して，関数の計算をしてくれる．

3.3. 近地点通過時の軌道長半径の計算

軌道長半径 $a$ を計算するためには，軌道力学におけるエネルギー保存則（Vis-viva方程式）

\begin{align} \frac{v^2}{2} - \frac{GM}{r} = - \frac{GM}{2a} \end{align}

を用いる．今回は，第一天体に対する軌道長半径を計算したいので，

\begin{align} \frac{v_1^2}{2} - \frac{1-\mu}{r_1} = - \frac{1-\mu}{2a} \end{align}

を $a$ について解くことで，軌道長半径を得ることができる．

\begin{align} a = \left(\frac{2}{r_1} - \frac{v_1^2}{1-\mu}\right)^{-1} \end{align}

ここで， $r_1$ と $v_1$ は第一天体に対する慣性座標系での位置・速度であり，

\begin{align} r_1 &= \sqrt{\left(x+\mu\right)^2 + y^2}\\ v_1 &= \sqrt{\left(v_x - y\right)^2 + \left(v_y + x + \mu\right)^2} \end{align}

と計算できる．CRTBPで扱っている状態量は回転座標系であるのに対して，Vis-viva方程式で用いる状態量は慣性座標系における状態量であるため，特に $v_1$ を計算する際には回転成分を補正する必要がある．

4. 近地点マップの描画プログラム

ここまでで説明したアンサンブル・シミュレーション機能（とイベント処理機能）を用いて，PCRTBPで伝播させた軌道の近地点を検出し， $\theta-a$ を描画するコードを以下に示す．数分くらいの計算時間で終わるはずである．冒頭で紹介した図はnum_trj = 50000で解いた結果（この場合，100分前後の計算時間が掛かる）であり，その実装は著者のGithub上のJupyterノートブックに記載している．

## 0. モジュールのインポート（並列計算の初期化）
using Distributed
addprocs(exeflags="--project=$(Base.active_project())")
@everywhere using DifferentialEquations
@everywhere using LinearAlgebra

## 1. ODEの問題設定
# 計算条件の定義
num_trj = 1000
@everywhere μ = 0.01215058426994 # 地球・月系のμ
@everywhere t_fin = 2.0pi * 100
x0_temp = [0.8369304414031459, -4.265566441114701, -4.076205405804159, -0.43142925920758124]
tspan = (0.0, t_fin)

# 平面円制限三体問題の運動方程式の定義
@everywhere function eom_pcrtbp!(dx, x, μ, t)
    # Preparation
    r1 = x[1:2] + [μ, 0.0]
    r2 = x[1:2] - [1.0 - μ, 0.0]
    a_cori = [2.0x[4], -2.0x[3]]
    a_cf = [x[1], x[2]]
    # dx/dt
    dx[1:2] = x[3:4]
    dx[3:4] = -(1.0 - μ) * r1 / (norm(r1)^3) - μ * r2 / (norm(r2)^3) + a_cori + a_cf
end

## 2. イベント処理の定義
# イベント処理#1: 半径方向速度から近地点を検知
@everywhere function condition(x, t, integrator) # Event when event_f(u,t) == 0
    r1_vec = x[1:2] + [μ, 0.0]
    return r1_vec ⋅ x[3:4] # 半径方向速度（＝地球からみた探査機の位置ベクトルと速度ベクトルの内積）
end
@everywhere affect!(integrator) = nothing # 保存をするだけなので，何も実行しない

# イベント処理#2: 地球 or 月に近づきすぎると計算を終了する（今回の例では使わないが，モンテカルロシミュレーションを解く際には有用）
@everywhere function condition2(x, t, integrator) # Event when event_f(u,t) == 0
    r1 = x[1:2] + [μ, 0.0]
    r2 = x[1:2] - [1.0 - μ, 0.0]
    return norm(r1) > 1.0e-3 && norm(r2) > 1.0e-3
end
@everywhere affect2!(integrator) = terminate!(integrator) # DifferentialEquations.jlが用意しているODEの計算停止関数

# イベント処理のインタフェース定義
cb = ContinuousCallback(condition, affect!, nothing, save_positions=(true, false))
cb2 = ContinuousCallback(condition2, affect2!, save_positions=(false, false))
cbs = CallbackSet(cb, cb2)

## 3. アンサンブル・シミュレーションの定義と実行
# 試行ごとに実行したい処理の定義
@everywhere function prob_func(prob, i, repeat)
    jacobi_constant = 3.15 # ヤコビ定数
    x0 = 3.0rand() # 初期位置のx成分を0〜3までの範囲でランダムに定義
    φ0 = 2.0pi*rand()
    v0 = sqrt(x0^2 + 2.0 * (1.0 - μ) / abs(x0 + μ) + 2.0 * μ / abs(x0 - 1.0 + μ) - jacobi_constant) # ヤコビ定数が一致するように速度を定義
    u0 = [x0, 0.0, v0 * cos(φ0), v0 * sin(φ0)]
    remake(prob, u0=u0) # 初期値を試行ごとにランダムに変更
end

# アンサンブルシミュレーションを定義し，問題を解く
prob = ODEProblem(eom_pcrtbp!, x0_temp, tspan, μ)
ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
sim = solve(ensemble_prob, Vern7(), EnsembleDistributed(), trajectories=num_trj, callback=cbs, saveat=t_fin);

Makie.jlを用いて，抽出された近地点通過時の $\theta-a$ を描画するコードとその結果を以下に示す．

## 4. Makieによる結果の描画
using CairoMakie

fig = Figure(; resolution=(1024, 1024))
ax = Axis(fig[1, 1])

for n = 1:num_trj
    x_all = hcat(sim[n].u...)'
    x_all = x_all[2:end-1, :] # 初期値と終端値を取り除く

    # θと軌道長半径の計算
    θ = atan.(x_all[:, 2], x_all[:, 1] .+ μ)
    r_1 = sqrt.((x_all[:, 1] .+ μ) .^ 2 + x_all[:, 2] .^ 2)
    v_1_sq = (x_all[:, 3] - x_all[:, 2]) .^ 2 + (x_all[:, 4] + x_all[:, 1] .+ μ) .^ 2
    a = 1.0 ./ (2.0 ./ r_1 - v_1_sq ./ (1.0 - μ))

    scatter!(mod2pi.(θ), a, color=:black, markersize=1.5)
end

xlims!(0.0, 2.0pi)
ylims!(0.2, 0.8)
fig

5. おわりに

本記事では，宇宙軌道力学分野で用いられる円制限三体問題を題材として，常微分方程式の数値解析ライブラリDifferentialEquations.jlの拡張機能の１つである並列計算/モンテカルロ・シミュレーション機能の使い方を紹介した．Julia言語の並列計算やモンテカルロ・シミュレーション機能がいかにパワフルであるか感じとってもらえただろうか．Julia言語の常微分方程式ソルバーDifferentialEquations.jlは，非常に効率の良いライブラリなので，宇宙軌道力学分野はもちろんのこと幅広い数値解析分野で幅広く活用できるだろう