next up previous
Next: 5 構造体とクラス Up: 計算天文学 II 第5回 常微分方程式の初期値問題(1) Previous: 3 ルンゲ・クッタ以外の方法

Subsections


4 線形多段階法

RK 法は、「一段階」法である。これはどういう意味かというと、微分方程式

\begin{displaymath}
\frac{dy}{dx} = f(x,y)
\end{displaymath} (6)

$x=x_i$ での近似解 $y_i$ があったとして、次の $x$の値 $x_{i+1} =
x_i + \Delta x_i$ での近似解 $y_{i+1}$ を求めるのに、$x_i, y_i $ と 微分方程式そのものだけで十分であるということである。中間のややこしい値 を計算する必要はあるが、それは RK 法のほうが勝手にやることで、使う方が 入力を与える必要はない。

これに対して、一段階ではない方法、つまり多段階法というのは、 $y_{i+1}$ を計算するのに、 $y_i$ 以外の情報、例えば $y_{i-1}$ や、そこで計算した 導関数の値 $f_{i-1}$ 、さらにもっと昔の情報を使うやりかたのことである。 これは、プログラムとしてはもちろんRK法に比べれば面倒になる。昔の値をとっ ておかないといけないし、また、一番最初に計算を始める時にどうするかとい う問題もあるからである。

さらに、「一段階でない」というだけなので、可能な計算法にあまりに多様な 可能性がある。例えば、以下のような方法が考えられる。


\begin{displaymath}
y_{i+1} = y_{i-1} + 2 {\Delta x} f_i
\end{displaymath} (7)

これは、陽的中点公式といわれるもので、式としては偏微分方程式の解法の時 に空間微分の1次の項に使った中心差分と同じ形になっている。従って、理屈 としては2次の精度をもった公式になっている。

この方法は、しかし、立派な名前までついているのにもかかわらず、実は使え ない公式である。

どのような問題があるのかを示すために、以下の線形方程式

\begin{displaymath}
y = -ky
\end{displaymath} (8)

に陽的中点公式を適用したらなにが起きるかを考えてみる。刻みは固定の $\Delta x $ とすると、離散化したものは
\begin{displaymath}
y_{i+1} = y_{i-1}+2 -ky_i{\Delta x}
\end{displaymath} (9)

になる。これは線形差分方程式なので、前にもやったように固有値を調べれば いい。今 $\alpha = k\Delta x$ として整理すると、固有方程式は
\begin{displaymath}
\lambda^2 + 2\alpha\lambda -1 = 0
\end{displaymath} (10)

となって、これは実解を2つ持つ。それらを $\lambda_1, \lambda_2$ とすれ ば、 $\lambda_1 \lambda_2 = 1$ なので $\alpha \ne 0 $ ならどちらかは必 ず絶対値が 1 より大きい。

ちょっと式を見ればわかるように、絶対値が 1 より大きい固有値は $\alpha>0$、つまり $k>0$ なら負である。したがって、必ず振動的に発散す ることになる。

なお、上のような、安定な線形微分方程式について振舞いを調べるというのが、 常微分方程式の安定性解析の基本になる。非線形な方程式では違うとかいろん なことがあるかもしれないわけだが、まあ、少なくとも線形安定でないと話に ならないし、それ以上のことは一般論としていうのは難しいからである。

4.1 アダムス法

さて、中点公式はともかくとして、使える線形多段階法はじゃあどんなものか というわけだが、これも無限にいろんな作り方がある。そのへんの詳しい話は そういう本に譲ることにして、ここではもっとも広く使われているアダムス法 について説明する。

原理は、いくつかのステップでの導関数(微分方程 式の右辺) $f$の値を憶えておいて、それを通る補間多項式を作り、それを積 分して解を求めようというものである。

\begin{figure}\begin{center}
\leavevmode
\epsfxsize 8 cm
\epsffile{lmm.eps}\end{center}\end{figure}

上の図に概念を示す。ここでは、ラグランジュ補間 (ニュートン補間)をして多項式を作る。で、その作った多項式を積分する。 例えば、点 $i$ から $i+1$まで積分するのに、点 $i-p$ から $i$までの関数 値を使うとすれば、$p$次の多項式で

\begin{displaymath}
P(x_j) = f_j = f(x_j, y_j) \quad (i-p \le j \le i)
\end{displaymath} (11)

を満たすものを作る。で、$i+1$での解 $y_{i+1}$
\begin{displaymath}
y_{i+1} = y_i + \int_{x_i}^{x_{i+1}} P(x)dx
\end{displaymath} (12)

で与えられる。 刻み $h$が定数であるとすれば、$p$を決めれば上の式を
\begin{displaymath}
y_{i+1} = y_i + h\sum_{l=0}^{p}a_{pl}f_{i-l}
\end{displaymath} (13)

の形に書き直せる。

簡単な例として、$ p = 1$の場合を考えてみよう。この時、補間多項式は一次 であって

\begin{displaymath}
P(x) = f_i - \frac{f_{i-1}-f_i}{h}(x-x_i)
\end{displaymath} (14)

となって、これを積分すれば、結局
\begin{displaymath}
y_{i+1} = y_i + \frac{h}{2}(3f_i - f_{i-1})
\end{displaymath} (15)

となる。

一般に、アダムス法では任意段数の公式が構成でき、その次数は段数に等しい ことがわかっている。これは、ルンゲ・クッタなどに比べればはるかによい性 質をもっているということでもある。

4.2 出発公式

アダムス法はいくらでも高次の公式が作れ、計算量もあまり多くないというこ とがわかっているが、必ずしも広く使われているというわけでもない。その理 由はいろいろあるが、一つは、

「どうやって計算を始めるべきかよくわからない」

ということである。つまり、初期値問題としてはもちろん $x_0$ における $y_0 $しか知らないのに、多段階法ではその前の時刻での解が必要になるわけ である。これに対する対応策はいくつかあるが、時間刻み一定の場合には、基 本的にはルンゲ・クッタなどの別な方法で解を求めておくというやりかたが普 通である。

というわけで結局プログラムを書く手間が2倍以上になるというのが、多段階 法の実用上の問題である。

4.3 陰的アダムス法

さて、前に述べた公式では、補間多項式を陽的に求めた。すなわち、時刻 $i$ とそれより以前の値だけを使っていた。これに対し、陰的な補間多項式、つま り $t_{i+1}$での関数の値を使った公式というものも考えられる。 刻み $h$が定数であるとすれば、$p$を決めれば前と同様に

\begin{displaymath}
x_{i+1} = x_i + h\sum_{l=-1}^{p-1}b_{pl}f_{i-l}
\end{displaymath} (16)

の形に書けることになる。また手をぬいて $ p = 1$の場合を考えれば、これは単 に台形公式
\begin{displaymath}
x_{i+1} = x_i + \frac{h}{2}(f_i + f_{i+1})
\end{displaymath} (17)

となる。

陰的公式の場合、例によってどうやって代数方程式を解くかが問題になる。通 常の方法は、

というものである。ただし、特に線形多段階法の場合、反復を繰り返さないこ とが多い。反復回数を1回とか2回に固定してしまうのである。なぜそ れでいいか、また、そもそもなぜ陰的公式を使うかというあたりはレポート課 題ということで。

なお、このやり方を、予測子・修正子法と呼ぶ。線形多段階法はほとんどこの 形で使われるため、線形多段階法のことをさして予測子・修正子法と呼ぶ人も いる。


next up previous
Next: 5 構造体とクラス Up: 計算天文学 II 第5回 常微分方程式の初期値問題(1) Previous: 3 ルンゲ・クッタ以外の方法
Jun Makino
平成16年11月13日