Return to 宇宙論パラメータと宇宙年齢

補足：SymPy と SPB で宇宙年齢のグラフを描く

宇宙年齢の表式の導出については，以下のページを参照。

ここでは，Python の SymPy Plotting Backends を使って宇宙年齢のグラフを描いてみる。

In [1]:
from sympy import *
from sympy.abc import *
# 虚数単位，円周率，ネイピア数
from sympy import I, pi, E
from spb import *

# グラフを SVG で Notebook にインライン表示
%config InlineBackend.figure_formats = ['svg']
In [2]:
# グラフを描くためではなくデフォルト設定のため
import matplotlib.pyplot as plt

config = {
    'axes.labelsize': 'x-large',
    'mathtext.fontset': 'cm'
}
plt.rcParams.update(config)

$\Omega_{\Lambda} = 0$ の場合

$$H_0 t_0 = -\frac{1}{\Omega_{\rm m} -1}+\frac{\Omega_{\rm m}}{(\Omega_{\rm m}-1)^{\frac{3}{2}} }
\tan^{-1}\sqrt{\Omega_{\rm m}-1} \quad \mbox{for}\ \ \Omega_{\rm m} > 1$$$$H_0 t_0 = \frac{1}{1-\Omega_{\rm m}}-\frac{\Omega_{\rm m}}{(1-\Omega_{\rm m})^{\frac{3}{2}} }
\tanh^{-1}\sqrt{1-\Omega_{\rm m}} \quad \mbox{for}\ \ \Omega_{\rm m} < 1$$

$\Omega_{\rm m} \rightarrow \Omega$ として…

In [3]:
# Omega > 1
def t1(x):
    return -1/(x-1) + x/((x-1)*sqrt(x-1)) * atan(sqrt(x-1))
# Omega < 1
def t2(x):
    return 1/(1-x) - x/((1-x)*sqrt(1-x)) * atanh(sqrt(1-x))
# Omega = 1
def t0(Omega):
    return 2/3

def t(Omega):
    if Omega > 1:
        return t1(Omega)
    elif Omega == 1:
        return t0(Omega)
    else:
        return t2(Omega)

$\Omega_{\Lambda} = 1 – \Omega_{\rm m}$ すなわち $k = 0$ の場合

$$H_0 t_0 = \frac{2}{3\sqrt{\Omega_{\rm m} -1}}\tan^{-1} \sqrt{\Omega_{\rm m} -1} \quad \mbox{for}\ \ \Omega_{\rm m} > 1$$$$H_0 t_0 = \frac{2}{3\sqrt{1-\Omega_{\rm m} }}\tanh^{-1} \sqrt{1-\Omega_{\rm m} } \quad \mbox{for}\ \ \Omega_{\rm m} < 1$$

$\Omega_{\rm m} \rightarrow \Omega$ として…

In [4]:
# Omega > 1
def T1(Omega):
    return 2/(3*sqrt(Omega-1)) * atan(sqrt(Omega-1))
# Omega < 1
def T2(Omega):
    return 2/(3*sqrt(1-Omega)) * atanh(sqrt(1-Omega))

def T(Omega):
    if Omega > 1:
        return T1(Omega)
    elif Omega == 1:
        return t0(Omega)
    else:
        return T2(Omega)

SymPy Plotting Backends で少し込み入った関数を plot() するときの小技

通常ならば

plot(T(Omega), (Omega, 0.001, 2))

とするところだが，関数 T(Omega) の定義が少し込み入っていて if 文がある場合，これだとエラーになる。

対策は，

  • 関数名のみにする（かっこや引数は省略する）
  • force_real_eval=True のオプションをつける

つまり，以下のようにすればよい。

plot(T, (Omega, 0.001, 2), force_real_eval=True)

オプションを設定して plot() する

In [5]:
var('Omega')
p = plot(
    (T, (Omega, 0.001, 2), 
     r"$\Omega_{\Lambda} = 1 - \Omega_{\rm m}$", 
     {'color':'red', 'linewidth':2}), 
    (t, (Omega, 0.001, 2), 
     r"$\Omega_{\Lambda} = 0$", 
     {'color':'black', 'linewidth':2}),
    force_real_eval=True, show=False
)

# 各種オプションの設定
ax = p.ax

#グラフの上下左右に目盛線を付ける
ax.yaxis.set_ticks_position('both')
ax.xaxis.set_ticks_position('both')
# tick の向きを内側に
ax.tick_params(which ="both",direction = 'in')

# 横軸の表示範囲
ax.set_xlim([0, 2])
ax.set_ylim([0.5, 1.5])

# x の主目盛
ax.set_xticks([0.2*i for i in range(1,11)])
# y の主目盛
ax.set_yticks([0.5+0.1*i for i in range(11)])

# 副目盛も表示
ax.minorticks_on()
# 副目盛には grid をつけない
ax.grid(False, which="minor")

# 主目盛の grid を細い点線で
ax.grid(True, which="major", linestyle='dotted');

# タイトル
ax.set_title(r"宇宙年齢 $t_0$ の密度パラメータ依存性")
# 軸のラベル
ax.set_xlabel(r"$\Omega_{\rm m}$")
ax.set_ylabel(r"$H_0\ t_0$")

# 凡例の表示サイズ
ax.legend(prop={"size":"x-large"});