185

第 16 章 常微分方程式の初期値問題 — 一段法
山口 ・・・(略)・・・計算学そのものはものすごく健康な応用数学の一つの流れだと思うんだけれども，今までに行われてきている数値解析というのは，その計算学がそのまま数学の影響を受けて，純粋な数値解析というような方向へ向かってしまって，一つの行きづまりを迎えている。例を挙げればルンゲクッタ法ですね。
ルンゲが考えたときには猛烈に健康だったと思うし・
・
・
（略）
・
・
・それがだんだん発展していくと，たとえばルン
ゲクッタ法のいまの研究というのは，いかなる常微分方程式にもこれだけの次数のルンゲクッタ公式を作ったら，これだけの精度が出る，と。現象という側面はもう欠け落ちてしまっているのです。
山口昌哉

編「数値解析と非線型現象」(日本評論社)

16.1 常微分方程式
常微分方程式とは導関数を含む方程式で，一般には以下のように表わすことが出来る。
!
dy d2 ydr y
ϕ x, y, , 2 , ..., r = 0 (x ∈ R, y ∈ Rn )(16.1)dx dx
dx
この n 次元 r 階常微分方程式を満足する一変数関数 y(x) を見つけ出すには，代数的な操作を行うことになるが，積分と同様，全ての常微分方程式に対して適用出来る公式というものは存在しない。
従って，一般には近似解を求めることしか出来ない。
本章では特に最高階数の導関数について陽的 (explicit) な

dr y edy d2 ydr−1 y= ϕ x, y, , 2 , ..., r−1
r
dx
dx dx
dx

!
(16.2)という形式で表現出来る常微分方程式のみ取り扱うことにする。
基本となるのは一階の一次元常微分方程式dy= f (x, y)dxである。2 階以上の常微分方程式については


 


 y1  y 
 y2  

dy/dx

 = 

..

 ..  .
 .    r−1r−1 d y/dxyr(16.3)

第 16 章 常微分方程式の初期値問題 — 一段法

186と置き換えることによって，r 階常微分方程式 (16.2) を 1 階の常微分方程式

  y1  

y2

 
 y2  

y3

 d  .  .
..
 ..  =  
dx  y 
yr
 r−1  


  e
yr
ϕ(x, y1 , y2 , ..., yr−1 )(16.4)と同一視出来る。従って，以降は n 次元 1 階常微分方程式dy= f(x, y)dx(16.5)

のみ考えることにする。

16.2 初期値問題と Lipschitz 条件
常微分方程式の解は無数に存在する。最も簡単な次の例題でそれを示す。
例題 16.2.1

1 次元の常微分方程式dy=ydxの解は右辺を左辺に移項して両辺を x について積分することによって得られる。この時，積分定数c ∈ R が残り，結局，解 y(x) はy(x) = c exp(x)となる。即ち，この解は無数に存在する。
解を一意に定めるためには，常微分方程式に対して条件を設定する必要がある。そのうち，変数 x が x = x0 と固定された地点での y(x0 ) = y0 を与えられた問題を，「初期値問題」(initial value

problem, IVP) と呼び，この値を初期値と呼ぶ。従って，常微分方程式の初期値問題とはdy= f(x, y)dx y(x0 ) = y0

(16.6)という形で与えられる。
この初期値問題の解が一意に決まるためには次の Lipschitz 条件を必要とする。
定理 16.2.2 (Lipschitz 条件)

x ∈ [x0 , α] と任意のベクトル z1 , z2 ∈ Rn に対し∥f(x, z1 ) − f(x, z2 )∥ ≤ L∥z1 − z2 ∥(16.7)

を満足する正の定数 L ∈ R が存在する時，常微分方程式の初期値問題 (16.6) は一意の解 y(x) を持つ。
先の例題でこれを確認してみる。

16.3. 差分からの導出: Euler 法, 中点法, 古典的 Runge-Kutta 法

187

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

 dy


 dx


 y(0)==y1

は L = 1 で Lipschitz 条件 (16.7) を満足する。この解はy(x) = exp(x)となり，一意に定まる。
本章では，Lipschitz 条件を満足する常微分方程式の初期値問題 (16.6) 用の数値解法について簡単に触れる。

16.3 差分からの導出: Euler 法, 中点法, 古典的 Runge-Kutta 法1 階常微分方程式の初期値問題を数値的に解くには，初期値を含むある閉区間 [x0 , α] を l 分割
P
し，各地点 xi = x0 + i−1
j=0 h j における解 y(xi ) の近似解 yi を逐次求める。この分割した小区間の幅
hi (i = 0, 1, ..., l − 1) を刻み幅 (Step Size) と呼ぶ。刻み幅は特段支障が無い限りは全て等間隔にしても良いが，解の変動が激しい箇所や，「硬い (Stiﬀ)」問題では部分的に小さくしたり，全体的に小さく取る必要が出てくる。
最も単純なアルゴリズムは，前章で示したように一階導関数を差分商に置き換えることですぐに得られる。例えば，前進差分商 ∆y/h=(y(x + h) − y(x))/h で dy/dx を置き換えれば

yi+1 := yi + hi f(xi , yi )(16.8)という漸化式がすぐに導出出来る。これを Euler 法 (前進 Euler 法，陽的 Euler 法) と呼ぶ。
また，同様に中心差分商 δy/(2h)=(y(x + h) − y(x − h))/(2h) で置き換えれば

yi+1 := yi−1 + 2hi f(xi , yi )(16.9)となる。これを中点法と呼ぶ。中点法の場合は，最初の y−1 = y(x0 − h−1 ) を初期値とは別に与える
か，y1 を Euler 法で求めた後，y2 から出発する必要がある。通常は後者を使用する。
これらをまとめると次のようなアルゴリズムになる。

アルゴリズム 35 (Euler 法)

1. 各離散点 xi (i = 0, 1, ..., l − 1) に対して以下を計算する。

yi+1 := yi + hi f(xi , yi )アルゴリズム 36 (中点法)

1. y1 を Euler 法で求めておく。
2. 各離散点 xi (i = 1, ..., l − 1) に対して
以下を計算する。
yi+1 := yi−1 + 2hi f(xi , yi )

第 16 章 常微分方程式の初期値問題 — 一段法

188

この Euler 法，中点法の打ち切り誤差 (離散化誤差) は先に示した通りそれぞれ刻み幅の 1 乗，2乗に比例して減少する。更に打ち切り誤差を減少させる方法として次の古典的 Runge-Kutta 法がある。

アルゴリズム 37 (古典的 Runge-Kutta 法)

1. 各離散点 xi (i = 0, 1, ..., l − 1) に対して以下を計算する。
(a) k1 , k2 , k3 , k4 を次の式から求める。


k1







 k2



k3





 k
4

= f(xi , yi )= f(xi + 12 hi , yi + 12 hi k1 )= f(xi + 12 hi , yi + 12 hi k2 )= f(xi + hi , yi + hi k3 )(b) yi+1 := yi + 16 hi (k1 + 2k2 + 2k3 + k4 )これは中点法の倍，即ち刻み幅の 4 乗に比例して打ち切り誤差を小さくできる。
例題 16.3.1

1 次元常微分方程式

 dy


 dx


 y(0)==

x+y
1

を Euler 法，中点法，古典的 Runge-Kutta 法で解いた数値解の相対誤差を図 16.1 に示す。

相対誤差
1

0.001

1e-06

Euler, h=1/10

Midpoint h=1/10

Euler, h=1/104

ERK44, h=1/10
MIdpoint, h=1/104

1e-09

1e-12

ERK44, h=1/104

1e-15
0.2

0.4

x

0.6

0.8

1

図 16.1: Euler 法, 中点法, 古典的 Runge-Kutta 法の相対誤差

16.4. 一般の Runge-Kutta 法

189

問題 16.3.1

1. 上の例題 16.3.1 の常微分方程式の初期値問題の解析解を求めよ。また， x = 1 における近似値を h = 1/2,1/4 として Euler 法, 中点法でそれぞれ求め，相対誤差も併せて求めよ。
2. 古典的 Runge-Kutta 法が刻み幅の 4 乗に比例することを確認せよ。[Hint: f に対して 2 変数の Taylor 展開を適用し，4 次までその係数が一致することを確認すればよい。]

16.4 一般の Runge-Kutta 法
前述のように，常微分方程式には代数的演算では求積不可能なものが存在する。従って，解の全体像を把握するためには離散化による数値解法が不可欠である。以下では常微分方程式の初期値問
題における数値解法，特に一段法に分類される Runge-Kutta 法全般について，数値例を交えつつ述べることにする。

m 段 Runge-Kutta 法のアルゴリズムは一般にアルゴリズム 38 のように表わすことができる。

アルゴリズム 38 (m 段 Runge-Kutta 法)

1. 区間 [x0 , α] を l 分割し，各離散点を x0 , x1 ,..., xl = αとする。離散点 xi と xi+1 の刻み幅を hi とおく。
2. 各離散点 xi (i = 0, 1, ..., l − 1) に対して以下を計算する。
(a) k1 , k2 , ..., km を次の式から求める。
P
k1 = f(xi + c1 hi , yi + hi · mj=1 a1 j k j )




P

m



 k2 = f(xi + c2 hi , yi + hi · j=1 a2 j k j )

..



.





 km = f(xi + cm hi , yi + hi · Pm am j k j )j=1(b) yi+1 := yi + hi ·Pm

j=1 w j k j

定数 c1 , c2 , · · · , cm , a11 , ..., am,m , w1 , w2 , · · · , wm を表 16.1 のように表形式1 で表わす。

1 ステップの局所離散化誤差が O(h p+1 ) 次のとき，この公式は m 段 p 次であると呼ぶ。
Runge-Kutta 法は係数 (表 16.1) の決め方によって大きく 3 つに分類できる。
1. 陽的 Runge-Kutta 法
c2
c3
..
.
cm

a21
a31
..
.
am1

a32
..
.
am2

···

am,m−1

w1

w2

···

wm−1

..

.

1 この係数表を創始者の名を取って，Butcher Table と呼ぶこともある。

wm

第 16 章 常微分方程式の初期値問題 — 一段法

190

表 16.1: m 段 Runge-Kutta 法の係数

c1

a11

a12

···

a1m

c2
..
.
cm

a21
..
.
am1

a21
..
.
am2

···
···

a2m
..
.
am,m

w1

w2

···

wm

c1

=

0

ai j

=

0 (i ≤ j)• k1 から km まで逐次的に計算できる。
2. 半陰的 Runge-Kutta 法
c1
c2
..
.

a11
a21
..
.

a22
..
.

..

cm

am1

am2

···

am,m

w1

w2

···

wm

ai j

=

.

0 (i < j)• k j ( j = 1, ..., m) を求めるには，次のような n 元の非線型方程式を解かねばならない。
k j = f(xi + c j hi , yi + hi

j−1Xa js k j + hi a j j k j )(16.10)s=1

3. 陰的 Runge-Kutta 法
c1
c2
..
.
cm

a11
a21
..
.
am1

a12
a22
..
.
am2

···
···
···

a1m
a2m
..
.
am,m

w1

w2

···

wm

• k1 , ..., km を求めるには，次のような mn 次元非線型連立方程式を解かねばならない。
P
k1 = f(xi + c1 hi , yi + hi · mj=1 a1 j k j )




Pm




 k2 = f(xi + c2 hi , yi + hi · j=1 a2 j k j )(16.11)

..



.





 km = f(xi + cm hi , yi + hi · Pm am j k j )j=1

16.5. 陰的 Runge-Kutta 法の内部反復

191

Runge-Kutta 法の公式は Butcher[2]，Lawson, Huta, Shanks, 田中 [35, 36, 37, 38] らが精力的に開発し，その性能評価を行っている。
陰的 Runge-Kutta 法は，陽的 Runge-Kutta 法に比較して次の点で優れている。

1. 少ない段数で高い次数 (最大次数 = 段数 × 2) の公式が実現できる。
2. 安定性に優れている。A-安定な公式は陽的 Runge-Kutta 法では実現不可能である
但し，陰的 Runge-Kutta 法は，非線型常微分方程式の場合は上で述べたように非線型方程式を解く必要があり，Newton 法等の反復解法 (以下，これを内部反復と称する) を使用することになる。
このため，アルゴリズムは陽的 Runge-Kutta 法より格段に複雑になる。しかし，特に安定性に優れているために，刻み幅が大きく取れ，Stiﬀ 問題向きの解法と言える。
線型常微分方程式であれば，陰的 Runge-Kutta 法でも内部反復を必要とせず，連立一次方程式を

1 度だけ解けばよい。また，刻み幅を小さくすれば，この連立一次方程式の解の存在も保証されている [19]。
以下，陰的 Runge-Kutta 法の内部反復について触れる。

16.5 陰的 Runge-Kutta 法の内部反復
逐次的に k1 から km まで計算できる陽的解法に比べ，陰的解法は非線形方程式 (16.11, 16.10) を解く必要がある。この非線型方程式を解く部分を内部反復と呼ぶ。

16.5.1 Newton 法を用いた内部反復
Newton 法のアルゴリズムを以下に示す。
アルゴリズム 39 (内部反復のための Newton 法)(0)
1. k(0)1 , ..., km を決める。

2. for l = 0, 1, 2, ...
(a)


 (l)
 (l+1)   (l) P k1 − f(xi + c1 hi , yi + hi mj=1 a1 j k(l)  k1  k1
) 
j


 (l)P (l+1)   (l)  k2 − f(xi + c2 hi , yi + hi mj=1 a2 j k(l))
  k2  k2j .  :=  .  − J −1 (k(l) , ..., k(l)
)

..
m 
1
 .   . 



.

 .   .  (l)Pm (l+1)   (l) 
(l) 
km
km
km − f(xi + cm hi , yi + hi j=1 am j k j )(b) 収束判定(l)(l)
ここで， J(k(l)
1 , k2 , ..., km ) は

 In − J11 −J21(l (l)(l)J(k1 , k2 , ..., km ) = 
..

.


−Jm1

−J12

···

−J1m

In − J22
..
.

···

−J2m
..
.

−Jm2

···

In − Jmm










第 16 章 常微分方程式の初期値問題 — 一段法

192
但し，

J pq = hi a pqmX∂f(xi + c p hi , yi + hi
a p j k(l)j ) ∈ Mn (R)
∂y
j=1

(l)
実際の計算では， J(k1(l) , k2(l) , ..., km) の逆行列を求めるのではなく，この行列を LU 分解しておき


 

P
 z1   k1(l) − f(xi + c1 hi , yi + hi mj=1 a1 j k(l)
) 
j


  (l) k − f(xi + c2 hi , yi + hi Pm a2 j k(l) ) 

z


2
j=1
  2
j

(l)

(l) 
J(k(l)
..
 ..  = 1 , k2 , ..., km ) .
 .  Pm  (l)(l) zm
km − f(xi + cm hi , yi + hi j=1 am j k j )という mn 次連立一次方程式を [z1 z2 ... zn ]T について解く。準 Newton 法を用いる場合は，LU 分解
は反復の前に一度だけやっておき，後は LU 分解された Jacobi 行列を用いて，後退代入だけ行う[16]。

16.6 陽的・陰的 Runge-Kutta 法の比較
一般に，陽的 Runge-Kutta 法は，陰的 Runge-Kutta 法に比べて安定領域が狭く，Stiﬀ な問題には向いていないと言われている。しかし，安定領域が狭いという件は，刻み幅を狭く取る必要があるということしか意味しない。前節で示した通り，陰的 Runge-Kutta 法では (16.11) 式のような非線型方程式を解く必要があるため，反復計算が増え，1step ごとの計算量は陽的 Runge-Kutta 法よりも多くなる。実際には，陽的・陰的のどちらがより少ない時間で実行できるのか？
もう一つの疑問は，同じ規格の浮動小数点数を用いた場合，陽的・陰的 Runge-Kutta 法の限界精度がどの程度違うのか，ということである。
これらの観点から，陽的 Runge-Kutta 法と陰的 Runge-Kutta 法とを比較してみる。

16.6.1 線型常微分方程式における陰的 Runge-Kutta 法のアルゴリズム
ここでは次のような定数係数線型常微分方程式を考える。

 dy



= Ay + g(x) = f(x, y)




 dx







y(x0 ) = y0

(16.12)

ここで，A ∈ Mn (R) である。
この場合，(16.10)，(16.11) 式はそれぞれ n 次，mn 次の連立一次方程式となる。
半陰的解法の場合(I − hi arr A)kr = f(xi + cr hi , yi + hi

r−1
X
j=1

ar j k j )(16.13)

16.6. 陽的・陰的 Runge-Kutta 法の比較
陰的解法の場合 I − hi a11 A
 −h a A
i 21


..

.

−hi am1 A

193−hi a12 A
..
.
..
.

···
..
.
..
.

···

−hi am,m−1


 

  k1   f(xi + c1 hi , yi )   .  ..
  ..  .
 
 =   .  .
..
−hi am−1,m A   ..  

 
 

I − hi amm Akmf(xi + cm hi , yi )−hi a1m A..
.

(16.14)16.6.2 数値実験 — 固定刻み幅における線型常微分方程式
ここでは数値実験により，陽的 Runge-Kutta 法，陰的 Runge-Kutta 法の精度の比較を，線型常微分
方程式に対して，固定刻み幅で計算することによって行う。計算は全て 2 進 53 桁倍精度で行った。
まず，ここで使用する陰的公式の係数を示す。

2 段 4 次 : Butcher[15]

1 段 2 次 [15]
1
2

√
3− 3
6√
3+ 3
6

1
2

1

3 段 6 次 : Butcher[15]
√
5− 15
10
1
2√
5+ 15
10

√
10−3 15452
9√
10+3 15
45
8
18

5
36√
10+3 1572√
25+6 15
180
5
18

1
4√
3+2 3
12
1
2

√
3−2 3
12
1
4
1
2

√
25−6 15180√
10−3 15
72
5
36
5
18

陽的公式の係数を示す。

2 段 2 次 : Optimal[15]
2
3

2
3
1
4

3
4

4 段 4 次:古典的 (classical)Runge-Kutta 法 (アルゴリズム 37)[15]
1
2
1
2

1
2

0

1
2

1

0

0

1

1
6

1
3

1
3

1
6

7 段 6 次 : Butcher[15]
1
3
2
3
1
3
1
2
1
2

1

1
3

0
1
12
1
− 16

0
9
44
11
120

2
3
1
3
9
8
9
8
9
− 11

0

1
− 12
3
− 16
− 38

− 38
− 34

63
44
27
40

18
11
27
40

1
2

0

− 16
11

4
− 15

4
− 15

11
120

第 16 章 常微分方程式の初期値問題 — 一段法

194

Butcher の係数は陰的・陽的にかかわらず様々な段数のものが存在するが，ここでは Jain[15] の値を使用した。
これらの係数を利用して次の２つの例題を解いてみる。この問題を取り上げたのは，どちらも同じ解析解でありながら，前者に比べて後者の Stiﬀ 性が格段に高く，比較対照しやすいためである。
例題 16.6.1

 dy
1



dx


 dy2


dx





=

−2y1 + y2 − cos x=

2y1 − 3y2 + 3 cos x − sin x
y1 (0) = 1 y2 (0) = 2(16.15)

例題 16.6.2

 dy
1



dx


 dy2


dx





=
=

−2y1 + y2 − cos x
1998y1 − 1999y2 + 1999 cos x − sin x(16.16)

y1 (0) = 1, y2 (0) = 2(16.15)，(16.16) とも三井 [20] の問題である。どちらも同じ解析解 (16.17) を持つ。
 y1 (x) = exp(−x)
 y2 (x) = exp(−x) + cos x(16.17)

しかし，後者は Stiﬀ になっており，陽的解法では刻み幅を小さくしないとオーバーフローを起こす。双方の計算の結果を表 16.2〜16.5 に示す。三井の数値実験と同じく，積分区間を [0, 20] とし，この区間を等間隔に分割し，最後の点 (xl = 20) での数値解の相対誤差を並べたものである。
計算例 16.6.1(表 16.2, 16.3) では両者とも目立った相違はない。計算例 16.6.2(表 16.4, 16.5) では，
安定領域外の刻み幅で陽的解法がオーバーフローしており，限界精度に達するのは陰的解法に比べて若干遅い。しかし，限界精度は両者の差はない。
問題 16.6.1
上の二つの線型常微分方程式 (16.15)(16.16) の解析解が (16.17) であることを確認せよ。

16.7 Rössler モデルの数値例
Rössler Model はカオス現象が見られる比較的簡単な３次元力学系の一つである [33]。

  

 x   −(y + z)




d    y  =  x + αy
dt   
β + z(x − µ)z(16.18)

ここでは α = β = 1/5 と固定して考える。µ を 3, 4, 5 と増やしていくと，この常微分方程式の解の周期が増加し，その運動はカオスになることが知られている [33]。

16.7. Rössler モデルの数値例

195

表 16.2: 計算例 16.6.1 : 陽的 Runge-Kutta 法

2-2:optimal
刻み幅

4-4:classical

7-6:Butcher

y1

y2

y1

y2

y1

y2

2

1/2
1/24

1.080e+06
1.639e+03

3.076e-02
1.113e-03

9.948e+04
2.415e+02

1.062e-03
2.634e-06

6.414e+03
9.081e-01

6.502e-05
9.221e-09

1/26
1/28
1/210

2.977e+022.373e+01
1.559e+00

6.351e-05
3.887e-06
2.417e-07

8.362e-01
3.168e-03
1.241e-05

9.178e-09
3.484e-11
1.347e-13

1.920e-04
8.164e-08
9.176e-08

1.951e-126.801e-16
1.632e-15

1/212
1/214

9.861e-02
6.182e-03

1.509e-08
9.425e-102.902e-07
7.094e-07

2.721e-15
3.401e-15

2.190e-07
6.971e-07

2.993e-15
3.129e-15

1/216
1/218

3.860e-04
2.616e-05

5.890e-11
3.712e-128.829e-07
1.531e-06

1.714e-14
2.272e-14

8.942e-07
1.550e-06

1.714e-14
2.299e-14

表 16.3: 計算例 16.6.1 : 陰的 Runge-Kutta 法

1-2

2-4:Butcher

3-6:Butcher

刻み幅

y1

y2

y1

y2

y1

y2

1/22
1/24
1/26

4.552e+04
1.459e+03
1.737e+02

8.924e-03
6.224e-04
4.046e-05

2.572e+03
1.434e+01
6.044e-02

2.806e-05
1.756e-07
7.790e-105.490e+00
1.266e-03
3.142e-07

5.031e-08
1.290e-11
3.023e-15

1/28
1/210
1/212

1.221e+01
7.849e-01
4.939e-02

2.556e-06
1.602e-07
1.002e-08

2.405e-04
9.966e-07
2.403e-07

3.145e-12
1.099e-14
3.263e-15

6.102e-08
9.086e-08
2.327e-07

0.000e+00
1.493e-15
3.263e-15

1/214
1/216

3.093e-03
1.938e-04

6.264e-10
3.915e-11

7.208e-07
9.044e-07

3.264e-15
1.714e-14

7.173e-07
8.921e-07

3.400e-15
1.687e-14

1/218

1.110e-05

2.427e-12

1.513e-06

2.244e-14

1.516e-06

2.244e-14

第 16 章 常微分方程式の初期値問題 — 一段法

196

表 16.4: 計算例 16.6.2 : 陽的 Runge-Kutta 法

2-2:optimal
刻み幅

4-4:classical

7-6:Butcher

y1

y2

y1

y2

y1

y2

2

1/2

NaN

NaN

NaN

NaN

NaN

NaN

4

1/2
1/26
1/28

NaN
NaN
NaN

NaN
NaN
NaN

NaN
NaN
NaN

NaN
NaN
NaN

NaN
NaN
NaN

NaN
NaN
NaN

1/210
1/212

6.776e-01
2.891e-04

7.096e-06
1.316e-08

1.079e-02
1.856e-05

1.089e-07
1.873e-10

2.089e-03
2.122e-07

2.108e-08
2.147e-12

1/214
1/216
1/218

3.396e-05
2.313e-06
1.485e-07

6.622e-10
3.946e-11
2.441e-12

6.203e-08
2.109e-09
2.909e-09

6.146e-13
2.312e-15
6.801e-16

1.263e-09
2.349e-09
2.373e-09

4.081e-16
0.000e+00
6.801e-16(“NaN”はオーバーフローしたことを表わす。)

表 16.5: 計算例 16.6.2 : 陰的 Runge-Kutta 法

1-2
刻み幅

2-4:Butcher

3-6:Butcher

y1

y2

y1

y2

y1

y2

2

1/2
1/24

7.115e+021.560e+01

1.459e-02
4.889e-04

1.297e+026.593e+00

1.098e-03
6.229e-05

4.998e-01
1.685e-03

4.257e-06
1.591e-08

1/26
1/28

1.174e+00
7.664e-02

3.055e-05
1.909e-06

5.851e-02
2.371e-04

5.798e-07
2.382e-09

5.933e-06
1.363e-08

5.879e-11
1.241e-13

1/210
1/212
1/214

4.841e-03
3.034e-04
1.898e-05

1.193e-07
7.458e-09
4.661e-109.316e-07
3.888e-09
1.621e-09

9.388e-12
3.671e-14
1.360e-16

1.797e-10
1.746e-10
1.773e-09

1.222e-15
1.360e-16
1.360e-16

1/216
1/218

1.186e-06
7.006e-08

2.913e-11
1.820e-12

2.162e-09
2.511e-09

0.000e+00
6.801e-16

2.121e-09
2.556e-09

0.000e+00
5.441e-16

16.7. Rössler モデルの数値例

197

0

0
y

5

y

5

-5

-5

-5

0

5
x
[0, 500], h = 1/8192
Explicit: 6th Order

-5

0

5
x
[0, 500], h = 1/8192
Implicit: 6th Order

図 16.2: 6 次の陽的 R-K 法 (左) と陰的 R-K 法 (右) の数値解 (µ = 4)

0

0
y

5

y

5

-5

-5

-10

-10
-5

0

5
10
x
[0,500], h = 1/8192
Explicit: 6th Order

-5

0

5
10
x
[0,500], h = 1/8192
Implicit: 6th Order

図 16.3: 6 次の陽的 R-K 法 (左) と陰的 R-K 法 (右) の数値解 (µ = 5)

第 16 章 常微分方程式の初期値問題 — 一段法

198

初期値を [1 0 0]T とし，積分区間を [0, 500] に設定して，6 次の陽的 Runge-Kutta，陰的 Runge-

Kutta 法を使用して (16.18) を計算した。このとき，刻み幅 h を h = 1/8192 とした時の結果を示す。
Fig.16.2 が µ = 4 の時の図，Fig.16.3 が µ = 5 の時の図である。
前述の通り，Rössler Model は µ が 3, 4, 5,... となるにつれて軌道周期が増加し，カオスになる。
数値解もそれに歩調を合わせて，同じ刻み幅であっても次第に精度が悪くなる。実際，µ = 5 の時は陽的・陰的 Runge-Kutta 法のどちらも同じ軌道を再現できていない。
まず，t = 500 の時の x の数値解を，7 段 6 次の陽的 Runge-Kutta 法と 3 段 6 次陰的 Runge-Kutta
法を使用して (16.18) を IEEE754 倍精度浮動小数点数を使用して計算した。その結果を以下に示す。

Explicit : 7 step (6th order), µ = 5
h = 1/4096
h = 1/8192
RN−4.8295189e + 00

−4.6293495e + 00

RZ
RP

−6.0714505e + 00
−4.2363258e + 00

−5.2090559e + 00
−3.5081142e + 00

RM

−5.1634064e + 00

−3.8688926e + 00

Implicit : 3 step (6th order), µ = 5
h = 1/4096
h = 1/8192
RN−3.8470654e + 00

−3.5698751e + 00

RZ
RP
RM

−6.1574318e + 00
−3.4452750e + 00
−3.4044998e + 00

−5.2671008e + 00
−4.6235796e + 00
−3.5246481e + 00

上の表は µ = 5 のときの y1 の値を，RN, RM, RP, RZ モード (第 2 章参照) でそれぞれ計算したものある。浮動小数点演算過程での丸めのモードを変更することによって，丸め誤差の大きさを精度を増やすことなく示すことができる [17]。全ての桁が異なっており，明らかに丸め誤差が数値解を覆い尽くしていることが分かる。
また刻み幅を縮小させても，丸め誤差の order には殆ど変化が無いことも同時に分かる。実際，
もっと小さい刻み幅で計算された数値解でも，order には明らかな変化は見られなかった。従って，
ここで現れた丸め誤差は，単なる計算量に比例する蓄積によるものでは無く，元の常微分方程式(16.18) の性質に依存するものである，と言える。このように，この積分区間 [0, 500] における Rössler
Model の数値計算においては，IEEE754 浮動小数点数倍精度の桁数では不足である。カオスになると言われている µ = 5.7 の場合も同様に，丸め誤差の影響が大きく，精度が必要であれば多倍長計算を行う必要がある。
参考までに z の相対丸め誤差と µ = 5.7 の時の z を描いたグラフを Fig.16.4 に示す。µ の値によって丸め誤差の変動が著しく異なっていることが分かる。µ = 5, 5.7 の時と µ = 4 の時とを比較してみると，単なる丸め誤差の蓄積以上に増大が激しく，微少な桁落ちが随所で起こり，それが丸め誤差の増大を引き起こしていることを予想させる。

16.7. Rössler モデルの数値例

199
Roessler Model: [0, 500], 128bit
Extrapolation 4steps, h=1/1024

Roessler Model: gmp 128bit(38 digits)
Extrapolation: 4 steps
Relative Round-off Error: h=1/1024\mu=5.020

1e-20\mu=5.7|Relative Error|

1e-25
1e-30

z 10
1e-35
\mu=4.0
1e-40
1e-45
1e-50

0
0

100

200

300

400

500

0

100

200

300

400

500

t

t

図 16.4: µ = 4, 5, 5.7 の時の z の相対誤差 (左) と µ = 5.7 の時の数値解 z(右)演習問題
1. 常微分方程式の初期値問題

 dy


 dx


 y(0)

=

x2 y

=

1

に対し，次の問に答えよ。

(a) この常微分方程式の解析解を求めよ。
g をそれぞれ求めよ。
(b) 刻み幅を h = 1/4 とする時，Euler 法と中点法で y(1) の近似解 y(1)g に含まれる相対誤差をそれぞれ求めよ。
(c) 上で求めた y(1)参考図書
常微分方程式を代数的に解くことを解説した書籍は山のようにあるので，どれか一冊は相性の良い物を選んで眺めておくと良い。版を重ねて未だに使用されているものとして
矢野健太郎・石原繁
微分方程式
裳華房, 1994 年を挙げておく。
数値解法を扱ったものは数多くあるが，和書としては
三井斌友・小藤俊幸
常微分方程式の解法
共立出版, 2000 年

第 16 章 常微分方程式の初期値問題 — 一段法200
三井斌友・小藤俊幸・齊藤善弘
微分方程式による計算科学入門
共立出版, 2004 年の 2 冊がよい。前者は線型常微分方程式の理論，Runge-Kutta 法，A 安定性について初学者にも分かりやすく述べており，入門書として優れている。
後者は，常微分方程式の離散解法理論，Hamilton 系の解法，遅延微分方程式・確率微分方程式の解法を扱っており，最新の話題満載の専門書である。
百科辞典的なものとしては

E.Hairer, S.P.Nørsett and G.Wanner
Solving Ordinary Diﬀerential Equations, I
Springer-Verlag, 1993 年
E.Hairer and G.Wanner
Solving Ordinary Diﬀerential Equations, II
Springer-Verlag, 1996 年があるので，日本語版でも出た暁には購入しておくのがよろしかろう。もし原書より安い値で店頭に並んだら，訳者が相当無理をしてコストカットをしたと思った方がよい。なお原書は plain TEX
で組んであり，いざ翻訳しようとして原著者の原稿ファイルを見た TEX 担当者は死ぬ思いをした，という噂を耳にしたことがある。