Return to 角径距離

補足：SymPy で角径距離のグラフを描く

角径距離の導出の詳細については，以下のページを参照。

ここでは，SymPy Plotting Backends を使って角径距離のグラフを描く。

モジュールの import と設定

In [1]:
from sympy.abc import *
from sympy import *
from spb import *

# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']

import matplotlib.pyplot as plt
plt.rcParams['mathtext.fontset'] = 'cm'

$\Omega_{\Lambda} = 0$ の場合の角径距離

\begin{eqnarray}
d_A
&=& \frac{2}{H_0 \Omega_{\rm m}^2 (1+z)^2} \left\{2 – \Omega_{\rm m} + \Omega_{\rm m} z – (2-\Omega_{\rm m}) \sqrt{1 + \Omega_{\rm m} z}\right\}
\end{eqnarray}

以下では $H_0$ を（$H_0 = 1$ として）省略し，表記の都合上 $\Omega_{\rm m} \rightarrow \Omega$

In [2]:
def dA(Omega, z):
    return 2/(Omega**2*(1+z)**2) * (2-Omega+Omega*z-(2-Omega)*sqrt(1 + Omega*z))

$\Omega_{\rm m} + \Omega_{\Lambda} = 1$ の場合の角径距離

\begin{eqnarray}
d_A
&=& \frac{1}{H_0 (1+z)} \int_0^z \frac{dz}{\sqrt{(1-\Omega_{\rm m}) + \Omega_{\rm m} (1+z)^3} }
\end{eqnarray}

解析的には積分できないので，.evalf() で数値積分に。（激遅だけど楽ちん。）

In [3]:
def dAL(Omega, z):
    return 1/(1+z)*integrate(1/sqrt((1-Omega) + Omega*(1+x)**3), (x, 0, z)).evalf()

$0 \leq z \leq 5$ までのグラフ

初めて dAL(0.3, z) を plot() する際にはいろいろ文句を言われるかもしれないが，なんとかグラフは描いてくれる。ものすごく時間がかかるので，辛抱。

In [4]:
%%time

plot((dAL(0.3, z), '$\Omega_{m}=0.3,\, \Omega_{\Lambda}=0.7$'), 
     (dA(0.3, z), '$\Omega_{m}=0.3,\, \Omega_{\Lambda}=0$'), 
     (dA(1, z), '$\Omega_{m}=1,\, \ \ \,\Omega_{\Lambda}=0$'), 
     (z, 0, 5), 
     xlim = (0, 5), ylim = (0, 0.5), title='角径距離',
     ylabel = '$H_0\, d_A(z)$');

/usr/local/lib/python3.8/dist-packages/spb/series.py:1087: UserWarning: NumPy is unable to evaluate with complex numbers some of the functions included in this symbolic expression: [hyper]. Hence, the evaluation will use real numbers. If you believe the resulting plot is incorrect, change the evaluation module by setting the `modules` keyword argument.
/usr/local/lib/python3.8/dist-packages/spb/series.py:262: UserWarning: The evaluation with NumPy/SciPy failed.
NameError: name 'hyper' is not defined
Trying to evaluate the expression with Sympy, but it might be a slow operation.


CPU times: user 12.9 s, sys: 1.41 s, total: 14.3 s
Wall time: 11.8 s

$ 0 \leq z \leq 20$ までのグラフ

$z$ の範囲を変え，若干の細かな設定を追加してみた。あいかわらず何か警告が出るが，気にせずに。

In [5]:
%%time

p = plot((dAL(0.3, z), '$\Omega_{m}=0.3,\, \Omega_{\Lambda}=0.7$'), 
         (dA(0.3, z), '$\Omega_{m}=0.3,\, \Omega_{\Lambda}=0$'), 
         (dA(1, z), '$\Omega_{m}=1,\, \ \ \,\Omega_{\Lambda}=0$'), 
         (z, 0, 20), 
         xlim = (0, 20), ylim = (0, 0.5), title='角径距離',
         ylabel = '$H_0\ d_A(z)$', show = False)

fig, ax = p.fig, p.ax

ax.grid(linestyle='dotted')
# x の主目盛を 5 刻みに
ax.set_xticks([i for i in range(0,21,5)])
# 副目盛も表示
ax.minorticks_on()
# 副目盛には grid をつけない
ax.grid(False, which="minor")
# x の副目盛を 1 刻みに
ax.set_xticks([i for i in range(1,20)], minor = True)
# y の副目盛を 1 刻みに
ax.set_yticks([0.05*i for i in range(1,10)], minor = True);

/usr/local/lib/python3.8/dist-packages/spb/series.py:1087: UserWarning: NumPy is unable to evaluate with complex numbers some of the functions included in this symbolic expression: [hyper]. Hence, the evaluation will use real numbers. If you believe the resulting plot is incorrect, change the evaluation module by setting the `modules` keyword argument.

CPU times: user 10.9 s, sys: 91.7 ms, total: 11 s
Wall time: 10.8 s

数値リストにして plot_list() を使ってグラフを描く

（NumPy を使わずに）dAL(Omega, z) の数値リストを作成し，plot_list() で数値データをグラフにしてみる。.evalf() は変わらないので list_dAL を求めるのは激遅のままだが，グラフを描く際に，なんだかんだと文句はつけられなくなるので，少しだけ精神的に楽かなと。（それにしても SymPy の数値積分は遅い。）

In [6]:
%%time

list_z = [0.1*i for i in range(201)]
list_dAL = [dAL(0.3, 0.1*i) for i in range(201)]

CPU times: user 2min 25s, sys: 206 ms, total: 2min 25s
Wall time: 2min 25s

参考：数値積分としてscipy.integrate.quad() を使う

.evalf() を使った SymPy の数値積分があまりにも遅すぎて精神衛生上よろしくないので，SymPy（のみ）を使うという本稿の方針からは逸脱するが，SciPy の数値積分 quad() を使ってみる。

quad() を使う際には以下のように，決め事にそって被積分関数を定義してやる。

  • 被積分関数は原則 $x$ の1変数関数の決め打ちで。
  • $x$ 以外のパラメータにも依存する場合は，$x$ の後に引数として並べる。
  • quad() を呼び出す際には，$x$ 以外のパラメータは args に入れる。
  • quad() は (積分値, 推定誤差) を返すので，積分値のみを知りたいときはそのように。
In [7]:
from scipy.integrate import quad
In [8]:
def f(x, Omega):
    return 1/sqrt((1-Omega) + Omega*(1+x)**3)

def dAL2(Omega, z):
    tmp = quad(f, 0, z, args=(Omega))
    return 1/(1+z)*tmp[0]
    # 別解
    # ans, err = quad(f, 0, z, args=(Omega))
    # return 1/(1+z)* ans
In [9]:
%%time

list_z = [0.1*i for i in range(201)]
list_dAL2 = [dAL2(0.3, 0.1*i) for i in range(201)]

CPU times: user 2.09 s, sys: 7.97 ms, total: 2.1 s
Wall time: 2.1 s

… ということで，SciPy の数値積分 quad() で計算した $d_A(\Omega, z)$ は約70倍高速であることがわかる。

$0 \leq z \leq 5$ までのグラフ

In [10]:
p1 = plot_list(
        list_z, list_dAL2, 
        '$\Omega_{m}=0.3,\, \Omega_{\Lambda}=0.7$', 
        xlim = (0, 5), ylim = (0, 0.5),
        title='角径距離',
        xlabel = '$z$', ylabel = '$H_0\, d_A(z)$', 
        show = False)
p2 = plot(
        (dA(0.3, z), '$\Omega_{m}=0.3,\, \Omega_{\Lambda}=0$'), 
        (dA(1, z), '$\Omega_{m}=1,\, \ \ \,\Omega_{\Lambda}=0$'), 
        (z, 0, 5), show = False)
p = p1 + p2

p.show();

$0 \leq z \leq 20$ までのグラフ

In [11]:
p1 = plot_list(
        list_z, list_dAL2, 
        '$\Omega_{m}=0.3,\, \Omega_{\Lambda}=0.7$', 
        xlim = (0, 20), ylim = (0, 0.5),
        title='角径距離',
        xlabel = '$z$', ylabel = '$H_0\, d_A(z)$', 
        show = False)
p2 = plot(
        (dA(0.3, z), '$\Omega_{m}=0.3,\, \Omega_{\Lambda}=0$'), 
        (dA(1, z), '$\Omega_{m}=1,\, \ \ \,\Omega_{\Lambda}=0$'), 
        (z, 0, 20), show = False)
p = p1 + p2

fig, ax = p.fig, p.ax

ax.grid(linestyle='dotted')
# x の主目盛を 5 刻みに
ax.set_xticks([i for i in range(0,21,5)])
# 副目盛も表示
ax.minorticks_on()
# 副目盛には grid をつけない
ax.grid(False, which="minor")
# x の副目盛を 1 刻みに
ax.set_xticks([i for i in range(1,20)], minor = True)
# y の副目盛を 1 刻みに
ax.set_yticks([0.05*i for i in range(1,10)], minor = True);