Subsections

• 1 線形化
• 2 係数行列は帯行列
• 3 帯行列を区分化して解く：漸化式
• 4 帯行列を区分化して解く：後退代入
• 5 Henyey法の利点

4 Henyey法

星の内部構造を解く問題は、非線形の境界値問題であることを既に述べた。 主系列近傍の恒星の内部構造は、フィッティング法で解けた。 1960年代に恒星進化の理論は、主系列から赤色巨星への進化を追うようになったが、 赤色巨星になると、表面での境界値に対する依存性が鈍くなって、 フィッティング法では上手く解くことが出来なくなっていた。 こういった状況をブレークしたのは、カリフォルニア大学の Henyeyであった。彼は、一種のrelaxation methodで代数方程式を解くことに 成功したのである 。 この方法は、一般の２点境界値問題に適用出来る方法で、 天文学の社会では、Henyey法と呼ばれている。 以下に、このHenyey法を解説しよう 。

問題とする微分方程式を

$${{dy_i}\over{dx}}=f_i(x, y_j) \qquad i,j=1,\cdots,4$$ (105)

とする。 独立変数$x$の領域を$N$メッシュに切って、微分方程式を差分方程式に書き換え、 これを数値的に解くことを考える。元の微分方程式を４階の微分方程式だとすると、 計$4(N-1)$本の差分方程式となる。これらは４つの変数$y_i$ $(i=1, \cdots, 4)$の $N$個のメッシュそれぞれでの値、合計 $4N$個の未知数、$y_i(x_n)$、の代数方程式と見なすことができる。 未知数の数と方程式の数が合わない？ このままでは、その通り。しかし、境界条件があることをお忘れなく。 内側の境界に境界条件が２本、外側の境界にも２本、合計４本の境界条件 の式があるので、方程式の数は合計$4N$本になって、未知数の数と同じになり、 代数方程式は原理的に解けるのである。 ただ、これは「原理的に」である。 この方程式が非線形であれば、そのまま解くのが容易でないのは明らかであり、 線形であったとしても、 実際に数値的に解けるかどうかは、 $N$による。実際的には、$N=1000$とすると、$4000$元の代数方程式

$$\mbox{\boldmath$A$}\mbox{\boldmath$x$}=\mbox{\boldmath$b$}$$ (106)

を解かねばならないことになり、これを直接、 $4000\times 4000$の行列 $\mbox{\boldmath$A$}$の逆行列を計算して

$$\mbox{\boldmath$x$}=\mbox{\boldmath$A$}^{−1}\mbox{\boldmath$b$}$$ (107)

とするのは、実際的ではない。 $4000\times 4000$もの大行列の逆行列 を求めるのが困難だからである。

従属変数$y_i$の$x=x_n$での値であることを表わすのに、上付き添え字$n$をつけて $y_i^n$で表わすことにする：

$$y_i^n=y_i(x_n).$$ (108)

微分を差分に置き換えると、(8.105)式は

$${y^{n+1}_i-y^n_i\over\Delta x} ={1\over 2}[f_i(y_j^{n+1})+f_i(y_j^n)]; \quad i,j=1,\cdots,4;\quad n=1,2,\cdots,N-1$$ (109)

となる。 但し、

$$\Delta x=x_{n+1}-x_n;\quad n=1, 2, \cdots, N-1$$ (110)

である。 中心と表面での境界条件を、それぞれ、

$$g_l(y_j^1)=0 ;\quad l=1, 2;\quad j=1,\cdots,4$$ (111)

と

$$g_l(y_j^N)=0 ;\quad l=3, 4;\quad j=1,\cdots,4$$ (112)

としよう。

1 線形化

このままでは、方程式系は非線形である。推定値からのずれを取り扱うことにして、 問題を線形化しよう。 即ち、真の解$y_i$についての推定値を$\bar y_i$として、 真の値からのずれを$\delta y_i$とする：

$$y_i^n=\bar y_i^n+\delta y_i^n;\quad i=1, \cdots, 4; \quad n=1, 2, \cdots, N.$$ (113)

(8.113)式を (8.109), (8.111), 及び(8.112)式に代入して、 $\delta y_i^n$の２次以上の項を無視すると、 次の$4N$本の式を得る：

$$\bar y_i^{n+1} +\delta y_i^{n+1} - \bar y_i^n - \delta y_i^n... ...n})}\over{\partial \bar y_j^{n}}} \right]; \quad i=1,\cdots,4,$$ (114)

$$g_l(\bar y_j^1) +\sum_{j=1}^4{{\partial g_l(\bar y_j^1)}\over{\partial\bar y_j^1}} =0; \quad l=1, 2,$$ (115)

$$g_l(\bar y_j^N) +\sum_{j=1}^4{{\partial g_l(\bar y_j^N)}\over{\partial\bar y_j^N}} =0; \quad l=3, 4.$$ (116)

整理すると、

$$\sum_{j=1}^4 A^n_{ij}\delta y^{n}_j + \sum_{j=1}^4 B^n_{ij... ...y^{n+1}_j +d^n_i=0;\quad i=1,\cdots,4;\quad n=1,2,\cdots,N-1,$$ (117)

$$\sum_{j=1}^4 G^1_{lj}\delta y^1_j+\bar g^1_l=0; \quad l=1,2,$$ (118)

それに

$$\sum_{j=1}^4 G^N_{lj}\delta y^N_j+\bar g^N_l=0; \quad l=3,4,$$ (119)

但し、ここに

$$A^n_{ij}={{\Delta x^n}\over{2}}F^n_{ij}+\delta_{ij},$$ (120)

$$B^n_{ij}={{\Delta x^n}\over{2}}F^{n+1}_{ij}-\delta_{ij},$$ (121)

$$\bar d^n_i =-(\bar y_i^{n+1}-\bar y_i^n)+{{\Delta x}\over{2}}[\bar f^{n+1}_i +\bar f^n_i]$$ (122)

であり、 $\delta_{ij}$は Kroneckerのデルタを表す。 また、

$$\bar f_i^n\equiv f_i(\bar y_j^n); \quad i,j=1,\cdots,4;\quad n=1,2,\cdots,N,$$ (123)

$$\bar g_l^1\equiv g_l(\bar y_j^1); \quad l=1,2;\quad j=1,\cdots,4,$$ (124)

$$\bar g_l^N\equiv g_l(\bar y_j^N); \quad l=3,4;\quad j=1,\cdots,4,$$ (125)

$$F^n_{ij}={\partial\bar f_i^n\over\partial\bar y_j^n}; \quad i,j=1,\cdots,4;\quad n=1,2,\cdots,N,$$ (126)

$$G^1_{lj}={\partial\bar g_l^1\over\partial\bar y_j^1}; \quad l=1,2;\quad i=1,\cdots,4,$$ (127)

$$G^N_{lj}={\partial\bar g_l^N\over\partial\bar y_j^N}; \quad l=3,4;\quad i=1,\cdots,4$$ (128)

で定義する。

2 係数行列は帯行列

方程式(8.117)-(8.119)は、$4N$個の未知変数$\delta y_i^n$に関する 線形方程式である。 一見したところ、この方程式系は大行列を要し、解くことが出来ないと思えるかも知 れない。しかし、係数行列の要素の多くはゼロである。実際、 $\mbox{\boldmath$A$}_n \equiv \{A_{ij}^n\}$、 $\mbox{\boldmath$B$}_n \equiv \{B_{ij}^n\}$、 $\mbox{\boldmath$G$}_1 \equiv \{G_{lj}^1\}$、 $\mbox{\boldmath$G$}_N \equiv \{G_{lj}^N\}$、 $\mbox{\boldmath$\delta$}\mbox{\boldmath$y$}_n \equiv \{\delta y_i^n\}$等と記して、方程式(8.117)-(8.119)を行列で表すと、

$$\left(\begin{array}{cccccc} \mbox{\boldmath$G$}_1 & 0 & 0 & ... ...h$d$}}_{N-1}\\ \bar{\mbox{\boldmath$g$}}_N \end{array}\right)$$ (129)

になり、左辺の$4N\times 4N$の係数行列の対角帯以外はゼロである。 Henyey法は、この特質を有効に利用するのである。

3 帯行列を区分化して解く：漸化式

まず、２本の中心での 境界条件の式(8.118)と、最初の差分方程式(8.117)のうちから２本 ($i=1$, 2)をセットにして考える。この４本の式は次の様に表される：

$$(\mbox{\boldmath$P$}_1)\pmatrix{\delta y^1_1\cr \delta y^1_... ...\cr \delta y^2_3\cr\delta y^2_4\cr}+(\mbox{\boldmath$q$}_1).$$ (130)

ここに

$$(\mbox{\boldmath$P$}_1)=\pmatrix{G^1_{11}&G^1_{12}&G^1_{13}... ...A^1_{13}&A^1_{14}\cr A^1_{21}&A^1_{22}&A^1_{23}&A^1_{24}\cr},$$ (131)

$$(\mbox{\boldmath$Q$}_1)=\pmatrix{0& 0& 0& 0\cr 0& 0& 0& 0\... ...13}&-B^1_{14}\cr -B^2_{21}&-B^2_{22}&-B^2_{23}&-B^2_{24}\cr},$$ (132)

$$(\mbox{\boldmath$q$}_1)= \pmatrix{-\bar g_1^1\cr -\bar g_2^1\cr -\bar d_1^1\cr -\bar d_2^1\cr}$$ (133)

である。 この4本の方程式は、8個の未知変数$\delta y^1_i$、 $\delta y^2_i$ $(i=1, \cdots, 4)$を含むから、 $4$個の変数は、残り4個の変数の線形従属として表すことが出来る。即ち、

$$\pmatrix{\delta y^1_1\cr \delta y^1_2\cr \delta y^1_3\cr\... ...\cr \delta y^2_3\cr\delta y^2_4\cr}+(\mbox{\boldmath$r$}_1),$$ (134)

但し、ここに

$$(\mbox{\boldmath$R$}_1) =(\mbox{\boldmath$P$}_1)^{-1}(\mbox{\boldmath$Q$}_1),$$ (135)

$$(\mbox{\boldmath$r$}_1) =(\mbox{\boldmath$P$}_1)^{-1}(\mbox{\boldmath$q$}_1)$$ (136)

であり、 $(\mbox{\boldmath$P$}_1)^{-1}$ は $(\mbox{\boldmath$P$}_1)$の逆行列を表す。 (8.134)式は、一旦 $\mbox{\boldmath$\delta$}\mbox{\boldmath$y$}_2$が求まれ ば、 $\mbox{\boldmath$\delta$}\mbox{\boldmath$y$}_1$が求められる ことを意味している。

次に、$n=1$の(8.117)式の残りの２本 $(i=3,4)$ と、$n=2$の(8.117)式の最初の２本 $(i=1, 2)$ のセットを考える。 このセットは

$$(\mbox{\boldmath$S$}_1)\pmatrix{\delta y^1_1\cr \delta y^1_... ...\cr \delta y^3_3\cr\delta y^3_4\cr}+(\mbox{\boldmath$q$}_2),$$ (137)

と表される。但し、ここで

$$(\mbox{\boldmath$S$}_1)=\pmatrix{A^1_{31}&A^1_{32}&A^1_{33}... ...&A^2_{42}&A^2_{43}&A^2_{44}\cr 0& 0& 0& 0\cr 0& 0& 0& 0\cr},$$ (138)

$$(\mbox{\boldmath$P$}_2)=\pmatrix{B^1_{31}&B^1_{32}&B^1_{33}... ...A^2_{13}&A^2_{14}\cr A^2_{21}&A^2_{22}&A^2_{23}&A^2_{24}\cr},$$ (139)

$$(\mbox{\boldmath$Q$}_2)=\pmatrix{0& 0& 0& 0\cr 0& 0& 0& 0\... ...13}&-B^2_{14}\cr -B^2_{21}&-B^2_{22}&-B^2_{23}&-B^2_{24}\cr},$$ (140)

それに

$$(\mbox{\boldmath$q$}_2) =\pmatrix{{-\bar d^1_3}\cr {-\bar d^1_4}\cr {-\bar d^2_1}\cr {-\bar d^2_2}\cr}$$ (141)

である。 (8.134)式を(8.137)式に代入して、

$$\pmatrix{\delta y^2_1\cr \delta y^2_2\cr \delta y^2_3\cr\... ...\cr \delta y^3_3\cr\delta y^3_4\cr}+(\mbox{\boldmath$r$}_2),$$ (142)

を得る。但し、ここで、行列 $(\mbox{\boldmath$R$}_2)$ と $(\mbox{\boldmath$r$}_2)$ は

$$(\mbox{\boldmath$R$}_2) =\left[(\mbox{\boldmath$S$}_1)(\mb... ... +(\mbox{\boldmath$P$}_2)\right]^{-1}(\mbox{\boldmath$Q$}_2),$$ (143)

$$(\mbox{\boldmath$r$}_2) =\left[(\mbox{\boldmath$S$}_1)(\mb... ...q$}_2)-(\mbox{\boldmath$S$}_1) (\mbox{\boldmath$r$}_1)\right]$$ (144)

で定義する。

(8.142)式は

$$\pmatrix{\delta y^n_1\cr \delta y^n_2\cr \delta y^n_3\cr\... ...y^{n+1}_4\cr}+(\mbox{\boldmath$r$}_n); \quad n=1,2,\cdots,N-1$$ (145)

と一般化した漸化式に書くことが出来る。但し、

$$(\mbox{\boldmath$R$}_n) =\left[(\mbox{\boldmath$S$}_{n-1})... ... +(\mbox{\boldmath$P$}_n)\right]^{-1}(\mbox{\boldmath$Q$}_n),$$ (146)

$$(\mbox{\boldmath$r$}_n) =\left[(\mbox{\boldmath$S$}_{n-1})... ...(\mbox{\boldmath$S$}_{n-1})(\mbox{\boldmath$r$}_{n-1})\right],$$ (147)

であり、

$$(\mbox{\boldmath$S$}_{n-1})=\pmatrix{A^{n-1}_{31}&A^{n-1}_{... ...}&A^{n-1}_{43}&A^{n-1}_{44}\cr 0& 0& 0& 0\cr 0& 0& 0& 0\cr},$$ (148)

$$(\mbox{\boldmath$P$}_n)=\pmatrix{B^{n-1}_{31}&B^{n-1}_{32}&... ...A^n_{13}&A^n_{14}\cr A^n_{21}&A^n_{22}&A^n_{23}&A^n_{24}\cr},$$ (149)

$$(\mbox{\boldmath$Q$}_n)=\pmatrix{0& 0& 0& 0\cr 0& 0& 0& 0\... ..._{14}\cr -B^{n}_{21}&-B^{n}_{22}&-B^{n}_{23}&-B^{n}_{24}\cr},$$ (150)

それに

$$(\mbox{\boldmath$q$}_n) =\pmatrix{{-\bar d^{n-1}_3}\cr {-\bar d^{n-1}_4}\cr {-\bar d^n_1}\cr {-\bar d^n_2}\cr}$$ (151)

である。

4 帯行列を区分化して解く：後退代入

(8.145)式は、一旦 $\mbox{\boldmath$\delta$}\mbox{\boldmath$y$}_{n+1}$ が求まれば、 $\mbox{\boldmath$\delta$}\mbox{\boldmath$y$}_n$が求められる ことを意味している。問題をいつまでも先送りしているだけの様に見えるかもしれな いが、そうではない。

一番外側のメッシュ点$n=N$では、 $n=N-1$での(8.117)式の残りの２本の式と 表面での２本の境界条件 (8.112)式とで

$$(\mbox{\boldmath$\Sigma$})\pmatrix{\delta y^{N-1}_1\cr \del... ...\cr \delta y^N_3\cr\delta y^N_4\cr} =(\mbox{\boldmath$w$}),$$ (152)

を得る。但し、ここで、

$$\mbox{\boldmath$\Sigma$}=\pmatrix{A^{N-1}_{31}&A^{N-1}_{32}... ...}&A^{N-1}_{43}&A^{N-1}_{44}\cr 0& 0& 0& 0\cr 0& 0& 0& 0\cr},$$ (153)

$$(\mbox{\boldmath$\Pi$})=\pmatrix{B^{N-1}_{31}&B^{N-1}_{32}&... ...G^N_{33}&G^N_{34}\cr G^N_{41}&G^N_{42}&G^N_{43}&G^N_{44}\cr},$$ (154)

それに

$$(\mbox{\boldmath$w$}) =\pmatrix{-\bar d^{N-1}_{3}\cr {-\bar d^{N-1}_{4}}\cr {-\bar g^N_{3}}\cr {-\bar g^N_{4}}\cr}$$ (155)

である。 $n=N-1$の(8.145)式を(8.152)式に代入すれば、 $4$次の正方行列 $(\mbox{\boldmath$Z$})$

$$(\mbox{\boldmath$Z$})\equiv \left[(\mbox{\boldmath$\Sigma$})(\mbox{\boldmath$R$}_{N-1}) +(\mbox{\boldmath$\Pi$})\right]$$ (156)

と 列ベクトル $(\mbox{\boldmath$c$})$

$$(\mbox{\boldmath$c$})\equiv [(\mbox{\boldmath$w$})-(\mbox{\boldmath$\Sigma$}) (\mbox{\boldmath$r$}_{N-1})]$$ (157)

を使って、

$$(\mbox{\boldmath$Z$}) \pmatrix{\delta y^N_1\cr \delta y^N_2\cr\delta y^N_3\cr\delta y^N_I\cr} =(\mbox{\boldmath$c$})$$ (158)

となるが、 ($Z$)の逆行列を (8.158)式にあてがえば、

$$\pmatrix{\delta y^N_1\cr \delta y^N_2\cr\delta y^N_3\cr\delta y^N_I\cr} =(\mbox{\boldmath$Z$})^{-1}(\mbox{\boldmath$c$}).$$ (159)

を得る。こうして、ここにきて、一番外側のメッシュ点での解 $\delta y^N_i$ $(i=1, \cdots, 4)$が求められるのである。一旦、$\delta y^N_i$が決まれば、 これを$n=N-1$に対する(8.145)式に代入すれば、 $n=N-1$についての解 $\delta y^n_i$ を得る。同様の操作を繰り返し、 解 $\delta y^n_i$ が順次得られる寸法だ。 $\delta y^n_i$が充分小さくなるまで補正を繰り返してや れば、逐次解を得ることが出来る。

5 Henyey法の利点

こうして、非線形常微分方程式の境界値問題を、推定した解の近傍で線形化して、解 くことが出来ることが判ったと思う。線形化した点が、Henyey法の真髄ではない。 線形化して得た、方程式(8.129)を、左辺の大きい係数行列の逆行列を求め るのが困難なところを、小行列に区分化して、解いていった点が本質なのである。 元の係数行列は$4N\times 4N$次であった。それに対して、Henyey法では$4\times 4$の区分行列を$N$回取り扱うのである。行列の演算は、行列の次元の二乗に比例する から、元の係数行列そのものの逆行列演算は$16N^2$に比例するのに対し、Henyey法で は、$16N$にしか比例しない。従って、$N$が大きくなるにつれて、Henyey法の優越性 が顕著となる。