常微分方程式の初期値問題の数値解法
桂田 祐史

1994 年 6 月〜2011 年 5 月 21 日, 2022 年 6 月 19 日
(なかなか整理するための時間が取れない。我ながらひどい出来だと思うのだが…正直に白状すると、私は良く理解できていない。
この文書と紛らわしい「常微分方程式の初期値問題の数値解法入門」 [1] という題の文書も用意した。そちらは古い授業資料を元にしていて、やさしい微分方程式の問題をとにかく解いてみよう、という内容である。)

目次
1

序

2

2

常微分方程式の初期値問題の復習
2.1 数学理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 簡単な数値解法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2.1 前進 Euler 法 (forward Euler’s rule) . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 後退 Euler 法 (backward Euler’s rule) . . . . . . . . . . . . . . . . . . . . . . .
2.2.3 Runge-Kutta 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3
3
4
5
5
5

3

基本的な概念・用語
3.1 多段法, 段数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 局所離散化誤差、公式の次数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 前進 Euler 法の収束証明 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5
6
6
8

4

典型的なスキーム (1) Runge-Kutta 法とその一族
4.1 歴史 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 定義と Stetter の行列表現 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 収束定理 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 特徴 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 次数と段数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 前進型公式の考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6.1 1 段 1 次 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6.2 2 段 2 次 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6.3 3 段 3 次 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6.4 4 段 4 次 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7 埋め込み型の公式, RKF45 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7.1 RKF45 公式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.7.2 RKF45 による刻み幅の自動調節 (書き直し版、工事中) . . . . . . . . . . . . .

9
9
10
11
11
11
11
11
12
12
12
13
13
14

1

実験例: 爆発する問題を RKF45 で計算 . . . . . . . . . . . . . . . . . . . . . .
RKF45 自作プログラム開発 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15
18

5

典型的なスキーム (2) 線形多段法
5.1 PC (予測子修正子法) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21
21

6

その他の方法
6.1 Taylor 法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2 補外法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22
22
22

7

数値的安定性
7.1 線形安定性解析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23
23

8

Stiﬀ problem (硬い問題)

25

9

参考書

26

4.7.3
4.7.4

10 おまけ — 実際的な誤差の推測

27

A 数値計算するための情報
A.1 はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2 C 言語+グラフィックス・ライブラリィ . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3 C 言語+gnuplot . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.4 Java . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.5 C++ の利点 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.6 Ruby . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27
27
28
28
28
29
29

B 眠りから覚めてみたら (2022/6)

30

序

1

この文書では、1 階正規形の常微分方程式の初期値問題に対する数値解法を扱う。すなわち

(1)
(2)

dx
= f (t, x) (t ∈ I)
dt
x(t0 ) = x0を満たす x = x(t) の近似解を求めることを考える。ここで I は t0 ∈ R を含む R の区間で、(
f : Rn+1 ⊃ Ω(開集合) → Rn 連続,(t0 , x0 ) ∈ Ωは与えられているとする。
注意 1.1 高階の方程式も 1 階にできることが多いので、1 階正規形の方程式は十分一般的であると考えられる。
微分方程式の解は関数であり、これは関数空間の要素としてとらえるのが (現代の数学では) 普通である。(大抵の場合、関数空間は無限次元空間で、問題を難しくしている。)

2

微分方程式は解析的1 に解けないことが多い。たとえ解けても便利でないことがある2 。
近似解法では、解の有限的な近似表現を求める。具体的には、• 「有限次元の関数空間の要素で近似する」が基本。
• 特に連続変数を離散変数に置き換えて近似する離散変数法が有力。
残念ながら常微分方程式の初期値問題に限っても万能の方法はない。
プロでない平均的ユーザーとしては、実際的にはとりあえず Runge-Kutta 法を使い3 、不満があれば他の方法を考える、くらいで良いだろう。
この講義では、基礎概念を簡単に説明した後で、

1. 刻み幅の自動調節 (adaptive stepsize control)
2. 硬い方程式 (stiﬀ problem)などの話題を紹介する。詳しいことを知りたい場合は、まず三井 [2] を見るとよい。
最近 (2004 年 2 月)、面白い本が出版された。三井・小藤・齊藤 [3] である。第 2 章「ハミルトン系の解法」, 第 3 章「遅延微分方程式の解法」, 第 2 章「確率微分方程式の解法」と章の名前を見れば一目瞭然、現在盛んに研究されている分野への入門ができる (ヒットだと思う)。

(「最近」が 7 年前か…) 定番本の翻訳が出た。ハイラー・ネルセット・ヴァンナー [4] とハイラー・
ヴァンナー [5] である。現時点での決定版か。

2

常微分方程式の初期値問題の復習

2.1

数学理論

1 時間くらいでさらっと説明するための説明は、桂田・佐藤 [6] に書いた。
• 局所解の存在を保証するには、f の連続性を仮定するだけで十分4 。
• 解の一意性を保証するには f の連続性だけでは不十分。
例 2.1 (無限個の解の分岐) 次の常微分方程式の初期値問題には無限個の解が存在する。
(
x′ = x1/2 (t > 0)x(0) = 0.
実際 ∀t0 ≥ 0 に対して


 0(t ≤ t0 )def.
2x(t) =
 (t − t0 ) (t > t0 )4は解である。
1「解析的に解く」とは、不定積分を取る、四則演算を施す、逆関数を取る、初等関数に代入するなどで求めたり —いわゆる求積法で解いたり — 、解を級数で表現したりすることを指す。
2
例えば無限級数で解を表したとして、その値を計算するのはそんなに簡単ではない。無限級数は有限的な表現とは言えない。
3
実は実際に色々解いてみて、最近考え方が変ってきた。やはり最低限 RKF45 くらいは使いたい。しっかり解説を書
くべきだと思う。
4
ただしこれは有限次元の場合で、無限次元の場合は連続性だけではダメである (反例がある)。連続の場合の存在証明
はいわゆるコンパクト性の議論を用いるので、空間の次元が有限次元であることは本質的な仮定となってしまうのは自然である。一方、「f = f (t, x) が連続かつ x につき局所 Lipschitz 条件を満たすならば、局所解が一意的に存在する」という定理は無限次元でもそのまま成り立つ。

3

• f が次に示す “変数 x に関する局所 Lipschitz 条件” を満たせば一意性が成り立つ。∀(t0 , x0 ) ∈Ω, ∃V : (t0 , x0 ) の近傍, ∃LV > 0 s.t.
∥f (t, x1 ) − f (t, x2 )∥ ≤ LV ∥x1 − x2 ∥ ((t, x1 ), (t, x2 ) ∈ V ).
f がこの条件を満たすための分かりやすい十分条件として、f が C 1 級であることがあげられる。すなわちf ∈ C 1 (Ω) =⇒ f : 局所 Lipschitz.
• 大域的な存在について。(t, x(t)) → ∂Ω または ∥x(t)∥ → ∞ となるまで左右に延長できる (延長不能解の存在定理)。
このうち、(t, x(t)) が f の定義域 Ω の境界に近付くという条件は分かりやすいが、∥x(t)∥ → ∞の方は見慣れない人もいるかも知れない。これについては次の例を見るとよい。
例 2.2 (爆発解)

(

x′
= x2
x(0) = 1の解はx(t) =

1
1−t(t ∈ (−∞, 1))

であり、

lim x(t) = ∞.
t↑1

• 次の一様 Lipshitz 条件を満たせば、爆発しない。
∃L > 0 ∀(t, x1 ), (t, x2 ) ∈ Ω

2.2∥f (t, x1 ) − f (t, x2 )∥ ≤ L∥x2 − x1 ∥.

簡単な数値解法(コンピューター・プログラミングの演習で、必ずと言って良いほど出会うポピュラーな数値解法を復習しよう。)
I = [a, b] とする。N ∈ N に対し、I を N 個の小区間に分ける:
a = t0 < t1 < t2 < · · · < tN = b.
このとき、各 tj における x の値 x(tj ) の近似値 xj を求めることを考える方法を離散変数法 (discretevariable method) と呼ぶ。hj := tj+1 − tj (j = 0, 1, · · · , N − 1) を刻み幅と呼ぶ。
区間の分割の仕方としては、例えば

hj ≡ h i.e. h =

b−a
,
N

tj = a + jh (j = 0, 1, · · · , N ).

のように等分割することが多い。以下この小節ではそれを仮定して説明するが、可変刻み幅に一般化することは容易である。

4

2.2.1

前進 Euler 法 (forward Euler’s rule)

微分係数 x′ (t) を前進差分商

x(t + h) − x(t)hで近似して作った漸化式

xj+1 = xj + hf (tj , xj ) (j = 0, 1, 2, · · · )で {xj }N
j=0 を計算する。
数学のどこにでも出て来る巨人 Leonhard Euler (1707–1783, スイスの Basel に生まれ、ロシアの St Peters-

burg にて没する) による。

2.2.2

後退 Euler 法 (backward Euler’s rule)

前進差分商のかわりに後退差分商

x(t + h) − x(t)hで近似して得られる漸化式

xj+1 = xj + hf (tj+1 , xj+1 ) (j = 0, 1, 2, · · · )で {xj }N
j=0 を計算する。
xj+1 を求めるために、方程式を解く必要がある (f によっては計算が面倒になることもある)。こ
ういう方法を一般に陰解法 (implicit method) と呼び、前進 Euler 法のように方程式を解かずに、求める値が直接計算できる方法を一般に陽解法 (explicit method) と呼ぶ。

2.2.3

Runge-Kutta 法

漸化式



k1 = hf (tj , xj )
 k = hf (t + h/2, x + k /2)2
jj1
 k3 = hf (tj + h/2, xj + k2 /2)
k4 = hf (tj + h, xj + k3 )1
(k1 + 2k2 + 2k3 + k4 )6で {xj }N
j=0 を計算する方法を (古典的、あるいは 4 次の) Runge-Kutta 法と呼ぶ。
xj+1 = xj +3

基本的な概念・用語
なるべく一般的な解法で説明しよう。簡単のため、刻み幅は一定である (hj ≡ h) とする。

5

3.1

多段法, 段数k 段法とは xj+k を定めるために
xj+k = a0 xj + a1 xj+1 + · · · + ak−1 xj+k−1 + hΦ(tj , tj+1 , · · · , tj+k , xj , xj+1 , · · · , xj+k )(3)≡ L(tj , tj+1 , · · · , tj+k , xj , xj+1 , · · · , xj+k , h)のように xj , xj+1 , · · · , xj+k を含んだ方程式を用いる数値解法のことである。
このような方程式のことをスキーム (scheme) と呼ぶ。
k−1
X
ここで a0 , · · · , ak−1 は
ai = 1 を満たす定数であり、Φ は f によって定まる、微分・積分など
i=0の無限小演算を含まない写像で、f ≡ 0 ならば Φ ≡ 0 となるものである。この (k, {ai }k−1i=0 , Φ) が一つの方法を特徴づける。
問

Euler 法、Runge-Kutta 法がこの形になっていることを確かめよ。
• 上の整数 k をスキームの段数 (step number) と呼ぶ。
• Euler 法、Runge-Kutta 法では、段数 k = 1 である。このときはL(tj , tj+1 , xj , xj+1 , h) = xj + hΦ(tj , tj+1 , xj , xj+1 )という形、つまり

xj+1 = xj + hΦ(tj , tj+1 , xj , xj+1 )というスキームになる。

• k ≥ 2 なるスキームを多段法 (multistep method) と呼ぶ。
• Φ が xj+k によらないように表される時、陽解法 (explicit method) であると呼び、そうでない場合を陰解法 (implicit method) と呼ぶ。陰解法では xj+k を求めるために、(一般には非線型の) 方程式を解かねばならないので、ほとんどの場合に反復法が必要になり、面倒である
が、次数が高くて安定性のよい方法が作れる。

3.2

局所離散化誤差、公式の次数

解法 (3) の、t における局所離散化誤差 (local truncation error) をτ (t, h) :=1[x(t + kh) − L(t, t + h, · · · , t + kh, x(t), x(t + h), · · · , x(t + kh))]h

で、また大域的離散化誤差 (global discretization error, global truncation error) をτ (h) :=

max |τ (t, h)|t∈[a,b−kh]で定義する。
これを用いて公式の次数 (order, 位数とも呼ぶ) を次のように定義する。
公式の次数が少なくとも m

def.

⇔

C m 級の一意解を持つ任意の初期値問題に適用した場合
τ (h) = O(hm ) (h → 0).
6

例 3.1 (ポピュラーな公式について) 前進 Euler 法, 後退 Euler 法は共に 1 段 1 次の公式であり、古典的 Runge-Kutta 法は 4 段 4 次の公式である。
注意 3.2 f があまり滑らかでないときなど、解がなめらかでない場合、次数 m の公式を用いてもτ (h) = O(hm ) は期待できない。
あらく言って m 次の公式とは、Taylor 展開して考えたとき、m 次の項まで一致するものであって、次のような性質を持つ:(i) (どういうわけか運良く) 第 j ステップまで誤差なく計算出来たとすると
x(tj+1 ) − xj+1 = O(hm+1 ) (h → 0).
(ii) 実際は誤差が累積するのでx(tN ) − xN = O(hm ) (h → 0).
この左辺を全離散化誤差 (total discretization error) と呼ぶ。
例 3.3 Euler 法は 1 次、古典的 Runge-Kutta 法は 4 次の公式である。そこで例えば x′ (t) = x(t ∈ (0, 1)), x(0) = 1 という初期値問題に適用した場合の累積誤差を表示したものが次の図である(横軸は区間の分割数 N , 縦軸は累積打ち切り誤差で、いずれも対数目盛)。
1

0.01

0.001

0.0001

1e-05

.1

1e-06

1e-07

1e-08

1e-09

0.01

1e-10

1e-11

1e-12

1e-13

0.001
1

10

100

1e-14

1000

1

10

100

1000

(左が Euler 法によるもの、右が Runge-Kutta 法によるもので、傾きがそれぞれ −1, −4 であることが分かる。)
実はこの問題の場合、 Euler 法では xn+1 = (1 + h)xn , Runge-Kutta 法では xn+1 = (1 + h + h2 /2 +h3 /3! + h4 /4!) xn となる。x(t + h) = eh x(t) = (1 + h + h2 /2 + h3 /3! + h4 /4! + · · · + hn /n! + · · · ) x(t)であるから、Euler 法は 1 次の項まで、Runge-Kutta 法は 4 次の項まであっていると言える。


収束のための条件

以下の三条件が成り立つならば、h → 0 とするとき、近似解が真の解に収束することが証明できる。

(i) 公式の次数が少なくとも 1 以上 (適合条件 (consistency) を満たす)。
(ii) Φ が xj , xj+1 , · · · , xj+k について Lipschitz 条件を満たす。
(iii) 初期値 x1 , x2 , · · · , xk−1 が適切に用意される。



7



3.3

前進 Euler 法の収束証明

もっとも簡単な前進 Euler 法の場合に近似解が厳密解に収束することを証明してみよう。



定理 3.4 (Euler 法の収束) 常微分方程式の初期値問題

dx= f (t, x) (t ∈ (a, b)),dtx(a) = x0

において、f は連続で、Lipschitz 条件∥f (t, x1 ) − f (t, x2 )∥ ≤ L∥x1 − x2 ∥を満たし、C 2 級の解が存在すると仮定するとき、Euler 法による近似解は微分方程式の解に収束
する。




x が C 級であるから、∀h ∃θ ∈ (0, 1) s.t.
2

1
x(t + h) = x(t) + x′ (t)h + x′′ (t + θh)h2 .
2x′ (t) = f (t, x(t)) であることに注意して、t = tn , t + h = tn+1 とおくと、
1x(tn+1 ) = x(tn ) + f (tn , x(tn ))h + x′′ (tn + θh)h2 .
2Euler 法の公式より
xn+1 = xn + f (tn , xn )hであるから、辺々引き算して(4)

1
x(tn+1 ) − xn+1 = x(tn ) − xn + (f (tn , x(tn )) − f (tn , xn )) h + x′′ (tn + θh)h2 .
2

これから、もしも x(tn ) = xn であったとすると、

1x(tn+1 ) − xn+1 = x′′ (tn + θh)h22を得る。これから Euler 法の局所離散化誤差は O(h) である (ゆえに Euler 法の次数は 1 である) ことが分かる。
さて、x が C 2 級であることから
M := max ∥x′′ (s)∥s∈[a,b]が存在する。
(4) から∥x(tn+1 ) − xn+1 ∥ ≤ ∥x(tn ) − xn ∥ + ∥f (tn , x(tn )) − f (tn , xn )∥ · |h| +≤ ∥x(tn ) − xn ∥ + L ∥x(tn ) − xn ∥ · |h| +1= (1 + L|h|) ∥x(tn ) − xn ∥ + M h2 .
2
ここで

en := ∥x(tn ) − xn ∥81
· M h22

1 ′′∥x (tn + θh)∥ h22

とおくと、1
en+1 ≤ (1 + L|h|)en + M h2 .
2e0 = 0, n|h| ≤ b − a に注意して、以下の補題 3.5, 3.6 を用いると(1 + L|h|)n − 1 1
M |h|· M h2 =[(1 + L|h|)n − 1]
L|h|
2
2L
M |h| nL|h|
M |h| L(b−a)≤(e− 1) ≤
(e
− 1).
2L
2L

en ≤





補題 3.5 数列 {en }n≥0 が、実定数 A, B (ただし A ≥ 0) に対して漸化不等式

en+1 ≤ Aen + B(n ≥ 0)を満たすならばen ≤ An e0 + (An−1 + An−2 + · · · + A2 + A + 1)B.
特に A ̸= 1 ならば

e n ≤ An e 0 +証明

An − 1B.
A−1

順に代入していくだけである。

e1 ≤ Ae0 + B,
e2 ≤ Ae1 + B ≤ A(Ae0 + B) + B = A2 e0 + (A + 1)B,e3 ≤ Ae2 + B ≤ A(A2 e0 + (A + 1)B) + B = A3 e0 + (A2 + A + 1)B
より明らかにen ≤ An e0 + (An−1 + An−2 + · · · + A + 1)B.
(なんか、Gronwall の数列版か？ )

補題 3.6 x > 0 に対して

(1 + x)1/x < e.


f (x) = log(1 + x) について、テイラーの定理より ∀x > 0, ∃θ ∈ (0, 1) s.t.

証明

1
1
−1
log(1 + x) = f (x) = f (0) + f ′ (0)x + f ′′ (θx)x2 = 0 + 1 · x + ·x2 < x.
2
2 (1 + θx)2
ゆえにlog(1 + x)1/x =log(1 + x)
< 1.
x

これから (1 + x)1/x < e を得る。

4

典型的なスキーム (1) Runge-Kutta 法とその一族

4.1

歴史

以下の説明は、一松 [7] による。

9

Carl Runge (1856–1927, ドイツの Bremen に生まれ、Göttingen にて没する) が 4 次の Runge-Kutta
公式を導いた (1895)。Heun (1900) と Martin Wilhelm Kutta (1867–1944, Upper Silesia ( 現在のポーランド) に生まれ、ドイツの Fürstenfeldbruck にて没する) (1901) は一般的に5 研究した。1940 年代の Gill の研究もあったが (Runge-Kutta-Gill の公式はかつては有名だった)、公式の本格的な再検討が始まったのは、Ceschino, Butcher, 田中正次等の 1960 年頃からの研究による。

4.2

定義と Stetter の行列表現

前節に解説したタイプの公式のうち、k = 1 の場合、すなわち 1 段法を Runge-Kutta 型公式という。いくつかの (t, x) について、右辺の f (t, x) を計算し、それらの重みつき平均によって、適当な次数の (=その次数までの真の解の Taylor 展開と一致するような) 公式を作っている。具体的には
Runge-Kutta 型公式の一般形は

s
X



x
=
x
+
h
µi k i

j
 j+1
i=1
!
s
X



= f ti + αi h, xj + h(i = 1, · · · , s)βiℓ kℓ kiℓ=1の形に書くことが出来る。ここで s を段数 (number of stages) と呼ぶ。段数とは、要するに 1 ス
テップ先に進めるために必要な f の計算回数である6 。
上の式は、k1 , k2 , · · · , ks についての s 元連立 1 次方程式を解いて、その重みつき平均を xj に足して xj+1 を求める、という手順で使うことになる。
段数 s と係数 {αi ; 1 ≤ i ≤ s}, {βij ; 1 ≤ i, j ≤ s}, {µi ; 1 ≤ i ≤ s} を選ぶとスキームが定まることになる。それら係数を並べた
α1 β11 · · · β1sα2 β21 · · · β2s.
..
..
. · · · ..
.

αsβs1 · · · βssµ1 · · · µsのような表を Stetter の行列表現 (Stetter’s notation) と呼ぶ7 。
公式が前進型 (陽的、explicit) であるとは

i≤j

=⇒

βij = 0が成り立つことを言う8 。このとき、k1 , k2 , · · · , ks の順に容易に計算できる。
公式が前進型でない場合、公式が陰的 (implicit) であるという。陰的な場合でも

i<j

=⇒

βij = 0が成り立つ9 場合は半陰的 (semi-implicit) であるという。
5

例えば、4 段 4 次の公式にはどういうものがありうるかとか、5 段 5 次の公式は存在しないようだとか。
1 段法なのに「段数」とは？ 英語では step, stage と区別があるのだが、日本語では同じ「段」となってしまって紛らわしい。残念ながら、この用語はもう定着してしまっている。
s
s
X
X
7
µi = 1 のような条件があるのだと思うが…
多分 αi =
βiℓ ,
6

8

ℓ=1

i=1

行列 (βij ) の対角線の上側と対角線上にある成分がすべて 0 であること。つまり (βij ) が狭義の下三角行列であるという条件である。
9(βij ) の対角線の上側にある成分がすべて 0 (対角成分自身は 0 でなくてもよい) ということ。つまり (βij ) が狭義下三角行列であるという条件である。

10

収束定理

4.3

(工事中)





定理 4.1 (Runge-Kutta 型公式の収束のための条件) Runge-Kutta 型公式が収束するためには
次の二条件が成立することが必要十分。
′(i) 適合性 (ρ(1) = 0, ρ (1) =

s
X

βj )

j=0

(ii) 安定


4.4



特徴

Runge-Kutta 型公式の特徴として、次のものがあげられる。
(1) 自己出発的 (self-starting) である。すなわち、多段法 (k ≥ 2) では計算の最初に必要となる x1 ,x2 , · · · , xk−1 を準備することなく、計算が開始できる。
(2) 計算の途中で刻み幅 (stepsize) h の変更が簡単である (→ adaptive stepsize control に便利)。

4.5

次数と段数

次数を大きくしようとすると、普通は f の値の計算回数 (段数, stage) が増えることになる。
公式の範囲を限定して、段数 s を与えたとき、可能な最高の次数 m を、到達可能次数 (到達可能位数) と呼ぶ (Butcher による)。
m の下からの評価は、少なくとも一つの m 次公式を作ることによって得られるが、m の上からの評価は s 元連立代数方程式の解の不存在証明なので、s が大きくなると急激に困難になる。
陽的 Runge-Kutta 型公式の場合には、s ≤ 8 までの場合に到達可能次数が分かっているという (一松 [7]10 )。


到達可能次数

段数 s
到達可能次数 m

1
1

2
2

3
3

4
4

6
5

7
6



9 10
7 8





ここで m ≥ 5 のとき s > m となることに注意しよう (この事実が 4 次の Runge-Kutta 法が人気のある理由の一つである)。
なお陰的 Runge-Kutta 型公式では、s 段で 2s 次を到達する公式が存在することが分かっている。

4.6

前進型公式の考察

4.6.1

1段1次

s = 1 の場合は、公式は前進 Euler 法だけしかない。

10

今ではこの本も古くなったので、もっと大きな s まで分かっているかも知れない。また別の文献に書いてあることと矛盾していたような覚えがある。

11

4.6.22段2次

s = 2, q = 2 の公式は Heun 法と総称される。
β11 = β12 = β22 = 0,となる。β21 =: β とおくと、µ1 = 1 −µ1 + µ2 = 1,

1
,
2β

µ2 =

β21 µ2 =

1
2

1
2β

となり、



k1 = hf (ti , xj ),k2 = f (ti + βh, xj + βk1 ) ,

1
のとき、
2


h
k1
k2 = hf ti + , xj +,2
2

1
xj+1 = xj + 1 −
2β


k1 +

1
k2 .
2β

特に β =

xj+1 = xj + k2

(中点公式, 全体を改良 Euler 法).

また β = 1 のとき、

1
xj+1 = xj + (k1 + k2 )2

k2 = hf (ti + h, xj + k1 ),4.6.3(台形公式, 全体を修正 Euler 法).

3段3次(省略する。一松 [7] を見よ。)

4.6.4

4段4次

以下は Kutta (1901) の研究だそうである (一松 [7] からの孫引き)。
前提条件
β11 = β12 = β13 = β14 = β22 = β23 = β24 = β33 = β34 = β44 = 0のもとで、パラメーターを

β21 =: α,β32 =: λ,β31 + β32 = β,β42 =: ν,とおくと、次の Kutta の条件式をえる。

µ1 + µ2 + µ3 + µ4 = 1,







1

 α 2 µ2 + β 2 µ3 + γ 2 µ4 = ,

3

1


αλµ3 + (αν + βµ)µ4 = ,


6






 α2 λµ3 + (α2 ν + β 2 µ)µ4 = 1 ,12

β43 =: µ,β41 + β42 + β43 =: γ1αµ2 + βµ3 + γµ4 = ,21
α 3 µ2 + β 3 µ3 + γ 3 µ4 = ,41
αβλµ3 + (αν + βµ)γµ4 = ,81
αλµµ4 = .
24

この一般解も知られていて (自由度 2 が残る)、そのうちすべての係数が ≥ 0 で、0 < α ≤ β ≤ γ ≤ 1
である「単調な」公式は、次の古典的 Runge-Kutta 公式 (Runge の原公式, 本来の Runge-Kutta 公式) しかない。

12





本来の Runge-Kutta 公式

Kutta のパラメーターで、
1
α=λ=β= ,2
すなわち Stetter の行列表現で



0
 1

 2

 0


 0

1
6

ν = 0,

0

0

0
1
2
0
1
3

0
0
1
1
3

µ = 1,0

γ = 1,




0 


0 
,

0 

1
6

スキームで

k1 = hf (ti , xj ),k2 = hf (ti + h/2, xj + k1 /2),
k3 = hf (ti + h/2, xj + k2 /2),
k4 = hf (ti + h, xj + k3 ),1
xj+1 = xj + (k1 + 2k2 + 2k3 + k4 ).
6



4.7

埋め込み型の公式, RKF45

4.7.1

RKF45 公式

E. Fehlberg は次の公式 RKF4511 を提案した (1970, [8])。


16
6656
28561
9
2
(5)
xj+1 = xj + h
k1 +
k3 +
k4 − k5 + k6 ,
135
12825
56430
50
55
ただし

k1 = f (tj , xj )11
k2 = f tj + h, xj + hk1
4
4


1
3
k3 = f tj + h, xj + h(3k1 + 9k2 )
8
32


12
1
k4 = f tj + h, xj +h(1932k1 − 7200k2 + 7296k3 )
13
2197



439
3680
845
k5 = f tj + h, xj + h
k1 − 8k2 +
k3 −
k4
216
513
4104



1
8
3544
1859
11
k6 = f tj + h, xj + h − k1 + 2k2 −k3 +k4 − k5
.
2
27
2565
4104
40
11

この公式は御覧の通り大変複雑である。筆者はこの公式で誤植のドジを犯したことがある (レポート課題のプリントで間違えた — 罪重…)。この公式を用いたプログラムを書く時は複数の資料でチェックすることをお勧めする。なお、オリジナルは E.Fehlberg; “Klassische Runge-Kutta-Formeln vierter und niedrigerer Ordnung mit Schrittweiten-Kontrolle
und ihre Anwendung auf Wrmeleitungs-probleme”, Computing, Vol.6, pp.61–71 (1970).

13

これは 6 段 5 次の公式であるが、それだけでなく


25
1408
2197
1
∗
(6)
xj+1 = xj + h
k1 +
k3 +
k4 − k5
216
2565
4104
5
という値を作ると、x∗j+1 は x(tj+1 ) に対して 4 次の近似値となる。この次数の差を利用してステッ
プ幅の自動調節をする方法を以下に説明する。
このアイディアのキーは、 s 段 m 次公式に、f の値を計算することなく (m − 1) 次公式を付随させるところにある。このような Runge-Kutta 型公式を、(m − 1) 次公式が m 次公式に埋め込まれているといい、埋め込み型 Runge-Kutta 法と呼ぶ。

4.7.2

RKF45 による刻み幅の自動調節 (書き直し版、工事中)

t = a から t = b まで、誤差の大きさの見積が、計算者が指定した値 ε より小さくなるように解くことを目標にする。このとき ε のことを許容誤差限界と呼ぶ。[a, b] を a = t0 < t1 < · · · < tN = b と分割して解くことにする。t = tj における値 xj まで定まったとする。tj+1 における 5 次公式, 4 次
公式による近似値をそれぞれ xj+1 , x∗j+1 とすると、
x∗j+1 − x(tj+1 ) = c′ h5 + O(h6 )xj+1 − x(tj+1 ) = ch6 + O(h7 ),となる。ただし h := tj+1 − tj とおいた。これから

xj+1 − x∗j+1 = −c′ h5 + O(h6 ).
通常は、xj+1 は x∗j+1 よりも格段に精度が良いと期待出来るので、この式の値が x∗j+1 の誤差の良い評価となっていると考えられる:
xj+1 − x∗j+1 = (xj+1 − x(tj+1 )) + x(tj+1 ) − x∗j+1 ≒ x(tj+1 ) − x∗j+1 .
さて、δj+1 := xj+1 − x∗j+1と置いておく。
目標は [a, b] で解いたときの許容誤差限界を ε とすることであるから、[tj , tj+1 ] で解いたときの許容誤差限界はhε · , H := b − aHとするのが妥当であろう。それゆえ
εhδj+1 ≤
H
であれば良いが、そうでない場合は、h が大きすぎて十分な精度が得られていないと考え、もっと hを小さくすることにする。h の代りに、h ≪ h なる h を用いて

tj+1 = tj + hとおき、これに対応した xj+1 , x∗j+1 を求め、δ j+1 := xj+1 − x∗j+1 とおく。

δj+1 ∼ ∥c′ ∥ |h|5 ,δ j+1 ∼ ∥c′ ∥ hとなると期待できるから、

δ j+1 ∼

h
|h|
14

!5
δj+1 .

5

今度は δ j+1 ≤ εh/H となって欲しいわけだが、この不等式の左辺に推定値を代入した不等式

h
|h|

!5

を h について解くと

δj+1 ≤
h≤hを得る。安全のために

εh
Hδj+1


h := 0.9h

εh
H

1/4

εh
Hδj+1

1/4

で h を定めることにする。

4.7.3

実験例: 爆発する問題を RKF45 で計算(かなりいい加減なプログラムであるが、あえてここに含める理由の一つは、公式中の係数に書き間違いがないことを確認できるようにするためである。)/** rkf.c --- RKF45 のサンプル (まだ未完成です)*/#include <math.h>#define RKF45_STAGE 6
static double rkf45_alpha[RKF45_STAGE] = {0.0, 1.0/4, 3.0/8, 12.0/13, 1.0, 1.0/2};
static double rkf45_beta[RKF45_STAGE][RKF45_STAGE-1] = {
{
0.0,
0.0,
0.0,
0.0},
{
1.0/4,
0.0,
0.0,
0.0},
{
3.0/32,
9.0/32,0.0,
0.0},{1932.0/2197, -7200.0/2197, 7296.0/2197,
0.0,
0.0},
{ 439.0/216,-8.0,
3680.0/513, -845.0/4104,
0.0},{ -8.0/27,
2.0,-3544.0/2565, 1859.0/4104, -11.0/40}};
static double rkf45_mu[RKF45_STAGE] = {16.0/135, 0.0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55};
static double rkf45_mu2[RKF45_STAGE] = {25.0/216, 0.0, 1408.0/2565, 2197.0/4104, -1.0/5};
double sum(int n, double *x){
int i;
double s = 0;
for (i = 0; i < n; i++) s += x[i];return s;}/* 簡単なチェック */void check_table(){15

int i;for (i = 0; i < RKF45_STAGE; i++)printf("%g %g\n", rkf45_alpha[i], sum(i, rkf45_beta[i]));printf("rkf45_mu の和=%g\n", sum(RKF45_STAGE, rkf45_mu));printf("rkf45_mu2 の和=%g\n", sum(RKF45_STAGE, rkf45_mu2));}
double dotprod(int n, double *x, double *y){
int i;
double s = 0;
for (i = 0; i < n; i++) s += x[i] * y[i];return s;}/** RKF45 で時刻 *t から、1 ステップ進む* 特に問題なければ *hh だけ進むが、単位長さあたりの許容誤差限界 eps_tol* (解説のε/H に相当) を達成するために、必要ならば刻み幅の自動調節をする。
*/
int rkf45(double *newx,
double *t, double x, double *hh, double hmin, double eps_tol,
double f(double, double)){
int i, s = RKF45_STAGE;
double k[RKF45_STAGE], newx1, newx2, h = *hh, error;
do {/* k1,k2,k3,k4,k5,k6 を計算する */for (i = 0; i < s; i++)k[i] = h * f(*t + rkf45_alpha[i] * h, x + dotprod(i, rkf45_beta[i], k));/* 5 次精度, 4 次精度の値を計算 */newx1 = x + dotprod(s, rkf45_mu, k);newx2 = x + dotprod(s, rkf45_mu2, k);/* 見積誤差 */error = fabs(newx1 - newx2);
if (error < eps_tol) {*newx = newx2;*hh = h;*t += h;return 0;}// 新 h = 0.8 h (ε h/ (H 誤差見積)^(1/4)
h = 0.8 * h * pow(eps_tol * h / error, 0.25);}
while (h > hmin);printf("h=%g<=hmin=%g\n", h, hmin);return 1;}
double f(double t, double x){
return x * x;}
int main(){
int n = 1000;
double a, b, x, h, hmin, t;
check_table();a = 0; b = 1; h = (b - a) / n; hmin = h / 100000;
x = 1;
t = 0;
while (t < 1.0) {16

if (t + h > 1.0) h = 1.0 - t;
if (rkf45(&x, &t, x, &h, hmin, 1e-13, f) != 0) {fprintf(stderr, "誤差が小さくならない\n");exit(1);}
printf("t=%25.16f, x=%25.16g\n", t, x);}
printf("%25.16g\n", x);
return 0;}

このプログラムでは爆発の例に出した初期値問題

dx= x2 ,
dtx(0) = 1

を解いているが、(実行してみると12 ) 爆発時刻 t = 1 のところで、刻み幅がどんどん小さくなってい
くのが分かる。
なお、具体的な時間刻み幅の制御について具体的に書いてある本は多くない (そのためもあって、私自身よく分かっていない部分もあり、上記のプログラムも完成品ではない)。その点で森 [9] は貴重である。そこに掲載されているプログラムのアルゴリズムは (FORTRAN プログラムを C 風に書き直して示せば) 以下のようなものである (と思う。読み違いが無ければ。)。なるほどと思うところもあるが、今一つ納得できていないところもある。

maxiter = 100;
epsmac = 1e-15;// t0 から tn まで解く// eps は許容絶対誤差、epsmac は計算機イプシロン程度の数
epsv = max(eps, epsmac); // はて、相対誤差ならともかく、絶対誤差で筋が通る?
// eb は単位時間当り許される誤差
eb = fabs(epsv / (tn - t0));//
t = t0;h = tn - t0;//
for (iter = 1; iter <= itermax; iter++) {// とにかく刻み幅 h で進んでみて、局所打ち切り誤差の推定もする
rkf(t, h, x0, &xn, &errmax);// 推定誤差が時間刻み幅 h あたりの許容誤差よりも小さいかどうかチェック
if (errmax <= fabs(eb*h)) { // 小さい場合// 時刻を進める
t += h;// 進んだ分だけ残り時間は減らす (最終時刻 - 現在時刻)
h = tn - t;// 残り時間がなければ (計算しきったら) 戻る
if (fabs(h) <= epsmac)
return;// まだ残り時間があれば続行
x0 = xn;
12

t = 0.7 くらいから、刻み幅の調節が実際に動き始め、t = 1 の手前でぐんぐん小さくなっていく。コンパクトなスペースで紹介するにはどうしたら良いか分からないので結果はここには載せない

17}else {// 推定誤差が時間刻み幅 h あたりの許容誤差よりも大きい場合// 推定誤差と許容絶対誤差の大小を比較
if (errmax > epsv)// 推定誤差が許容絶対誤差より大きければ刻み幅を一気に 10% に縮める
h *= 0.1;
else {
h *= 0.9 * sqrt(sqrt(eb * fabs(h) / errmax));
if (fabs(h) < epsmax)
break;}}}}
fprintf(stderr,"rkf45(): t=%g で困難に出会いました。計算を中断します。", t);return 1;(なお、この手の問題のシミュレーションについては、伊理・藤野 [10] にも何か書いてあったような気がする。陳蘊剛さんの爆発のシミュレーションはどうだったかな…そのうちに一度まじめに調べてみよう。まあ、これは独り言みたいなものです。)

4.7.4

RKF45 自作プログラム開発

testrkf3.c/** testrkf3.c*/#include <stdio.h>#include <math.h>#include "rkf3.h"void myf(double t, double *x, double *f){f[0] =
x[1];f[1] = - x[0];}int main(){
int i, n;double a = 0, b = 1;
double x[2], newx[2], localerror, error, diff, lastdiff, lasterror;
double k[2][RKF45_STAGE], work[2];double h;printf("真の値=%20.15f\n", cos(b));
printf("nx
誤差 推定誤差\n");for (n = 1; n <= 10000; n *= 2) {x[0] = cos(a); x[1] = - sin(a);h = (b - a) / n;error = 0;
for (i = 0; i < n; i++) {18

bf_rkf45(2, myf,a + i * h, x, h, newx, &localerror,k, work);x[0] = newx[0];x[1] = newx[1];error += fabs(localerror);}diff = fabs(x[0] - cos(b));
printf("%5d %18.15f %7.1e %7.1e", n, x[0], diff, error);
if (n != 1) {printf(" %4.1f 分の 1 %4.1f 分の 1", lastdiff / diff, lasterror / error);}printf("\n");lastdiff = diff; lasterror = error;}
return 0;}

rkf3.h/** rkf3.h*/#define RKF45_STAGE (6)typedef void function(double t, double *x, double *f);
int bf_rkf45(int d, function f,
double t, double *x, double h,double *nextx, double *error,
double k[][RKF45_STAGE], double *work);

rkf.c
/*
* rkf3.c --- RKF45 (多次元, 何次元でも動くように)*/#include <stdio.h>#include <math.h>#include "rkf3.h"/* fabs() */static double rkf45_alpha[RKF45_STAGE] = {0.0, 1.0/4, 3.0/8, 12.0/13, 1.0, 1.0/2};
static double rkf45_beta[RKF45_STAGE][RKF45_STAGE-1] = {
{
0.0,
0.0,
0.0,
0.0},
{
1.0/4,
0.0,
0.0,
0.0},
{
3.0/32,
9.0/32,0.0,
0.0},{1932.0/2197, -7200.0/2197, 7296.0/2197,
0.0,
0.0},
{ 439.0/216,-8.0,
3680.0/513, -845.0/4104,
0.0},{ -8.0/27,
2.0,-3544.0/2565, 1859.0/4104, -11.0/40}};
static double rkf45_mu5[RKF45_STAGE] = {16.0/135, 0.0, 6656.0/12825, 28561.0/56430, -9.0/50, 2.0/55};
static double rkf45_mu4[RKF45_STAGE] = {25.0/216, 0.0, 1408.0/2565, 2197.0/4104,};-1.0/5,static double rkf45_mu_diff[RKF45_STAGE] = {

19

0.0

-1.0/360, 0.0, 128.0/4275, 2197.0/75240, -1.0/50, -2.0/55};
static double sqr(double x) { return x * x; }#include "rkf_utilities3.c" /* dotprod(), etc. *//** d 次元の ODE x’(t)=f(t,x) を解く。
* 時刻 t で値が x だとして、h だけ時間を進めたときの値 nextx を計算する。
*
* 注意:
*
x と next x は同じメモリー領域を指さないようにしておくこと。つまり*
bf_rkf45(d, f, t, x, h, x, err, k, w);*のようにして呼び出すと結果は壊れてしまう。
*--- こういう仕様にした方が全体に無駄が少なくなると考えた。
** work は問題の次元 d だけの長さの double 型の作業用領域*/
int bf_rkf45(int d, function f,
double t, double *x, double h,double *nextx, double *error,
double k[][RKF45_STAGE], double *work){
int i, j, s = RKF45_STAGE;/* x と nextx が重ならないことをチェック */if (x == nextx) {fprintf(stderr, "bf_rkf45(): x と nextx は同じメモリー領域ではいけない！ ");*error = 1e+10;
return -1;}/* k1,k2,k3,k4,k5,k6 を計算する */for (i = 0; i < s; i++) {double *xx = nextx; /* ちょっとの間使わせてもらう (お行儀が悪いけど) */double *kk = work;for (j = 0; j < d; j++)xx[j] = x[j] + dotprod(i, rkf45_beta[i], k[j]);f(t + rkf45_alpha[i] * h, xx, kk);for (j = 0; j < d; j++)k[j][i] = kk[j] * h;}/* 次のステップの値を 5 次公式で計算する＆ 4 次公式の誤差を推測する */for (j = 0; j < s; j++) {nextx[j] = x[j] + dotprod(s, rkf45_mu5, k[j]);work[j] = dotprod(s, rkf45_mu_diff, k[j]);}/* 推定誤差のノルム */*error = 0;
for (j = 0; j < d; j++) *error += sqr(work[j]);*error = sqrt(*error);/* もし 4 次の値が欲しければ */for (j = 0; j < d; j++)work[j] -= nextx[j];return 0;}

testrkf1 の実行結果
oyabun% ./testrkf3
真の値=
0.540302305868140

20

n
x
誤差 推定誤差
1 0.541185897435897 8.8e-04 1.4e-03
2 0.540325560014864 2.3e-05 8.2e-05 38.0 分の 1 16.8 分の 1
4 0.540302920658938 6.1e-07 5.0e-06 37.8 分の 1 16.2 分の 1
8 0.540302323044084 1.7e-08 3.1e-07 35.8 分の 1 16.1 分の 1
16 0.540302306371086 5.0e-10 2.0e-08 34.2 分の 1 16.0 分の 1
32 0.540302305883314 1.5e-11 1.2e-09 33.1 分の 1 16.0 分の 1
64 0.540302305868605 4.7e-13 7.6e-11 32.6 分の 1 16.0 分の 1
128 0.540302305868154 1.4e-14 4.8e-12 32.3 分の 1 16.0 分の 1
256 0.540302305868140 4.4e-16 3.0e-13 32.5 分の 1 16.0 分の 1
512 0.540302305868139 6.7e-16 1.9e-14 0.7 分の 1 16.0 分の 1
1024 0.540302305868140 4.4e-16 1.2e-15 1.5 分の 1 15.9 分の 1
2048 0.540302305868141 1.3e-15 7.7e-17 0.3 分の 1 15.2 分の 1
4096 0.540302305868142 2.3e-15 9.1e-18 0.6 分の 1 8.4 分の 1
8192 0.540302305868144 4.1e-15 5.7e-18 0.6 分の 1 1.6 分の 1
oyabun%testrkf1 の実行結果の分析 区間 [0, 1] を 1 等分した場合 (h = 1) でもそれなりに解けているのがすごい。h が極端に大きいときと極端に小さいときを除くと、h を半分にするごとに誤差は 32 分の1 程度になっていて、公式が 5 次の精度であることを裏付けている。一方、推定される誤差は 4 次
公式と 5 次公式の差を使っているので、4 次公式の誤差を推定したものなので、h が半分になるごとに 16 分の 1 程度になっている。

典型的なスキーム (2) 線形多段法

5(この節は工事中というか、大工事が必要と考えて欲しい。)
k 段法の公式のところで、Φ が線形の場合、すなわち
d0 xj + d1 xj+1 + · · · + dk xj+k = h (β0 fj + β1 fj+1 + · · · + βk fj+k )を線形多段法 (linear multistep method) と言う。ただし fj := f (tj , xj ).

(explicit(陽的)def.

⇔ βk = 0

def.

implicit(陰的) ⇔ βk ̸= 0
def.

ρ(λ) = d0 + d1 λ + · · · + dk λk の根 λ1 , · · · , λk について
|λi | ≤ 1 (i = 1, · · · , k),|λi | = 1 となる λi は単根。

という条件が安定条件である。
多段法には以下の特徴がある。

(1) 出発にあたって、未知の値 x1 , · · · , xk−1 を何らの方法で求めねばならない。
(2) 計算の途中で stepsize h を変更するのが難しい。
(3) ステップあたりの f の計算回数が少ないままで、高次の公式が作れる。
(4) 安定性に注意が必要。

5.1

PC (予測子修正子法)

(準備中)
21

6

その他の方法

6.1

Taylor 法

f の Taylor 展開を利用する方法である。簡単で

d
dx

k
f の計算しやすい f に対して有効である。

例 6.1 (Bessel 関数) これは一松 [7] に載っていた話。Bessel 関数 Jn は、Bessel の微分方程式

x2 y ′′ + xy ′ + (x2 − n2 )y = 0(7)を満たす。この方程式を一回微分すると

x2 y ′′′ + 3xy ′′ + (x2 − n2 + 1)y ′ + 2xy = 0.
k ≥ 2 に対して、(7) を (Leibniz の公式を用いて) k 回微分すると
x2 y (k+2) +k·2x·y (k+1) +k(k − 1)k(k − 1)·2·y (k) +xy (k+1) +k·1·y (k) +(x2 −n2 )y (k) +k·2x·y (k−1) +·2·y (k−2) = 0.
2
2

整理して

x2 y (k+2) + (1 + 2k)xy (k+1) + (x2 − (n2 − k 2 ))y (k) + 2kxy (k−1) + k(k − 1)y (k−2) = 0 (k ≥ 2)を得る。(ちなみにこの式は、k = 1 のとき k(k−1)y (k−2) = 0, k = 0 のとき 2kxy (k−1) +k(k−1)y (k−2) =(k)
0 とみなせば、任意の k ≥ 0 について成立する。) それゆえ ak := Jn (x0 ) とおくと、



n2
1

 a2 + a1 + 1 − 2 a0 = 0,


x0
x0






n2 − 1
2
3
a1 + a0 = 0,
a3 + a2 + 1 −
2

x0
x0
x0







2k + 1n2 − k 2
k(k + 1)2k
 ak+2 +ak+1 + 1 −ak + ak−1 +ak−2 = 0
2
x0
x0
x0
x20
という漸化式を得る13 。ゆえに a0 = Jn (x0 ), a1 = Jn ′ (x0 ) が分かれば、順に a2 , a3 , · · · が計算できる。ゆえに x0 の回りの望む階数までの Taylor 展開

Jn (x) ≒

m
X
ak

k!
k=0

(x − x0 )kが計算できる。適当な x1 に対して、

Jn (x1 ) ≒

m
X
ak
k=0

k!

(x1 − x0 ) ,k′

Jn (x1 ) ≒mXak(x1 − x0 )k−1
(k − 1)!
k=1(k)を計算して、再び上と同様にして {Jn (x1 )}k=0,1,2,··· が計算できるので、x1 を中心とする Jn の Taylor
展開が計算できる。以下これを繰り返すことで、広い範囲に渡って Jn の値が計算できる。

6.2

補外法

補外法 (extrapolation method) と呼ばれる一群の方法がある。特に、Richardson 補外に基づ
く Bulirsch-Stoer 法が有力であるとか。“Numerical Recipes” に Bulirsch-Stoer 法のプログラムが載っている。
(Stoer-Bulirsch [11] に載っているのがそう？ あまり詳しくないようだが。)
13

一松先生の本には、最後の項を x20 で割るのを忘れるという、軽い誤植があります。

22

数値的安定性

7

今まで区間 [a, b] を固定して、分割数 N → +∞ とした (h = (b − a)/N ) ときの収束を考えた。応用上は長時間解を追跡したい場合がある。この場合、刻み幅 h を固定して n → +∞ (t → ∞) としても「変なこと」が起こらないようにしたい。
注意 7.1 (「安定」に関する注意) 常微分方程式の数値解法の話では「安定」という語が何度も出て来たが、状況により違う意味で使われることが多い。まとめておくと

1. 微分方程式の (平衡点, あるいは周期解の) 安定性
2. 数値解法の安定性
3. 数値的安定性(離散化とは無関係)(離散化に起因する不安定性が起こらない; 丸め誤差とは無関係)(丸め誤差の増幅に起因する)

ただし、それぞれの場合にさらに細かい分類がある。
簡単なテスト問題に適用したときの数値的安定性を調べる、というのが基本的な作業方針 (テスト
問題の枠を離れたときにどうなるかの保証はないわけで、そういう意味ではあまりいばれないが)。

7.1

線形安定性解析λ ∈ C を Reλ < 0 なる定数として(
def.
x′ (t) = λx(t) (t ∈ I = (a, +∞))
x(a) = x0
について適用してみる。この解は x(t) = x0 eλt であり、|x(t)| = |x0 | eRe λt であるから、 lim x(t) = 0t→∞である。ところが、数値解法ではしばしば∃h1 > 0, ∃h0 ∈ (0, h1 ) ∀h > h1s.t.

lim |xj | = +∞j→∞

 0 < ∀h < h0

lim |xj | = 0j→∞ということが起こる。

1) 中点則の場合 (leapfrog method とも言うらしい). これはいわゆる中心差分近似
x′ (t) =を用いて作った公式

x(t + h) − x(t − h)+ O(h2 )2h(xj+1 = xj−1 + 2hfj ,x1 は適当に定めるのことを指す。これは任意の h > 0 に対して limj→∞ |xj | = +∞ を満たす。
中心差分近似に基づく中点公式 xn+1 = xn−1 +hf (tn , xn ) で x′ = λx, x(0) = 1 を解く (λ = −4.35)。
左側は x1 = eλh としたもの、右側は x1 = x0 + hf (t0 , x0 ) としたもの。

23

Midpoint rule for dx/dt=-4.35 x, x(0)=1. (h=0.001)
3

Midpoint rule for dx/dt=-4.35 x, x(0)=1. (h=0.001, appro. X1)
2000

2

1500

1000

1

500

0

0

-1
-500

-2

-3
0

-1000

1

2

3

4

5

-1500
0

1

2

3

4

5

2) R-K 型公式の場合. 公式の段数を s, 次数を m とすると、
xj+1 = R(λh)xjのように書ける。ここで R(z) はR(z) =s 次以下の多項式,s 次以下の多項式|R(z) − ez | = O(|z|m ) (|z| → 0)

なる z の有理式である (特に公式が explicit の場合には、R(z) は z の多項式になる)。例えば
z2 z3 z4
Euler 法の場合 R(z) = 1 + z, 古典的 Runge-Kutta 法の場合、R(z) = 1 + z + + + . この2
3! 4!
R(z) に対して
def.
R = {z ∈ C; |R(z)| < 1}
を絶対安定領域といい、λh ∈ R =⇒ lim xj = 0
j→+∞がなりたつ。絶対安定領域が左半平面 {z ∈ C; Re z < 0} を含むような公式を A-安定 (A-stable)という。A-安定な公式では、Re λ < 0 なるとき、任意の h > 0 に対して数値解が 0 に収束する。
次に示すのは Runge-Kutta 法の絶対安定領域である。カラーでお見せできないのが残念ですが。

8
6
4

2

2
0
-3

0
-2
-1
-2

0
1

なお、レベル 1 の等高線が負軸と交わる点の座標は (巾根で書けはするのだけれど...)、−2.7852935634052816235297591900854630347647578967169 · · · .
なお、後述の硬い方程式の項を参照せよ。

24

3) LM 法の場合. 公式(α0 − zβ0 )xj + (α0 − zβ1 )xj+1 + · · · (α0 − zβk )xj+k = 0
に対して特性方程式を(α0 − zβ0 ) + (α0 − zβ1 )ξ + · · · (α0 − zβk )ξ k = 0
で定義し、その根を ξℓ (z) (ℓ = 1, · · · , k) とする。この場合
def.

R = {z ∈ C; |ξℓ (z)| < 1}
とおくとλh ∈ R =⇒ lim xj = 0.
j→∞

Stiﬀ problem (硬い問題)

8

以下は工学的見地からの説明である (新しい言葉に説明がついているが、数学的な定義になっていないものが多い)。
じていすう
1
時定数 (time scale) とは、解が に減衰するのに必要な時間のことである14 。 x(t) が τ を正の定
e
数として


t
x(t) = exp −τのように表されているのならば τ が時定数である。


定義 8.1 (硬い方程式のあいまい定義 — その 1) ある安定な微分方程式が解くべき全区間に比較して極めて小さい時定数をもつ、指数関数的に減衰する解を一つの特解として持つ時、その方
程式は stiﬀ である (硬い方程式である) と呼ばれる。






定義 8.2 (硬い方程式のあいまい定義 — その 2) 一つの問題の中に時定数が大きな部分と、小さな部分がある時、その方程式は stiﬀ である (硬い方程式である) と呼ばれる。


こういう問題を数値解法で解く場合、数値的安定性の要請から小さな時定数にあわせて h を選ぶと、なかなか計算が進まない。これは数値解法に特有の困難である (もとの問題自身は安定である)。
例 8.3 λ1 ≪ λ2 < 0 なる定数 λ1 , λ2 に対して、方程式
!
!
d
λ1 0
x1
=
dt x2
0 λ2
は硬い。

A-stable な方法があればよいが、14

半減期というものに似ている。

25





定理 8.4 (Germund Dahlquist)(i) 陽的 Runge-Kutta 法は A-stable になり得ない。

(ii) 陽的線型多段法は A-stable になり得ない。
(iii) 陰的線型多段法も 3 次以上の公式は存在しない。
(iv) 2 次 A-stable 陰的線型多段法の中では台形則が「最適」である。


のような「否定的な」結果がある。そこで対応策として、• A-stable はあきらめ、別の (もう少し弱くした) 安定性の概念を考える。Stiﬀ stability (S-安定ギ ア性) の概念。Gear の方法。

• A-stable な陰的 Runge-Kutta 法はかなり次数の高い公式が作られている。
などが考えられる。

参考書

9[1] 笠原 [12] は、明治大学数学科の常微分方程式論の講義で、長くテキストとして採用されていた。
数学科の教科書として標準的な内容である。
[2] ハイラー・ネルセット・ヴァンナー [13], ハイラー・ヴァンナー [14] は、常微分方程式の初期値問題の数値解法に関する定評のある文献の、まっとうな翻訳である。著者の Hairer のサイト
http://www.unige.ch/~hairer/software.html にプログラムが掲載されている。
この本とこのサイトの使い方を説明する方が生産的なのかもしれない (2019/3/6 久しぶりに加筆)。
[3] 三井 [15] は、常微分方程式の数値解析の専門家が書いた数少ない和書であり、特に理論的な解説をきちんとしてある本として和書ではユニークであった。これをコンパクトにしたものに同じ著者の三井 [2] がある。これは三井 [16] の単行本化である。コンパクトではあるが書かれたのが新しいだけの利点があると思う。
[4] 戸川 [17] には、色々な Runge-Kutta 型公式 (made in Japan も多い！ ) のプログラムがたくさん載っている。
[5] 渡部・名取・小国 [18] は、題名からするとプログラムのみの本のようであるが、そうではなく、常微分方程式の項は杉原正顕氏によるすぐれた解説がある。
[6] William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery, Numerical
Recipes (in C), Cambridge University Press, 邦訳 技術評論社.
Bulirsch–Stoer 法が最善の方法であると信じている著者たちのプログラムが載っている (日本語で読める本としては珍しいと思われるので貴重15 )。
[7] Geometric Numerical Integration は何と訳すのかな。Hairer-Lubich-Wanner [19], 初版の書評
McLachlan [20] (2019/3/6 加筆)
15

この本は非常に有名だが、叩く人も多くて「悪評」だけ読んでいると敬遠してしまいそう。しかし参考になるところは多い。実際この本よりも役立たずの本は世の中にたくさんある。つまり「この本一冊で OK だ」という捉え方をする
人 (欧米では多いらしい) に警告を発するために批判されているらしい。

26

Hairer-Hairer [21] も見てみよう。
やはり http://www.unige.ch/~hairer/software.html に Fortran, MATLAB のプログラムがある。

おまけ — 実際的な誤差の推測

10

急速に α に収束する列 {xn }n∈N があるとき、
def.

εn = ∥xn − α∥で定義される誤差の大きさについて、十分先の番号 n に対してはεn+1 ≪ εnが成り立つので、∥xn − xn+1 ∥ ≤ ∥xn − α∥ + ∥xn+1 − α∥ = εn + εn+1 ≒ εn .
よって

∥xn − xn+1 ∥を xn の誤差の大きさ εn の見積りとすることが出来る。
急速に収束しない列の場合はどうか？ 例えば n を分割数とした時の差分法の解 u(n) などでは、この仮定が成り立たないと思われる。そういう場合は例えば
def.

jxj = u(2 )とすることによって、同じテクニックが使える。つまり例えば n = 512 のときの近似解と n = 1024の時の近似解の差の大きさを、n = 512 の時の近似解の誤差の大きさの見積りとすることが出来る。

A

数値計算するための情報

A.1

はじめに

常微分方程式の初期値問題をコンピューターで解くとき、次の二点が問題になることが多い、と思う。

1. 解の可視化が必要になる可能性が高いが、それをどう実現するか。
2. 問題の次元が高い場合、ベクトルを扱える言語を使うことが望ましいが、何を選択するか。
問題の次元が 1 や 2 であれば、言語や処理系は何でも良いような気がする。それこそ「仮称 (十進)BASIC」16 などを使うと、とりあえずグラフィックスが使えるので気軽に試せて良いと思われる。
ずっと以前の相場では、FORTRAN 77 くらいのレベルの低いプログラミング言語を使ってプログ
ラムを苦労して書いて、可視化は XY プロッターや、それと同程度の機能を提供するグラフィック
ス・ライブラリィを利用する、ということだったと思う (例えば、森 [22] などを見よ)。やり方がま
ずいと、対話性の低いプログラムが出来やすいと思われる。(脱線かも知れないが、ふと、グリック[23] の中のエピソードが浮かんできた。)
今ではもう少し上手いことが出来ると思われるが、案外、それと同程度のことしかやっていない場
合が多いかもしれない。例えば、最近の日本の大学の風潮として、C 言語でプログラムを書くという場面が多そうだが、それは FORTRAN 時代とあまり変らない気がする。
実は、身の回りで常微分方程式の初期値問題を真剣に解く状況に遭遇したことがないので、良く分からないが、「気軽な試み」はいくつかしている。
16

http://hp.vector.co.jp/authors/VA008683/

27

A.2

C 言語+グラフィックス・ライブラリィ

実に色々ある。以前だったら、全部実際に試して比較してみているところであるが、今はそういう暇がない。

• GLSC17
2 変数関数のグラフの鳥瞰図が気軽に描ける点が「受けている」理由なのかと思われる。
「GLSC の紹介」18• PGPLOT19
正直に白状すると、本当の意味で使った経験がない。世の中の評価は高いようなので、これが手堅いのか？ という気がしているのだけど。
• EGGX/ProCALL20
X Window System 用。ある時期使ってみたことがある。気に入っている点も多い。
• GrWin21
Windows 用である。マウスなどが使える。
• nxgraph.h22 — 「計算物理のための C/C++言語入門」23

A.3

C 言語+gnuplot

UNIX 風の環境では、gnuplot のインストールが簡単で、C からパイプを使ったプロセス間通信でgnuplot をコントロールすることが出来る。gnuplot ですべて済むとも思えないが、それで済む場合
は、それを使うのが楽であるし、便利であると思われる。
「C から gnuplot を呼び出す」24
なお、以前の gnuplot は「簡単なグラフを簡便に書くためのツール」であったが、今では gnuplot
自身で Runge-Kutta 法などの計算を行なわせることが可能である。例えば、「gnuplot で微分方程式25を解く」 が参考になる。

A.4

Java

あまり自信はないが、ほどほどの効率で、対話性の高いプログラムを書くためには Java が便利ではないか？ と思っている。いわゆる GUI が比較的簡単に実現できる。研究室で作成したプログラム
例をあげておく。
!
!
!
d x
a b
x
=
のシミュレーション。http://nalab.mind.meiji.ac.jp/~mk/labo/
•
dt y
c d
y
java/prog/ODE1.html

17

ftp://ftp.st.ryukoku.ac.jp/pub/ryukoku/software/math/
http://nalab.mind.meiji.ac.jp/~mk/labo/howto/intro-glsc/intro-glsc.html
19
http://www.astro.caltech.edu/~tjp/pgplot/20
http://www.ir.isas.jaxa.jp/~cyamauch/eggx_procall/21
http://spdg1.sci.shizuoka.ac.jp/grwinlib/
22
http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/nxgraph.html
23
http://cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/
24
http://nalab.mind.meiji.ac.jp/~mk/labo/howto/intro-gnuplot/node21.html
25
http://www.ss.scphys.kyoto-u.ac.jp/person/yonezawa/contents/program/gnuplot/diff_eq.html
1828

• 渦糸のシミュレーション。http://nalab.mind.meiji.ac.jp/~mk/labo/java/uzu.html(2003 年頃)
• オレゴネーター http://nalab.mind.meiji.ac.jp/~mk/labo/report/open/2002-naitou-prog/
Oregonator.html(これも 2003 年頃)
残念ながらここで紹介することは出来ないのだが、実は他所の研究室で、非常に対話性の高い、見事なプログラムを見たことがある。

(2019 年筆 Java について、個人的に 2000 年頃感じていた期待はいつの間にかしぼんでしまった。
一つの理由として、Internet 上で、主にセキュリティの観点から、Java が利用しにくい状況になったことがある。もう一つは、たとえ Java アプレットで気軽にシミュレーションが出来ても、それだけでは「イマイチぱっとしない見もの」に過ぎず、意外と教育 (学習) 効果が上がらないらしい、と分かったことがあげられる。)

A.5

C++ の利点

C++ は C と大差ないように思えるかも知れないが、ベクトルを利用しやすくなっているので、案外大きな差があるのではないか？ と想像している。この点で、牧野 [24] に興味を持っている。
グラフィックスについても、何か面白い工夫がありうると思っているのだが…
Eigen というのを知った (2013/12, http://nalab.mind.meiji.ac.jp/~mk/knowhow-2013/node32.
html)。これを使うと良いかも。
久しぶりに卒業研究で常微分方程式のシミュレーションをする学生が現れた (複数)。C++ & Eigen& GLSC というラインアップでプログラムを作ってもらった (2019/3/6 加筆)。内容によっては GLSC3Dを使うと良いのかもしれない (今回は投影図で済ませた)。

A.6

Ruby

(これは 2011 年頃に書いた。)
前項に続き、牧野センセがらみになるが、Ruby を使うのも面白い。
牧野 [25] では、Ruby+NArray26 +PGPLOT27 +Ruby/PGPLOT28 +rongo29 というライン・アップで常微分方程式のシミュレーション＆可視化をしている。
そのココロは、• 簡単な常微分方程式ならば、効率は大した問題ではないから、インタープリター である Ruby
くらいでも十分だ。
• むしろベクトルなどが気軽に扱えるのは便利だ。NArray を使えば効率も十分高いし。
• rongo のように言語を拡張してグラフィックス機能を持たせたように出来るのは、嬉しい。
というところだろうか (私のいいかげんなまとめが信用できない人は、上の本を読んで下さい) 。

• Rubyist Magazine「数値計算と可視化」30
26

http://narray.rubyforge.org/index.html.ja
http://www.astro.caltech.edu/~tjp/pgplot/28
http://raa.ruby-lang.org/list.rhtml?name=ruby-pgplot
29
http://jun.artcompsci.org/softwares/rongo/30
http://jp.rubyist.net/magazine/?0006-RLR
27

29

• Ruby for Scientiﬁt Computing31
Ruby のような言語から、各種のライブラリィを簡単に呼び出すというのは面白い。
Scipy32 なんてのもある。

B

眠りから覚めてみたら (2022/6)

現在、現象数理学科というところに所属しているが、学部のゼミの様子が数学科の時と違っていて、試行錯誤 (悪戦苦闘) している。色々あって、このところゼミでは最初に常微分方程式の初期値問題
をやっている。という訳で、長いこと放っておいたものをひっぱり出す羽目になっている。
今さら C 言語でプログラムを書く気はなくて、かと言って Java という気にもなれず、Python, Julia,
Ruby (or Crystal) のどれがいいかなあ、と迷っているが、学生は C 言語よりは Python という人が多いので、何となく Python に流れている。

参考文献[1] 桂田祐史：常微分方程式の初期値問題の数値解法入門, http://nalab.mind.meiji.ac.jp/~mk/labo/text/num-ode.pdf (1995〜2011, 2021 小改訂).
[2] 三井斌友：常微分方程式の数値解法, 岩波書店 (2003), 「微分方程式の数値解法 I」岩波講座応
用数学 (1993) の単行本化.
こ と う としゆき

よしひろ

[3] 三井斌友, 小藤俊幸, 齊藤善弘：微分方程式による計算科学入門, 共立出版 (2004).
[4] E. ハイラー, S. P. ネルセット, G. ヴァンナー：常微分方程式の数値解法 I 基礎編, シュプリン
ガー・ジャパン (2007).
[5] E. ハイラー, G. ヴァンナー：常微分方程式の数値解法 II 発展編, シュプリンガー・ジャパン (2008).
[6] 桂田祐史, 佐藤篤之：力のつく微分積分 — 1 変数の微積分, 共立出版 (2007).
ひとつまつしん[7] 一 松 信：数値解析, 朝倉書店 (1982/10).
[8] Fehlberg, E.:
Klassische Runge-Kutta-Formeln vierter und niedrigerer Ordnung mit
Schrittweiten-Kontrolle und ihre Anwendung auf Wärmeleitungsprobleme, Computing, Vol. 6,
pp. 61–71 (1970).
[9] 森正武：FORTRAN 77 数値計算プログラミング, 岩波書店 (1986, 1987).
よりたけ[10] 伊理正夫, 藤野和建：数値計算の常識, 共立出版 (1985).
[11] Stoer, J. and Bulirsch, R.: Introduction to numerical analysis, second edition, Springer-Verlag(1992), (R. Bartels, W. Gautschi, and C. Witzgall による翻訳である。).
こうじ[12] 笠原晧司：微分方程式の基礎, 数理科学ライブラリー, 朝倉書店 (1982).
31
32

http://pub.cozmixng.org/~the-rwiki/rw-cgi.rb?cmd=view;name=Ruby+for+Science
http://www.scipy.org/

30

[13] E. ハイラー, S. P. ネルセット, G. ヴァンナー：常微分方程式の数値解法 I 基礎編, シュプリン
ガー・ジャパン (2007), 三井斌友監訳.
[14] E. ハイラー, G. ヴァンナー：常微分方程式の数値解法 II 発展編, シュプリンガー・ジャパン (2008),
三井斌友監訳.
たけとも[15] 三井斌友：数値解析入門, 朝倉書店 (1985).
[16] 三井斌友：微分方程式の数値解法 I, 岩波講座 応用数学, 岩波書店 (1993).
[17] 戸川隼人：UNIX ワークステーションによる科学技術計算ハンドブック 基礎篇 C 言語版, サイエ
ンス社 (1992), http://www.mscom.or.jp/~tog/anna/ でプログラムが公開されていたけれど、
今はない？
[18] 渡部力, 名取亮, 小国 力監修：Fortran 77 による数値計算ソフトウェア, 丸善 (1989 (1982?), 常微分方程式の項は杉原正顕氏によるすぐれた解説がある。
[19] E. Hairer, C. L. and Wanner, G.: Geometric Numerical Integration: Structure-Preserving
Algorithms for Ordinary Diﬀerential equations, Springer Series in Computational Mathematics, Springer (2010/3/11), ISBN: 978-3642051579, http://www.unige.ch/~hairer/software.
html.
[20] McLachlan, R.: Review: Featured Review: Geometric Numerical Integration: StructurePreserving Al gorithms for Ordinary Diﬀerential Equations. By E. Hairer, C. Lubich, and G.
Wanner. Springer-Verlag, New York, 2002. $84.95. xiii+515 pp., hardcover. ISBN 3-540-43003-2.,
SIAM Review, Vol. 45, No. 4, pp. 817–821 (2003).
[21] Hairer, E. and Hairer, M.: GniCodes — Matlab Programs for Geometric Numerical Integration(2003), http://www.unige.ch/~hairer/preprints/gnicodes.ps.gz.
[22] 森正武：FORTRAN 77 数値計算プログラミング, 岩波書店 (1987).
[23] ジェイムズ・グリック：カオス — 新しい科学をつくる, 新潮文庫, 新潮社 (1991), James Gleick,
Chaos — making a new science, Viking Penguin (1987) の上田 睆亮監修, 大貫昌子による訳であるが、 せっかくの翻訳が入手しにくくなっている。中古をいとわなければ入手できる。
[24] 牧野淳一郎：パソコン物理実地指導, 共立出版 (1999).
[25] 牧野淳一郎：とんでる力学, 丸善 (2005).

31