§1. はじめに
こんにちは、@Inuverse44 です。
ソフトウェア開発をはじめて1年と半年程度経ちました。大学院時代は宇宙論を研究していたのですが、今回記事にしたのは、過去の自分に向けて、当時知っておけばもっと楽に研究の数値計算ができた方法についてです。
宇宙論という分野でもっとも一般に聞き馴染みのある言葉はおそらく「ビッグバン」です。実際は時空のダイナミクスや、重力の拡張、銀河構造形成など、対象は様々なわけですが、どの対象にしても往々にして「モデル」をいろいろ切り替えながら議論するわけです。
同じようなロジックで計算をするけれども、いろいろなモデルを切り替えながら計算したい。私はそれだけのためにファイルを作成したり関数を作成していました(そのような煩雑な方法しか知らなかったわけです)。今回はその当時の悩みを解消できる実装方法の「ストラテジパターン」を用いて、宇宙論の基礎的な計算である宇宙年齢を、いくつかのモデルに対して計算してみようと思います。
§2. ストラテジパターン
そもそもソフトウェア開発において、デザインパターンという概念があります。これは、長年の経験則から、いくつかの便利な実装パターンのことを指します。
ストラテジパターンは、いくつか似たような計算があるときに便利です。
詳しくは後述しますが、今回の例では、
- 物質優勢宇宙での宇宙年齢
- 放射優勢宇宙での宇宙年齢
- $\Lambda$ CDMモデルでの宇宙年齢
- ...
のように、計算ロジックはほぼ同じだけれども、設定が異なる時に威力を発揮するのがストラテジパターンです。
より詳しいことについては、下記記事などを参考にしてみてください。
今回の例では§3で解説されています。
§3. 宇宙年齢の計算
宇宙論におけるもっとも基礎的だといってよい式であるフリードマン方程式からスタートします:
H^2
= H_0^2
\left(
\frac{\Omega_\mathrm{m, 0}}{a^3}
+ \frac{\Omega_\mathrm{r, 0}}{a^4}
+ \Omega_\mathrm{\Lambda, 0}
+ \frac{\Omega_\mathrm{k, 0}}{a^2}
\right)
\tag{1}
詳細は解説しませんが、
- $a = a(t)$:スケール因子(ざっくりと宇宙の大きさ)
- $H(t) = \frac{\mathrm{d}a/\mathrm{d}t}{a} = \dot{a}/a$:ハッブルパラメタ(ざっくりと宇宙の膨張率)
- $H_0 = H(t_0)$:現在のハッブルパラメタ(ざっくりと宇宙の膨張率)
- $\Omega_\mathrm{m, 0}$:現在の物質に密度に関するパラメタ
- $\Omega_\mathrm{r, 0}$:現在の放射エネルギーに関するパラメタ
- $\Omega_\mathrm{r, 0}$:現在の真空エネルギーに関するパラメタ
- $\Omega_\mathrm{r, 0}$:現在の宇宙の曲率に関するパラメタ
です。
さて、(1)式を用いると、宇宙年齢は、
\begin{align}
t_\mathrm{age}
&=
\int_0^{t_\mathrm{age}} \mathrm{d}t
= \int_0^1 \frac{\mathrm{d}a}{\dot{a}}
= \int_0^1 \frac{\mathrm{d}a}{aH}
= \frac{1}{H_0} \int_0^1 \frac{\mathrm{d}a}{a \frac{H}{H_0}}
\\
&=
\frac{1}{H_0} \int_0^1 \frac{\mathrm{d}a}{a E}
\tag{2}
\\
&=
\frac{1}{H_0}
\int_0^1
\frac{\mathrm{d}a}{
\left[
\frac{\Omega_\mathrm{m,0}}{a}
+ \frac{\Omega_\mathrm{r, 0}}{a^2}
+ \Omega_\mathrm{\Lambda, 0} a^2
+ (
1
- \Omega_\mathrm{m, 0}
- \Omega_\mathrm{r, 0}
- \Omega_\mathrm{\Lambda, 0}
)
\right]^{1/2}
}
\tag{3}
\end{align}
となります。途中で導入した$E$は無次元パラメタで$E = H/H_0$です。明示的に書けば、
E(a)
= \left[
\frac{\Omega_\mathrm{m,0}}{a^3}
+ \frac{\Omega_\mathrm{r, 0}}{a^4}
+ \Omega_\mathrm{\Lambda, 0}
+ \frac{(
1
- \Omega_\mathrm{m, 0}
- \Omega_\mathrm{r, 0}
- \Omega_\mathrm{\Lambda, 0}
)}{a^2}
\right]^{1/2}
\tag{4}
です。
※ここでは、$a(0) = 0$とします。
詳しく知りたい場合は例えば、
などを参考にしてください。
§4. 実装
宇宙年齢を計算する基礎的な問題では、先の節に登場した(3)式の$\Omega_\mathrm{m, 0}$のみnonzeroであったり、$\Omega_\mathrm{r, 0}$のみnonzeroのような場合を計算しますし、一般の場合では全てのパラメタに具体的な数字を入れて数値計算します。
しかしながら、それぞれのパターンでロジックを書くのは避けたいです。我々が結局解きたいロジックというのは、より抽象的な計算の
\frac{1}{H_0} \int_0^1 \frac{\mathrm{d}a}{a E}
の部分で、具体的な値を代入した$\Omega_\mathrm{m, 0}$や$\Omega_\mathrm{r, 0}$は、色々のパターンで付け替えできるような状態であったほうが嬉しいです。
ストラテジパターンはこの抽象の部分と具体的な付け替えできる部分を同時に提供してくれる書き方になります。
抽象的な方から見てみましょう。一般の数値計算をする人たちはCやC++、Fortranが主流でしょうが、業務で使用している言語がKotlinのため、数値計算の速さ度外視のKotlinで実装しています。
◼︎抽象的なクラス
interface CosmologyModelStrategy {
val name: String
/**
* Calculates the dimensionless Hubble parameter E(a) = H(a) / H0.
*/
fun E(a: Double): Double
/**
* Returns the integrand for the age calculation: 1 / (a * E(a)).
* This corresponds to dt * H0 = da / (a * E(a)).
*/
fun integrand(a: Double): Double {
if (a <= 0.0) return 0.0 // Singularity at a=0 handled: lim a->0 is 0 for standard models
val e = E(a)
return if (e == 0.0) 0.0 else 1.0 / (a * e)
}
}
interfaceでfun E(a: Double): Doubleとしているので、これを継承するクラスは、この関数に具体的な実装を施すことを強要されます。
次に取り替えられる2つの具体的なモデルを用意しておきましょう。数式とそれに対応するコードを示しておきます。
◼︎具体的なモデル1(物質優勢宇宙)
物質優勢宇宙の時、$\Omega_\mathrm{m, 0} = 1$であり、それ以外は(近似的に)$\Omega_\mathrm{r, 0}, \Omega_\mathrm{\Lambda, 0}, \Omega_\mathrm{k, 0} = 0$です。よって、$E$は
E(a) = \frac{1}{\sqrt{a^{3}}}
です。また、計算しなければならない積分は
t_\mathrm{age}
= \frac{1}{H_0}
\int_0^1 \sqrt{a} \, \mathrm{d}a
となります。これは解析的に計算すると$t_\mathrm{age} = \frac{2}{3} H_0^{-1}$となりますが、これは数値積分との整合性の確認のために利用しましょう。
コードは下記に対応します。
class MatterDominatedFlatImpl : CosmologyModelStrategy {
override val name = "Matter Dominated Flat"
override fun E(a: Double): Double {
if (a <= 0.0) return Double.MAX_VALUE
return sqrt(a.pow(-3))
}
}
◼︎具体的なモデル2($\Lambda$CDMモデル)
もう一つのモデルは138億年という数字を出してくれる有名なモデルです。先とは違い、$\Omega_\mathrm{i, 0}$の値は観測で決められ、全てnonzeroとします。つまり、とくべき式は(4)を代入して(3)の数値積分なのです。
諸々の値は
の結果を使用しています。
class LambdaCdmImpl(
private val Om: Double = CosmologyConstants.OMEGA_MATTER_NOW,
private val Or: Double = CosmologyConstants.OMEGA_RADIATION_NOW,
private val Ol: Double = CosmologyConstants.OMEGA_VACUUM_NOW
) : CosmologyModelStrategy {
private val Ok = 1.0 - Om - Or - Ol
override val name = "LambdaCDM (Om=$Om, Or=$Or, Ol=$Ol)"
override fun E(a: Double): Double {
if (a <= 0.0) return Double.MAX_VALUE // Effectively infinite H at big bang
val term1 = Om * a.pow(-3)
val term2 = Or * a.pow(-4)
val term3 = Ol
val term4 = Ok * a.pow(-2)
return sqrt(term1 + term2 + term3 + term4)
}
}
これらを計算するためには、どの具体的なモデルを使うのかを選択するためのクラスを設定する必要があります。
class Calculator(
private var strategy: IntegrateStrategy, // 積分モデルのインスタンス
var model: CosmologyModelStrategy // 宇宙論モデルのインスタンス
) {
// 積分を実行する
fun run(start: Double, end: Double, eps: Double = 1e-6): Double {
// ↓ここが非積分関数
return strategy.run(start, end, eps) { a -> model.integrand(a) }
}
// 宇宙年齢
fun calculateAge(hubbleTimeGyr: Double = CosmologyConstants.HUBBLE_TIME_GYR): Double {
// 積分して宇宙年齢を計算
val integral = run(0.0, 1.0)
return hubbleTimeGyr * integral
}
}
このコードでは、数値積分もストラテジパターンで書いてあります。このCalculatorクラスに自分の計算したいモデルのインスタンスを入れることで、簡単に計算を走らせることができます。
もし、さらに計算したいモデル(例えば放射優勢宇宙の場合など)があれば、定義したinterfaceに倣って、具体的なモデルを追加で書き下せばOKです。
main関数では下記のような形でモデルを交換しています。ただただ、calculator.modelに自分の計算したいモデルのインスタンスを事前に代入しておくことで、計算が実現可能です。
val lambdaCdm = LambdaCdmImpl()
val calculator = Calculator(SimpsonImpl(), lambdaCdm)
println("Model: ${lambdaCdm.name}")
println("H0: $HUBBLE_NOW km/s/Mpc")
var age = calculator.calculateAge()
println("Age: $age Gyr")
println()
val edsModel = MatterDominatedFlatImpl()
calculator.model = edsModel
println("Model: ${edsModel.name}")
age = calculator.calculateAge()
println("Age: $age Gyr")
val hubbleTime = HUBBLE_TIME_GYR
val theoreticalEdsAge = (2.0 / 3.0) * hubbleTime
println("Theoretical Age (2/3 * 1/H0): $theoreticalEdsAge Gyr")
println("Difference: ${kotlin.math.abs(age - theoreticalEdsAge)} Gyr")
最終的にこのコードを実行すると
と物質優勢宇宙の場合には理論値との絶対誤差が$O(10^{-6})$で、$\Lambda$CDMモデルの宇宙時間がほぼ想定値の 13.8 Gyr、つまり138億年が計算できました。
§5. まとめ
モデルを複数使うような系を頻繁に使うような分野の場合、数値計算にストラテジパターンを実装しておくと、数値計算活動がよりスムーズに進むのではないかと思います。やっていることは、(2)式のように、抽象的な計算ロジックを定義すること(interface)と具体的なモデルを複数定義したこと、そしてどのモデルを使っているのかをコンテキストを記録するクラスを実装しています。
ソースコード
