🌓

Julia言語の常微分方程式ソルバーを使って円制限三体問題の軌道計算をしてみる〜イベント処理編〜

2023/01/14に公開

Julia

1.はじめに

本記事では，宇宙軌道力学分野でよく登場する円制限三体問題を題材にして，常微分方程式 (Ordinary Differential Equation; ODE) ライブラリDifferentialEquations.jlの拡張機能である「イベント処理」を紹介する．「イベント処理とは何ぞや？」という説明も後ほど行うが，軌道を計算する上で大活躍する機能である．DifferentialEquations.jlの入門（と境界値問題ライブラリBoundaryValueDiffEq.jlの使い方）に関しては，以下の記事にまとめている．

本記事で紹介する拡張機能を用いると，下記のような美しい図が描けるようになる（正確には下記の図を書くためには，並列計算機能も用いているが，その機能については，また別の記事で紹介する）．ちなみに，軌道設計の専門家は下記の図を用いて，「どうすれば効率よく月重力の影響を使いながら，軌道のエネルギーを増やしていくか？」ということを考えたりする $^1$ ．

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

2. 平面円制限三体問題の運動方程式

今回の例題では， $z$ 軸の運動を無視した 平面CRTBP（Planar Circular Restricted Three-Body Problem; PCRTBP） を用いる．PCRTBPで得られた軌道はCRTBPでも存在するが，CRTBPで存在する全ての軌道がPCRTBPで存在する訳ではないので，PCRTBPを扱う際には注意が必要である．座標系の取り方と無次元化の方法は通常のCRTBPと同様である．宇宙機の位置ベクトルを $\boldsymbol{r}=\begin{bmatrix}x & y \end{bmatrix}^{\top}$ ，速度ベクトルを $\boldsymbol{v}=\begin{bmatrix}v_x & v_y \end{bmatrix}^{\top}$ ，運動方程式を規定するパラメータを $\mu=\frac{m_2}{m_1+m_2}$ としたとき，CRTBPにおける宇宙機の運動方程式は，以下のような式で計算できる

\left\{ \begin{align*} \dot{\boldsymbol{r}} &= \boldsymbol{v}\\ \dot{\boldsymbol{v}} &= -\frac{1-\mu}{\|\boldsymbol{r}_1\|^3}\boldsymbol{r}_1 -\frac{\mu}{\|\boldsymbol{r}_2\|^3}\boldsymbol{r}_2 + \boldsymbol{a}_{\textrm{coli}} + \boldsymbol{a}_{\textrm{cf}} \end{align*} \right.

ただし，

\begin{align*} \boldsymbol{r}_1 &= \begin{bmatrix}x+\mu & y\end{bmatrix}^{\top}\\ \boldsymbol{r}_2 &= \begin{bmatrix}x-(1-\mu) & y\end{bmatrix}^{\top}\\ \boldsymbol{a}_{\textrm{coli}} &= \begin{bmatrix}2v_y & -2v_x\end{bmatrix}^{\top}\\ \boldsymbol{a}_{\textrm{cf}} &= \begin{bmatrix}x & y\end{bmatrix}^{\top} \end{align*}

である． $m_1$ と $m_2$ は，それぞれ第一天体と第二天体の質量を表すが，太陽＝地球系のCRTBPでは $\mu=3.036\times 10^{-6}$ ，地球＝月系のCRTBPでは $\mu=1.215\times 10^{-2}$ である．Juliaの関数で書くと，以下のようになる．

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

3. イベント処理機能の使い方

イベント処理 (Event Handler) 機能とは，ODEソルバーが数値解を計算する途中で"ある条件"を満たしたタイミングで"ある処理"を実行するような機能を指す．少し分かりづらいので，具体例を用いて，イベント処理機能が使われるシチュエーションを示そう．ODEソルバーは通常，（その途中結果がどうであれ）与えられた初期時刻から終端時刻まで計算し続ける．しかし，途中結果次第では，計算を停止したくなる状況は往々にある．例えば，「探査機が地球に突っ込んだタイミングでODEの計算を停止したい」とかである．この「〜したタイミング」というのは事前に知り得ないため，そのままでは地球に突っ込んだ後もODEソルバーが計算を続けてしまうことになる．イベント処理機能を使用すれば，探査機が地球に突っ込んだタイミングを検知して，ODEの計算を停止することができるようになる．その他にも，「衛星が近月点を通過するタイミング（＝条件）での位置・速度を保存したい（＝処理）」などにも対応が可能である．

DifferentialEquations.jlのイベント処理機能の仕様は，下記のウェブサイトに記載されている．本記事では，下記の仕様の説明を補いながら，イベント処理機能の説明を行う．

3.1. イベント処理機能の基本的な使い方

まずイベント処理機能には，イベント検知の"条件"を規定する関数condition()と条件を満たしたタイミングで実行される"処理"を規定する関数affect!()の２つを用意する必要がある．

condition()は，以下の様な形式で書かれたユーザー定義関数であり，戻り値がゼロになったタイミングでイベント処理機能が実行される．

function condition(u, t, integrator) # Event when event_f(u,t) == 0
    # ここに条件を記述（戻り値がゼロになったタイミングでイベント処理機能が実行される）
end

ここで，引数のuは状態量（＝従属変数），tは時刻（＝独立変数），integratorは積分に関わる情報を保存している構造体である（構造体の定義はリンク先を参照のこと）．

affect!()は，以下の様な形式で書かれたユーザー定義関数であり，条件が満たされたタイミングでここに記述した処理が実行される．

function affect!(integrator)
   # ここに処理を記述（条件が満たされたタイミングで，ここに記述した処理が実行される）
end

引数のintegratorを直接書き換えて，ODEの計算を行うことが可能であるため，Julia言語の慣例に則って，関数名の後ろに!を付けている（あくまで慣例だけなので!をつけたことによって処理結果が変わることはない）．ODE計算を停止したい場合は，

affect!(integrator) = terminate!(integrator)

と事前に用意されているterminate!()関数を使えば良い．

用意した関数をContinuousCallbackといったイベント処理を定義するための関数に渡して，インタフェースを定義します．それを以下のようにODEProblemのsolve()関数の引数として渡せば，イベント処理が実行可能です．

cb = ContinuousCallback(condition, affect!)

# ODEの問題定義・計算
prob = ODEProblem(eom_pcrtbp!, x0, tspan, μ)
sol = solve(prob, Vern7(), callback=cb)

例えば，PCRTBPの軌道伝播プログラムにおいて，第一天体あるいは第二天体に接近した場合に計算を停止するには，以下のようなプログラムを書けば良い．実行してみると，終端時刻3.0piより前に計算が停止していることが確認できる．

using DifferentialEquations

# 条件設定
μ = 0.01215058426994 # 地球・月系を仮定
x0 = [0.8369304414031459, -4.265566441114701, -4.076205405804159, -0.43142925920758124]
tspan = (0.0, 3.0pi)

# イベント処理機能の実装
function condition(u, t, integrator) # Event when event_f(u,t) == 0
    r1 = u[1:2] + [μ, 0.0]
    r2 = u[1:2] - [1.0 - μ, 0.0]
    return norm(r1) > 1.0e-3 && norm(r2) > 1.0e-3
end
affect!(integrator) = terminate!(integrator)
cb = ContinuousCallback(condition, affect!)

# ODEの問題定義・計算
prob = ODEProblem(eom_pcrtbp!, x0, tspan, μ)
sol = solve(prob, Vern7(), callback=cb, reltol=1e-10, abstol=1e-10)

イベント処理なしの場合（左）とイベント処理ありの場合（右）の結果の比較

3.2. イベント検知条件のゼロに近づく方向を限定する方法

イベント処理機能は，条件を規定する関数condition()がゼロになったタイミングで実行されるわけだが，関数がゼロに向かう方向（正からゼロに近づくのか？負からゼロに近づくのか？）を限定して実行することも可能である．使い所次第では，地味に使える機能である．

使い方はイベント処理のインタフェースを定義するContinuousCallback()関数を呼ぶときに，affect!だけではなく，第三引数にaffect_neg!を渡せば良い．

cb = ContinuousCallback(condition,affect!,affect_neg!)

第二引数に与えたaffect!は負からゼロに近づく（Upcrossing）ときに実行される処理を表し，第三引数に与えたaffect_neg!は正からゼロに近づく（Downcrossing）ときに実行される処理を表す．図で表すと下図のような感じである．

デフォルトでは，affect_neg!=affect!となっているため，第三引数を与えない状態では，両方とも検知され，同じ処理が実行される．片方のイベント検知を無視したい場合は，nothingを与えると良い．すなわち，以下のように書けば良い．

# 負からゼロに近づく（Upcrossing）ときのみ検知をしたい場合
cb_up = ContinuousCallback(condition,affect!,nothing)
#　正からゼロに近づく（Downcrossing）ときのみ検知をしたい場合
cb_down = ContinuousCallback(condition,nothing,affect_neg!)

軌道力学の問題で，この機能を上手く使うと，例えば，近地点を検出することが可能となる．条件を規定する関数condition()に，地球中心の軌道の半径方向速度を与える．近地点では半径方向速度がゼロになるわけだが，遠地点でも同様に半径方向速度がゼロになる．この２つを区別して，近地点を検出したい場合は，半径方向速度が負からゼロに近づく（Upcrossing）ときのみを検知すれば良い．文書のみで説明したが，具体的な実装を４章で示す．

3.3. イベント検知時の状態量を保存する方法

イベント検知をしたタイミングで，そのときの状態量uを保存してほしいようなシチュエーションがよくある．実は何もしなくても，デフォルトで，イベントを検知したタイミングの状態量uが，ODEの計算結果(＝solve()の戻り値)に保存されている．正確には，イベント検知時刻tの状態量uが２回保存される（イベント処理実行直前と直後の２回分）．これを変更したい場合は，イベント処理のインタフェースを定義するContinuousCallback()関数を呼ぶときに，引数save_positionsを渡せば良い．イベント処理実行直前のみ保存してほしい場合はsave_positions=(true,false)，直後のみ保存してほしい場合はsave_positions=(false,true)とすれば良い．

このままでは，イベント検知をしたタイミング以外の状態量も保存されてしまう．もしイベント検知をしたタイミングだけの状態量を抽出したい場合は，ODEProblemのsolve()関数に引数saveatを渡せば良い．そうすると，solve()の戻り値に，「saveatで指定された時刻の状態量」と「イベント検知時の状態量」が保存されるようになる．saveatにODEProblemの初期時刻と終端時刻を与えると，比較的簡単に「イベント検知時の状態量」が区別できる．（もっと良い方法があるかもしれないので，もし知っている方がいれば，是非コメント欄等で教えてほしい．）

例えば，イベント検知時に何もしなくて処理をしなくて良いので，ただ状態量を保存したい場合は，以下のような処理を実行して，sol[:,2:end-1]を取り出せば良い．

# condition関数は何かしら定義されているものとする
affect!(integrator) = nothing
cb = ContinuousCallback(condition, affect!, save_positions=(true, false))

# ODEの問題定義・計算
prob = ODEProblem(eom_pcrtbp!, x0, tspan, μ)
sol = solve(prob, Vern7(), callback=cb, saveat=[0.0, 3.0*pi])) # saveat=[tspanの最初, tspanの最後]

3.4. 複数のイベント処理を実装する方法

ODEを解くときに２つ以上のイベント処理を実装したい場合は，イベント処理のインタフェースを定義した上で，CallbackSet関数を使って，まとめれば良い．例えば，condition()で検知したときにaffect!()の処理を実行し，condition2()で検知したときにaffect2!()の処理を実行したい場合は，以下のように書けば良い．

cb = ContinuousCallback(condition, affect!)
cb2 = ContinuousCallback(condition2, affect2!)
cbs = CallbackSet(cb, cb2)

4. 平面円制限三体問題上の軌道の近点抽出

ここまでで説明したイベント処理機能を用いて，PCRTBPで伝播させた軌道から近地点の状態量を抽出するコードを以下に示す．近地点情報を抽出したい場合は，半径方向速度が負からゼロに近づく（＝Upcrossing）ときにイベント検出をするようにすれば良い．

using DifferentialEquations

# イベント処理#1: 半径方向速度から近地点を検知
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
affect!(integrator) = nothing # 保存をするだけなので，何も実行しない

# イベント処理#2: 地球 or 月に近づきすぎると計算を終了する（今回の例では使わないが，モンテカルロシミュレーションを解く際には有用）
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
affect2!(integrator) = terminate!(integrator) # DifferentialEquations.jlが用意しているODEの計算停止関数

# イベント処理のインタフェース定義
cb = ContinuousCallback(condition, affect!, nothing, save_positions=(true, false)) # 近地点(Upcrossing)を検出したい場合
# cb = ContinuousCallback(condition, nothing, affect!, save_positions=(true, false)) # 遠地点(Downcrossing)を検出したい場合
cb2 = ContinuousCallback(condition2, affect2!, save_positions=(false, false))
cbs = CallbackSet(cb, cb2)

# 計算条件の設定
μ = 0.01215058426994 # 地球・月系のμ
x0 = [1.1000137986067107, 0.0, 0.41017439921470306, -0.2165455493337665]
tspan = (0.0, 20.0*pi)

# 計算の実行
prob = ODEProblem(eom_pcrtbp!, x0, tspan, μ)
sol = solve(prob, Vern7(), callback=cbs, reltol=1e-10, abstol=1e-10, saveat=[0.0, 20.0*pi])

# 近地点の位置・速度は以下の値を参照すれば良い
sol[:, 2:end-1]

Makie.jlを用いて，抽出された近地点位置を描画するコードとその結果を以下に示す．

using CairoMakie
Makie.inline!(true)

# プロット用の軌道伝播
sol_forplot = solve(prob, Vern7(), callback=cbs, reltol=1e-10, abstol=1e-10)

# 描画
fig = Figure()
ax = Axis(fig[1, 1]; xlabel=L"x", ylabel=L"y")
scatter!(-μ, 0.0; markersize=20, color=:blue)
scatter!(1.0 - μ, 0.0; markersize=20, color=:gray)
lines!(sol_forplot[1, :], sol_forplot[2, :]; linewidth=1, color=:black)
scatter!(sol[1, 2:end-1], sol[2, 2:end-1]; markersize=20, color=:red)
ax.aspect = DataAspect()
fig

PCRTBP上の軌道の近地点抽出結果

5. おわりに

本記事では，宇宙軌道力学分野で用いられる円制限三体問題を題材として，常微分方程式の数値解析ライブラリDifferentialEquations.jlの拡張機能の１つであるイベント処理機能の使い方を紹介した．イベント処理機能とは，ODEソルバーが数値解を計算する途中で"ある条件"を満たしたタイミングで"ある処理"を実行するような機能を指す．このイベント処理機能を使って，円制限三体問題における近地点の検出を行う方法を示した．ここまで来れば，あとは計算時間を掛ければ，冒頭に示したような「近地点通過時の位相 $\theta$ と軌道長半径 $a$ をプロットした図」が書ける．普通にやると，かなりの計算時間が掛かるため，並列計算機能を使うことを推奨するが，この機能についても別の記事で改めて紹介したい．

参考文献

Shane D. Ross and Daniel J. Scheeres, "Multiple Gravity Assists, Capture, and Escape in the Restricted Three-Body Problem," SIAM Journal on Applied Dynamical Systems, Vol.6, Iss.3, pp.576-596, 2007. https://doi.org/10.1137/060663374