二分法
非線形方程式 $f(x) = 0$ の解を二分法で求めてみよう。
-
$f(x_{min}) \cdot f(x_{max}) < 0$ となるような $x_{min},\ x_{max}$ の対を見つける。
-
$x_0 = (x_{min} + x_{max}) / 2$ をとり、$f(x_0)$ を計算する。
-
if $f(x_{min}) \cdot f(x_{0}) < 0$ → 解は $x_{min}$ と $x_0$ の間にある。$x_{max} = x_0$ とおき、ステップ2に戻る。
if $f(x_{min}) \cdot f(x_{0}) > 0$ → 解は $x_0$ と $x_{max}$ の間にある。$x_{min} = x_0$ とおき、ステップ2に戻る。
if $f(x_{min}) \cdot f(x_{0}) = 0$ → 解は $x_0$
こうして、必要な精度 $\Delta x = |x_{max} - x_{min}|$ までステップ2, 3 を繰り返す。この方法を用いて、次の2つの方程式の解をそれぞれ求めなさい。
$$\sin x - \cos x = 0\ \ \ \ \ (0 < x < \pi)$$
$$x - 2 \ln (x+1) = 0\ \ \ \ \ (x > 0)$$
宇宙の晴れ上がり
ビッグバン時において高温高圧の状態にあった宇宙の物質は、宇宙膨張にしたがって次第に冷えていく。やがて温度が数千度のなった時点で、電子と陽子は
$$\mathrm{e}^- + \mathrm{p}^+ \rightarrow \mathrm{H} + \gamma$$
という反応を起こして中性化する。水素の電離度を $x$ とすると、電子、陽子、中性水素の数は、それぞれ全粒子数 $n$ を用いて
$$n_{\mathrm e} = n_{\mathrm p} = nx,\ \ \ \ \ n_{\mathrm H} = (1-x)n$$
と表せる。ここで $x$ は、次の中性化の反応の式（サハの電離公式）を用いて求められる。
$$\frac{n_{\mathrm e} n_{\mathrm p}}{n_{\mathrm H}n} = \frac{x^2}{1-x} = \frac{(2\pi m_{\mathrm e}kT)^{3/2}}{n(2\pi h)^3}e^{-B/kT} \tag{1}$$
さらに、赤方偏移 $z$ における宇宙膨張の式や温度降下の式
$$n=1.12\times 10^{-5}\Omega_Bh^2(1+z)^3\ \ \mathrm{cm}^{-3},\ \ \ \ \ T = 2.736(1+z)\ \ \mathrm{K}$$
を用いると、式 (1) は数値的に次のように書ける。
$$\log \left[\frac{x^2}{1-x}\right] = 20.99 - \log\left[\Omega_B h^2(1+z)^{3/2} \right] - \frac{25050}{1+z} \tag{2}$$
式 (2) を解いて、電離度 $x$ を $z$ の関数として求め図示せよ。特に $x = 1/2$ のとき $z$ はいくらか、また温度は何度か。なおここで $\Omega _B h^2 = 0.02$ とする。
参考文献：
・P. J. E. Peebles, "Principles of Physical Cosmology", pp.165-167
・シリーズ現代の天文学2「宇宙論I」第3章
多項式によるフィッティング
与えられた $I$ 個のデータ点 $(x_1, y_1), \cdots, (x_I, y_I)$ を、最小自乗法を用いて直線 $f(x) = a + bx$ でフィットしてみよう。それには
$$g(a,b) = \sum_{i=1}^I(y_i - a - bx_i)^2$$
で定義された関数の値が、極小をとるような $a$ と $b$ を求めればよい。すなわち
\begin{array}{l}
0 &=& \partial g / \partial a &=& -2 \sum (y_i - a - bx_i),\\
0 &=& \partial g / \partial b &=& -2 \sum x_i(y_i - a - bx_i)
\end{array}
ここで $S_x = \sum x_i$, $\ S_y = \sum y_i$, $\ S_{xx} = \sum x_i^2$, $\ S_{xy} = \sum x_i y_i$ とおくと、上式は
\begin{array}{l}
a I + b S_x &=& S_y, \\
a S_x + b S_{xx} &=& S_{xy}
\end{array}
と書ける。これを解いて、$a$ と $b$ を求めればよい。
自分で適当なサンプルを作って、最小自乗法でフィットせよ。また N 次式
$$y = a_0 + a_1x + \cdots + a_Nx^N$$
でのフィッティングの公式を求めて、やはり適当なサンプルのフィッティングに用いてみよ。
球対称降着流（Bondi Flow）
球対称（外部）重力場の中の、定常でバロトロピック（$P = P(\rho)$）な流れを考える。基本方程式は、連続の式とベルヌーイの式である。
\displaylines{
4\pi r^2 \rho u = - \dot M = \mathrm{const.} \\
\frac{u^2}{2} + h(\rho) - \frac{GM}{r} = 0
}
ここで、$\dot M$ は質量流量（内向きを正とする）、$h(\rho)$ は等温気体の場合
$$h(\rho) \equiv \int \frac{dP}{\rho} = c_{\infty}^2 \ln\left(\frac{\rho}{\rho_{\infty}} \right)$$
となる。ここで $c_{\infty} = \sqrt{P_{\infty}/\rho_{\infty}}$, $\ P_{\infty}$, $\ \rho_{\infty}$ はそれぞれ無限遠での音速、圧力、密度である。次に、中心からの距離、速度（内向きを正）、密度、ガス降着率をそれぞれ
\begin{array}{l}
x &\equiv& r/r_B,\ \ \ \ \ r_B\ \ \ \equiv\ \ \ GM / c_{\infty}^2 \\
v &\equiv& |u| / c_{\infty} \\
\alpha &\equiv& \rho / \rho_{\infty} \\
\lambda &\equiv& \dot M / (4\pi r_B^2 \rho_{\infty} c_{\infty})
\end{array}
と規格化し、これらの規格化変数を使って上式を書き換えると
\displaylines{
x^2 \alpha (x)v(x) = \lambda \\
\frac{v(x)^2}{2} + \ln [\alpha(x)] - \frac{1}{x} = 0
}
となる。与えられた $\lambda$ について、上式を $\alpha (x)$ と $v(x)$ について解き、グラフに描け。特に $\lambda_e = e^{3/2}/4$ を境に、解のふるまいがどのように変わるか調べよ。
参考文献：
・Shapiro & Teukolsky, "Black Holes, White Dwarfs, and Newtron Stars", Sect. 14.3
・現代天文学講座6「恒星の世界」6.4章
銀河ポテンシャル中の星の運動
銀河ポテンシャル（isochrone potential）
$$\Phi(r) = - \frac{GM}{b + \sqrt{b^2 + r^2}}\ \ \ \ \ \ \ \ (b = \mathrm{const.})$$
の中にある質量 $m$ の星の運動を考える。任意の初期条件の下で運動方程式を解き、その軌道を図示せよ。円柱座標 $(r, \phi, z)$ を用いて $z=0$ 面での運動を考え、力のベクトル $\bf{F}$ がポテンシャル $\Phi(r)$ を使って $\bf{F} = -\nabla \Phi$ で表されることを用い、動径 $r$ と方位角 $\phi$ それぞれの方向について運動方程式をたてて解けばよい。また星の全エネルギーが $E(<0)$ で与えられるとき、星の軌道の動径方向の振動の周期は、角運動量によらず一定値
$$T_r = \frac{2\pi G M}{(-2E)^{3/2}}$$
をとることを、数値的に確認せよ。
参考文献：
・Binney & Tremaine, "Galactic Dynamics", pp.103-110
星の内部構造（Lane-Emden 方程式）
星の内部構造は、力学平衡（圧力と重力のつり合い）の式と連続の式（質量保存の式）を連立して解くことで求められる。
\frac{dP}{dr} = -\frac{GM(r)}{r^2}\rho,\ \ \ \ \ \frac{dM(r)}{dr} = 4\pi r^2 \rho \\
\frac{1}{r^2}\frac{d}{dr}\left(\frac{r^2}{\rho}\frac{dP}{dr} \right) = -4\pi G\rho
さらにポリトロピックの仮定（$P \propto \rho ^r$）をして、次の変数変換を行う
\rho = \lambda \theta ^n,\ \ \ \ \ P = K\rho ^{1 + (1/n)} = K\lambda^{1 + (1/n)}\theta^{n+1},\ \ \ \ \ r = \alpha \xi = \left[\frac{(n+1)K}{4\pi G}\lambda^{(1/n)-1} \right]^{1/2} \xi\ \ \ \ \ \ \ \ (\lambda = \mathrm{const.})
これらの変数を使うことにより、Lane-Emden 方程式
$$\frac{1}{\xi ^2}\frac{d}{d\xi}\left(\xi ^2 \frac{d\theta}{d\xi} \right) = -\theta ^n$$
が得られる。この2階常微分方程式を、$n=0, 1, 2, 3, 4, 5$ の各場合について、それぞれ初期値
$$\theta = 1,\ \ \ \ \ \ \frac{d\theta}{d\xi} = 0,\ \ \ \ \ \xi = 0$$
のもとで解き、星の内部における密度、温度分布等を図示せよ。なお $n = 0, 1, 5$ の場合については解析解が見つかっているので、教科書等で調べて数値的に求めた解をチェックせよ。
参考文献：
・Chandrasekhar, "Stellar Structure and Stellar Atmosphere", pp84-95
・チャンドラセカール著「星の構造」
重力レンズ
一般相対性理論によると、質量 $M$ の天体による強い重力場中を伝搬する光の経路は、この天体を原点とした極座標 ($r, \phi$) を用いて
$$\frac{d^2}{d\phi ^2}\left(\frac{1}{r} \right) + \frac{1}{r} = \frac{3GM}{c^2}\frac{1}{r^2}$$
という式で表される。任意の初期値 $(r_0, \phi_0)$ のもとでこの式を積分することにより、光の経路を調べよ。ここで
$$u = \frac{1}{r},\ \ \ \ \ v = \frac{du}{d\phi}$$
と変数変換をすると解きやすくなる。相対性理論を考慮すると、$u$ と $v$ の初期値はそれぞれ
$$u_0 = \frac{1}{r_0},\ \ \ \ \ v_0 = \frac{(1-r_g/r_0)^{1/2}}{r_0 \tan \phi_0}\ \ \ \ \ \left( r_g \equiv \frac{2GM}{c^2}\right)$$
となる。また衝突パラメータ $b$ が、臨界値
$$b_c = (27)^{1/2}\frac{GM}{c^2}$$
を下回ると、光は天体から外に出てこられないことを示せ。
参考文献：
・Shapiro & Teukolsky, "Black Holes, White Dwarfs, and Newtron Stars", Sect. 12.5
・福江純著「ブラックホールの世界」7章
降着円盤
活動銀河核（Active Galactic Nuclei）や近接連星系の中心には、ブラックホールや中性子星、白色矮星などの高密度天体（コンパクト星）があると考えられている。その中心星の重力に引き寄せられたガスは、中心星の周りを速度
$$v_{\phi} = \sqrt{\frac{GM}{r}}$$
でケプラー回転しながら、ゆっくりと中心星へと落ちていく。こうして星の周りに降着円盤と呼ばれるガス円盤ができあがる。また星形成の際に中心星の周りに形成される原始惑星系円盤（惑星形成の場）も、同様の過程を経て中心星へと降着する。この円盤の時間進化を解いてみよう。基本方程式は、円盤の面密度を
$$\Sigma (r) \equiv \int \rho(r,z)dz$$
として、
\frac{\partial \Sigma}{\partial t} = \frac{1}{2\pi r}\frac{\partial \dot M}{\partial r}\\
\dot M = 6\pi \sqrt r \frac{\partial}{\partial r} \left(\sqrt r \nu \Sigma \right)
となる。ここで、$\dot M \equiv -2\pi r \Sigma v_r$ は質量降着率（$v_r$ はガスの動径方向の速度成分）、$\nu$ は粘性である。上の2式を組み合わせると、拡散方程式
$$\frac{\partial \Sigma}{\partial t} = \frac{3}{r}\frac{\partial}{\partial r}\sqrt r \frac{\partial}{\partial r}\left(\sqrt r \nu \Sigma \right)$$
になる。さらに物理量を、次のような無次元の変数に変換しよう。
$$x \equiv \sqrt{\frac{r}{r_{\mathrm out}}},\ \ \ \ \ \sigma \equiv \frac{\Sigma}{\Sigma_0},\ \ \ \ \ \tau \equiv \frac{t}{t_0},\ \ \ \ \ t_0 \equiv \frac{4r_{\mathrm out}^2}{3 \nu}$$
ここで $r_{\mathrm out}$ は円盤の外縁の半径（$0 < x \leq 1$）、$\Sigma _0$ は任意の定数、また簡単のため $\nu$ は定数とした。上の無次元変数を用いると、最終的に基本方程式 (3) は
$$\frac{\partial \sigma}{\partial \tau} = \frac{1}{x^3}\frac{\partial ^2}{\partial x^2}(x\sigma)$$
となる。次の初期条件（境界条件）のもとで上式を解き、時間がたつにつれ、円盤中のガスは次第に中心星へと落ちていくことを示せ。
\sigma(x,\tau) = \left\{
\begin{array}{l}
1.0\ \ \ \ \ \ \ \ at\ \ |x-0.5| \leq 0.1 \\
0.0\ \ \ \ \ \ \ \ at\ \ |x-0.5| > 0.1
\end{array}
\right.
また、ガスの持つ角運動量（$x\sigma(x)$）が内向きに輸送されるか、それとも外向きに輸送されるかを調べよ。
参考文献：
・Kato, Fukue & Mineshige, "Black-Hole Accretion Disks", Sect. 3.1
・柴田一成他著「活動する宇宙」7章