next up previous
Next: About this document ... Up: 8 常微分方程式の数値解法:境界値問題 Previous: 3 現実的な恒星内部構造モデル

Subsections


4 Henyey法

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

問題とする微分方程式を

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

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

を解かねばならないことになり、これを直接、 $4000\times 4000$の行列 $\mbox{\boldmath$A$}$の逆行列を計算して
\begin{displaymath}
\mbox{\boldmath$x$}=\mbox{\boldmath$A$}^{−1}\mbox{\boldmath$b$}
\end{displaymath} (107)

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

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

\begin{displaymath}
y_i^n=y_i(x_n).
\end{displaymath} (108)

微分を差分に置き換えると、(8.105)式は
\begin{displaymath}
{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
\end{displaymath} (109)

となる。 但し、
\begin{displaymath}
\Delta x=x_{n+1}-x_n;\quad n=1, 2, \cdots, N-1
\end{displaymath} (110)

である。 中心と表面での境界条件を、それぞれ、
\begin{displaymath}
g_l(y_j^1)=0 ;\quad l=1, 2;\quad j=1,\cdots,4
\end{displaymath} (111)


\begin{displaymath}
g_l(y_j^N)=0 ;\quad l=3, 4;\quad j=1,\cdots,4
\end{displaymath} (112)

としよう。

1 線形化

このままでは、方程式系は非線形である。推定値からのずれを取り扱うことにして、 問題を線形化しよう。 即ち、真の解$y_i$についての推定値を$\bar y_i$として、 真の値からのずれを$\delta y_i$とする:
\begin{displaymath}
y_i^n=\bar y_i^n+\delta y_i^n;\quad i=1, \cdots, 4;
\quad n=1, 2, \cdots, N.
\end{displaymath} (113)

(8.113)式を (8.109), (8.111), 及び(8.112)式に代入して、 $\delta y_i^n$の2次以上の項を無視すると、 次の$4N$本の式を得る:
\begin{displaymath}
\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,
\end{displaymath} (114)


\begin{displaymath}
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,
\end{displaymath} (115)


\begin{displaymath}
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.
\end{displaymath} (116)

整理すると、
\begin{displaymath}
\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,
\end{displaymath} (117)


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

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

但し、ここに
\begin{displaymath}
A^n_{ij}={{\Delta x^n}\over{2}}F^n_{ij}+\delta_{ij},
\end{displaymath} (120)


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


\begin{displaymath}
\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]
\end{displaymath} (122)

であり、 $\delta_{ij}$は Kroneckerのデルタを表す。 また、
\begin{displaymath}
\bar f_i^n\equiv f_i(\bar y_j^n);
\quad i,j=1,\cdots,4;\quad n=1,2,\cdots,N,
\end{displaymath} (123)


\begin{displaymath}
\bar g_l^1\equiv g_l(\bar y_j^1);
\quad l=1,2;\quad j=1,\cdots,4,
\end{displaymath} (124)


\begin{displaymath}
\bar g_l^N\equiv g_l(\bar y_j^N);
\quad l=3,4;\quad j=1,\cdots,4,
\end{displaymath} (125)


\begin{displaymath}
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,
\end{displaymath} (126)


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


\begin{displaymath}
G^N_{lj}={\partial\bar g_l^N\over\partial\bar y_j^N};
\quad l=3,4;\quad i=1,\cdots,4
\end{displaymath} (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)を行列で表すと、
\begin{displaymath}
\left(\begin{array}{cccccc}
\mbox{\boldmath$G$}_1 & 0 & 0 & ...
...h$d$}}_{N-1}\\
\bar{\mbox{\boldmath$g$}}_N
\end{array}\right)
\end{displaymath} (129)

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

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

まず、2本の中心での 境界条件の式(8.118)と、最初の差分方程式(8.117)のうちから2本 ($i=1$, 2)をセットにして考える。この4本の式は次の様に表される:
\begin{displaymath}
(\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).
\end{displaymath} (130)

ここに
\begin{displaymath}
(\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},
\end{displaymath} (131)


\begin{displaymath}
(\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},
\end{displaymath} (132)


\begin{displaymath}
(\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}
\end{displaymath} (133)

である。 この4本の方程式は、8個の未知変数$\delta y^1_i$$\delta y^2_i$ $(i=1, \cdots, 4)$を含むから、 $4$個の変数は、残り4個の変数の線形従属として表すことが出来る。即ち、
\begin{displaymath}
\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),
\end{displaymath} (134)

但し、ここに
\begin{displaymath}
(\mbox{\boldmath$R$}_1)
=(\mbox{\boldmath$P$}_1)^{-1}(\mbox{\boldmath$Q$}_1),
\end{displaymath} (135)


\begin{displaymath}
(\mbox{\boldmath$r$}_1)
=(\mbox{\boldmath$P$}_1)^{-1}(\mbox{\boldmath$q$}_1)
\end{displaymath} (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)式の残りの2本 $(i=3,4)$ と、$n=2$の(8.117)式の最初の2本 $(i=1, 2)$ のセットを考える。 このセットは

\begin{displaymath}
(\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),
\end{displaymath} (137)

と表される。但し、ここで
\begin{displaymath}
(\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},
\end{displaymath} (138)


\begin{displaymath}
(\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},
\end{displaymath} (139)


\begin{displaymath}
(\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},
\end{displaymath} (140)

それに
\begin{displaymath}
(\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}
\end{displaymath} (141)

である。 (8.134)式を(8.137)式に代入して、
\begin{displaymath}
\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),
\end{displaymath} (142)

を得る。但し、ここで、行列 $(\mbox{\boldmath$R$}_2)$ $(\mbox{\boldmath$r$}_2)$
\begin{displaymath}
(\mbox{\boldmath$R$}_2)
=\left[(\mbox{\boldmath$S$}_1)(\mb...
...
+(\mbox{\boldmath$P$}_2)\right]^{-1}(\mbox{\boldmath$Q$}_2),
\end{displaymath} (143)


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

で定義する。

(8.142)式は

\begin{displaymath}
\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
\end{displaymath} (145)

と一般化した漸化式に書くことが出来る。但し、
\begin{displaymath}
(\mbox{\boldmath$R$}_n)
=\left[(\mbox{\boldmath$S$}_{n-1})...
...
+(\mbox{\boldmath$P$}_n)\right]^{-1}(\mbox{\boldmath$Q$}_n),
\end{displaymath} (146)


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

であり、
\begin{displaymath}
(\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},
\end{displaymath} (148)


\begin{displaymath}
(\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},
\end{displaymath} (149)


\begin{displaymath}
(\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},
\end{displaymath} (150)

それに
\begin{displaymath}
(\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}
\end{displaymath} (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)式の残りの2本の式と 表面での2本の境界条件 (8.112)式とで

\begin{displaymath}
(\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$}),
\end{displaymath} (152)

を得る。但し、ここで、
\begin{displaymath}
\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},
\end{displaymath} (153)


\begin{displaymath}
(\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},
\end{displaymath} (154)

それに
\begin{displaymath}
(\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}
\end{displaymath} (155)

である。 $n=N-1$の(8.145)式を(8.152)式に代入すれば、 $4$次の正方行列 $(\mbox{\boldmath$Z$})$
\begin{displaymath}
(\mbox{\boldmath$Z$})\equiv
\left[(\mbox{\boldmath$\Sigma$})(\mbox{\boldmath$R$}_{N-1})
+(\mbox{\boldmath$\Pi$})\right]
\end{displaymath} (156)

と 列ベクトル $(\mbox{\boldmath$c$})$
\begin{displaymath}
(\mbox{\boldmath$c$})\equiv
[(\mbox{\boldmath$w$})-(\mbox{\boldmath$\Sigma$})
(\mbox{\boldmath$r$}_{N-1})]
\end{displaymath} (157)

を使って、
\begin{displaymath}
(\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$})
\end{displaymath} (158)

となるが、 ($Z$)の逆行列を (8.158)式にあてがえば、
\begin{displaymath}
\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$}).
\end{displaymath} (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法の優越性 が顕著となる。


next up previous
Next: About this document ... Up: 8 常微分方程式の数値解法:境界値問題 Previous: 3 現実的な恒星内部構造モデル
Jun Makino
平成15年6月5日