Subsections

• 2.1 RK法の場合
• 2.2 線形多段階法の場合
• 2.3 ニュートン補間とアダムス公式の一般的導出

2 刻み幅調節と埋め込み型公式

ルンゲクッタにしても線形多段階法を使うにしても、実際に問題を解こうとい う上ではいろいろ考えないといけないことがある。もっとも重要なことは、刻 み幅 $h$ をどうやって決めるかということである。

実際問題としては、解いていくにつれて導関数の値が大きく変わるといったこ とは普通に起きる。このために、刻み幅が一定のままで計算していくというの は非常に無駄なことが多い。

無駄を減らすためには、局所誤差が小さい（解が滑らかな）ところは刻みを大 きく、逆に解が急にかわるところでは刻みを小さくしてやればよい。このよう な方法を可変刻み (variable step) 、あるいは適応刻み (adaptive step) と いう。

2.1 RK法の場合

例えば Runge-Kutta 法の場合、 $y_n$ から $y_{n+1}$ にいくのに前に計算 した情報とかなにかを使うわけではないので、刻み幅は任意にとることができ る。問題は、どうやって刻み幅を決めるかである。

これにはいろいろな考え方があるが、通常使われるのは、局所離散化誤差を推 定して、それがある値以下になるようにするという方法である。

このための一つの考え方は、以下のようなものである。

1. まず、ある幅 $h$ で積分する。
2. 次に、その同じ幅を、 $h/2$ ずつ 2回に分けて積分する。
3. すると、2回に分けて積分したほうがより正確な答になっていると期待 できるので、二つの解の差が誤差の推定値（大きめではあるが）となっている。

誤差の推定値が求まったら、普通は次のステップを調節する。局所離散化誤差 の次数がわかっているから、それに応じてステップサイズを調整すればよい。 次数プログラムによっては、現在のステップでの誤差が予定した値よりも大き ければステップを小さくして計算し直すようになっているものもある。

実装（プログラムを作ること）が容易であることもあって、この方法はわり と広く使われている。が、かなり無駄が多い方法であるというのも確かである。 というのは、単に誤差の推定のためだけに、 50% の無駄な計算をしているか らである。まあ、 50% はたいしたことないという考え方もあるが、計算に何 日も掛かるというような状況なら 50% は無視できないので、まあ、なんとか したくなるのが人情というものである。

というわけで、もうちょっとうまい方法はないかということで、以下のような ことを考えた人がいる。

一般に、 RK型の公式では、最終的な値は

$$y_{n+1} = y_n + \sum_{i=1}^s b_ik_i$$ (1)

の形になっている。ここで、 $k_i$ は全部そのまま使って、 $b_i$ を別の $b'_i$に置き換えた

$$y' = y_n + \sum_{i=1}^s b'_ik_i$$ (2)

で、局所離散化誤差が元の公式よりも大きいが、むやみと大きくはない（例え ば、一次次数が低い）ものがあれば、元の公式との差を誤差の推定に使うこと ができる。

この形の公式のことを埋め込み型 embedded 公式という。最初に提案した人の 名前をとって Runge-Kutta-Fehlberg 公式ということも多い。これもいろいろ 提案されているが、前回にも紹介した Dormand-Prince 公式はすべて埋め込み型であり、こ れが最近はもっとも広く使われている。精度が高いものが必要な時は、 8次の公式が使われる。これらは

http://www.unige.ch/math/folks/hairer/software.html

から入手できる (Fotran と C がある)。

2.2 線形多段階法の場合

線形多段階法の場合、 RK と違って誤差の推定は容易である。例えば、アダム ス法で予測子・修正子ペアで使うなら、予測子と修正子の差は誤差の（オーダー としては）よい推定値になっている。したがって、この部分については面倒な ことは特にない。

問題は、実際に刻みを変えるほうである。もっとも一般的な方法は、「その場 で公式を導出する」というものである。これについて以下簡単に説明する。

2.3 ニュートン補間とアダムス公式の一般的導出

アダムス法を使うには、要するに多項式補間の系数が出せればよい。以下、ニュー トン補間を使う方法を説明する。

補間したい関数 $f(x)$ が、標本点 $(x_0, x_1, ... x_n)$ で、関数値 $(f_0, f_1, ... f_n)$ をとるとする。これを、以下のような形の多項式

$$P(x) = a_0 + a_1(x-x_0) + a_2 (x-x_0)(x-x_1) ... + a_n(x-x_0) \cdots (x-x_{n-1})$$ (3)

という形に書くことを考える。これを、 $a_p$までをとったものは $x_0 \cdots x_{p}$ までを通る最低次補間多項式になっているように構成すること にする。で、低い次数から順に作っていく。付け加える新しい項は、それまで に使った点すべてで0になるので、新しい点で標本と一致するように $a_i$ の 値を決めればいい。

まず、$n=0$ については、 $a_0 = f_0$ とすればいいのは明らかであろう。 次に$n=1$ であるが、

$$P(x_1) = f_0 + a_1 (x_1 - x_0) = f_1$$ (4)

から

$$a_1 = \frac{f_1-f_0}{x_1 - x_0}$$ (5)

となる。同様に、 $a_2$を求めると、

$$a_2 = \frac{f[0,2]-f[0,1]}{x_2-x_1}$$ (6)

ただし、

$$f[i,j] = \frac{f_i - f_j}{x_i - x_j}$$ (7)

である。

さて、段々式が繁雑になるが、もうちょっとの辛抱である。上の形は、

$$a_2 = \frac{f[2,1]-f[1,0]}{x_2-x_0}$$ (8)

というふうにも書き直せることに注意しよう。これは、$(1,2)$を通る補間式 の一次の係数と、 $(0,1)$を通る補間式の係数の差みたいなものになっている。

ここで、天下りに、 $k$ 階差商 (divided difference) $D_{k,l} \quad (0\le l \le n-k)$ というものを導入する。これは以下のように定義される。

$$D_{0,l} = f_l \quad (0\le l \le n)$$ (9)

$$D_{k,l} = \frac{D_{k-1,l} - D_{k-1,l+1}}{x_l - x_{l+k}}$$ (10)

実は、このように定義された $D$ が、

$$D_{k,0} = a_k$$ (11)

を満たす、求める係数であることを示すことが出来る。

証明には、$D_{k,l}$が、 $x_l$から始まる $k+1$点を通る補間多項式の係数 であるということを用いる。定義により

$$D_{k,0} = \frac{D_{k-1,0} - D_{k-1,1}}{x_0 - x_{k}}$$ (12)

である。 ここで、帰納 法を使って証明することにすると、

$$P_{k-1}(x) = D_{0,0} + D_{1,0}(x-x_0) + \cdots + D_{k-1,0}(x-x_0)\cdots(x-x_{k-2})$$ (13)

は、 $x_0, \cdots x_{k-1}$ を通る $k-1$次補間多項式であり、また、

$$P_{k-1}'(x) = D_{0,1} + D_{1,1}(x-x_1) + \cdots + D_{k-1,1}(x-x_1)\cdots(x-x_{k-1})$$ (14)

は、 $x_1, \cdots x_{k}$ を通る $k-1$次補間多項式であるとしてよい。こ の2つの線形結合によって、 $x_0, \cdots x_{k}$ を通る補間多項式を作るに は、

$$P_k = \frac{(x-x_0)P_{k-1}' - (x-x_k)P_k}{x_k - x_0}$$ (15)

とすればよい。ここで、$D_{k-1,0}$と$D_{k-1,1}$は、それぞれ $P_{k-1}$と $P_{k-1}'$ の最高次の係数であることに注意すると、 $P_k$の最高次の係数(これは $a_k$に等しい) は、$D_{k,0}$で与えられることがわかる。

さて、

$$Q_k(x) = D_{0,0} + D_{1,0}(x-x_0) + \cdots + D_{k,0}(x-x_0)\cdots(x-x_{k-1})$$ (16)

という多項式を考えると、これは $x_0, \cdots x_{k-1}$ で $P_k$と一致し、 さらに最高次の係数も一致しているので、結局 $Q_k = P_k$ であることがわ かる。つまり、 $Q_k(x)$ が求める補間多項式であり、 $D_{k,0}$がその最高 次の係数ということになる。

と、こんな風に差商を作っておくと、これから多項式の係数を機械的に計算し て積分して、、、、ということができる。さらに、差商の形でデー タをとっておくと、上の漸化式を使って$D$を順番に更新することができる。 これは Krogh 型の公式と呼ばれ、かなり広い範囲の問題について、もっ とも効率が良い公式であることがわかっている。

なお、上の方法で時間刻みを変えられるなら、出発公式は不必要になることに 注意しておく。つまり、最初は1次の予測子、2次の修正子から出発して、ステッ プ毎に次数をあげていけばいい。ステップサイズは次数を考慮して決めておく。