お話：数値解析 第 1 回
収束の種類
長田直樹

はじめに

1

今日では、科学・技術にとどまらず、日常生活に問題

表 1: 反復法の例
方法

密着する多くの分野でも、数値計算が重要な役割を非線型方程式

ニュートン法果たしている。

n 次代数方程式

ワイエルシュトラス法

数値計算とは、モデル化などにより得られた数学

連立 1 次方程式

ヤコビ法の問題の近似解 (数値解) を計算機によって求めるこ固有ベクトルベキ乗法とをいう。数値計算を支える数学的理論が数値解析である。

表 2: 離散化法の例
問題
方法

数値解析の研究や数値計算を行う上で役に立ちそうな数値解析の話題を、毎号読み切りで 1 年間連載微分する。

定積分

複合台形公式

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

オイラー法

本連載のサポートページを差分法

http://www.cis.twcu.ac.jp/~osada/rikei/におく。数値例に使ったプログラムのいくつかは、こ数列の収束の速さ3のページで公開していく。

ここでは、数列 {xν } は実数列または複素数列とする。Rn 、Cn 、あるいはバナッハ空間の数列に対し数値計算と数列2ても、絶対値を適当なノルムに置き換えることにより、ほぼ同様の議論を展開できる。

数値計算において、数列の極限を求めることに帰着させる場合がしばしばある。数列は、級数の部分和、反復法、離散化法、などにより得られる。表 1

物理的速さとの対比

3.1に反復法の例、表 2 に離散化法の例を掲げる。離散化法では、刻み幅を規則的に小さくする (例えば物理的な速さは、単位時間当りの移動距離で定義

1, 1/2, 1/4, 1/8, . . .) ことにより、数列が得られる。

される。数列の収束の場合に当てはめると、移動距離に相当するのは得られた有効桁数 (誤差の絶対値数列とは自然数全体の集合 N から集合 X への写像である。X が実数全体の集合 R のときは実数列、複素数全体の集合 C のときは複素数列、Rn のときは実ベクトル列などという。の常用対数をとり符号を − にしたもの)、単位時間は反復 1 回と考えられる。
近似値が x で真値が α のときの有効桁数は− log10 |x − α| なので、xν における収束の速さを近似値を数列の極限値として求めるアルゴリズム
では、何らかの条件を満たしたときに収束したものvνと見なし、計算を終了する。これらのアルゴリズムを比較する場合、収束の速さは重要な要素になる。
今回は収束の速さの分類とそれぞれの特徴を述べ= − log10 |xν+1 − α| + log10 |xν − α|¯¯¯ xν+1 − α ¯¯(1)= − log10 ¯¯xν − α ¯により定義する [1, p.51]。

る。

1

3.2

数列の収束の分類

数列 {xν } は α に収束し、誤差の比の極限値

ρ = limν→∞

1 点反復列の収束

4xν+1 − αxν − αα に収束する数列 {xν } が、関数 φ(x) によって、(2)が存在するとする。このとき、|ρ| ≤ 1 が成り立つ。

[証明せよ。] 0 < |ρ| < 1 のとき {xν } は α に収束率(あるいは縮小率)ρ で線型収束するという。ρ = 0 のときは超線型収束、ρ = 1 のときは対数収束あるいxν+1 = φ(xν ) と表される場合を考える。φ(x) を 1点反復関数 [3, p.8] という。
φ(x) が α で微分可能であれば
lim

ν→∞

xν+1 − α= φ′ (α)
xν − αとなる。このとき、数列 {xν } が α に収束率 φ′ (α) で線型収束することと 0 < |φ′ (α)| < 1 は同値である。

は劣線型 (sublinear) 収束するという。

定理 1 反復関数 φ(x) は α の近傍で C p 級 (複素関数数値計算は多くの分野にまたがっており、用語も分野によって異なることがしばしばある。対数収束は収束の加速法、劣線型収束は主として非線型方程式の反復解法の分野で用いられている。の場合は正則) であるとする。xν+1 = φ(xν ) によって生成される数列 {xν } に対し、次は同値である。

1. φ′ (α) = . . . = φ(p−1) (α) = 0, φ(p) (α) ̸= 0

日本語として定着してない用語には、括弧書きで英語をつけてある。

2. α に p 次収束し、漸近誤差定数が |φ(p) (α)|/p! と数列 {xν } が収束率 ρ で線型収束するとき、収束なる。の速さの極限は証明テーラーの定理

1lim vν = − log10 |ρ| = log10
ν→∞|ρ|

φ(xν )= α+となる。したがって、ρ が 0 に近いとき収束は速く、p−1 (k)
X
φ (α)
k=1k!

(xν − α)k + RνRν → 0 (ν → ∞)1 に近いとき収束は遅い。

より導ける。[証明を完成せよ。]

3.3

収束次数と漸近誤差定数¤

定理 1 より、1 点反復列の収束次数は自然数になる。

数列 {xν } が α に線型収束または超線型収束する
ものとする。実数 p ≥ 1 と零でない極限値|xν+1 − α|
ν→∞ |xν − α|p

C = lim

数値例

5関数 f (x) = (x − 1)2 (x − 2) = x3 − 4x2 + 5x − 2(3)の零点 (f (α) = 0 となる α を f (x) の零点という。)
が存在するとき、{xν } は α に p 次収束するという。 を、収束の種類の異なる 4 通りの反復法で求める。

C は漸近誤差定数と呼ばれる。
|f (xν )| < 1.0 × 10−15 となったときに反復を終了
数列 {xν } が α に p 次収束するとき、正数 M と自 する。
然数 N が存在してコンパイラは gcc 4.1.2 を用い、倍精度 (10 進 16|xν+1 − α| ≤ M |xν − α|p(ν ≥ N )桁弱) 計算を行った。非線型方程式の反復解法につ(4)いては、回を改めて説明する。となる。
収 束 次 数 p の 近 似 値 は 、連 続 す る 3 項

5.1

xν−1 , xν , xν+1 から
log10 |(xν+1 − α)/(xν − α)|log10 |(xν − α)/(xν−1 − α)|3 次収束の例

例 1 初期値を x0 = 2.1 にとりハレー法(5)f (xν )f ′ (xν ),xν+1 = xν −f (xν )f ′′ (xν )1−2(f ′ (xν ))2により与えられ、計算収束次数 (computational orderof convergence) と呼ばれる。[(3) から (5) を導け。]
2

ν = 0, 1, . . . (6)表 3: ハレー法：x0 = 2.1
f (xν )

νxν

2.10

1.2 × 10−1

00.90

−1.1 × 10−1

122.0021
2.000000026

2.1 × 10−3
2.6 × 10−8

120.950.97−2.9 × 10−3−7.3 × 10−432.000000000000000

0.0

ν

xν0

中略2122

表 4: ニュートン法：x0 = 2.1
f (xν )

νxν

02.10

0.12

122.0152.000450.016
0.000453452.00000041
2.00000000000033
2.000000000000000

0.999999948
0.999999974−74.1 × 10
3.3 × 10−13
0.000.900−1.1 × 10−2

120.9470.960−2.8 × 10−3
−1.6 × 10−3中略

収束は非常に速く x3 が真値を与えている。倍精度

1000000
2000000

0.99999977
0.99999989−5.3 × 10−14−1.3 × 10−14

4000000

0.999999943
中略0.999999968−3.3 × 10−15で計算しているため、表 3 からは xν+1 − 2 と xν − 27273219の関係は読み取れない。

例 4 初期値を x0 = 0.9, 係数 A = 1/f ′ (x0 ) として

例 2 初期値を x0 = 2.1 としたニュートン法ν = 0, 1, . . .

簡易ニュートン法(7)xν+1 = xν − Af (xν ),による結果は表 4 である。

ν = 0, 1, . . .

(8)を適用すると表 6 のようになる。

|xν+1 − 2| ≤ 2|xν − 2| ,2|f (x7273219 )| < 1.0 × 10−15 となり ν = 7273219

ν = 0, 1, 2, 3.

で反復が終了する。収束は極めて遅い。反復回数が2 倍になると誤差が半分になっている。つまり誤差を満たしている。反復を 1 回繰り返すごとに有効桁xν − 1 が 1/ν の定数倍で近似できることが推察で

数 (= − log10 | 誤差 |) が 2 倍になっている。

5.3−9.9996 × 10−16対数収束の例

5.42 次収束の例

f (xν )
xν+1 = xν − ′,f (xν )

−2.7 × 10−15−6.8 × 10−16

表 6: 簡易ニュートン法：x0 = 0.9
ν xνf (xν )で計算した結果を表 3 に示す。

5.2

表 5: ニュートン法：x0 = 0.9f (xν )きる。

線型収束の例例 3 初期値を x0 = 0.9 としてニュートン法により計算した結果を表 5 に示す。

(xν+1 − 1) ; (1/2)(xν − 1),

63 次収束6.1

二宮市三先生パソコンが普及する以前、数値計算はプログラム

ν = 0, . . . , 20

言語 (主として FORTRAN) でプログラムを書いたとなっているので線型収束である。反復を 1 回繰りり、数値計算ライブラリを利用することによりなされ返すごとに誤差が半分になっている。ていた。プログラム言語で利用する標準関数や数値計算ライブラリは数学ソフトウェアと呼ばれていた。

3数学ソフトウエア開発の大御所に二宮市三先生(名古屋大名誉教授) がおられる。名大在職中に表 7: 8 倍精度計算：x0 = 2.1

NUMPAC という汎用数値計算ライブラリーを作られた方である。二宮先生は定年退職後、Borland C++
上に IEEE754 規格を拡張した 4 倍精度と 8 倍精度

ニュートン法

ハレー法

νxν

xν

01

2.10
2.015

2.10
2.0021232.03 45
2.06 412.07 26
2.022 52456

2.012 33
2.024 222.048 102.066 41のシステムを完成させた。([2, p.6] 参照)
毎年 6 月頃に数値解析シンポジウム (今年は、http://nas2008.akita-pu.ac.jp/) が 開 催 さ れる。数値計算の実務家や数値解析の研究者など 100人以上が、二泊三日で泊まり込み熱心な討論を行う。
本誌の読者にはお馴染の一松信先生も、常連の参加者である。
昨年 (2007 年) のシンポジウムで二宮先生と歓談のおり、先生が筆者に 8 倍精度システムを勧められとなることが見て取れる。

たので、喜んで頂くことにした。数日後に圧縮されx 軸に反復回数、y 軸に log10 |xν −2| をとり、折れたシステムが郵送で届いた。インストールと動作確線で結んだグラフを示す。ハレー法は破線、ニュー認だけして、そのままにしてあったのだが、本原稿トン法は一点鎖線で表している。の執筆に 8 倍精度システムが役立つことになった。

0−106.2

8 倍精度システム−20

二宮先生の 8 倍精度システムは、8 倍精度浮動小−30

数点型 octo を符号なし整数 (32 ビット)8 個からなる配列 (32 バイト =256 ビット) により構成してあ−40る。符号に 1 ビット、指数に 15 ビット、仮数に 240−50ビット (仮数部の先頭ビットが 1 になるように指数を−60決める) で、有効桁数は 10 進約 72 桁 (; log10 2240 )である。

−700123456

実務上のほとんどの数値計算は倍精度で十分であるが、3 次以上の高次収束の様子を見るときなどに8 倍精度での計算では、2 次収束 (ニュートン法)と 3 次収束 (ハレー法) の違いがよく分かる。

は、多倍長計算は重宝である。

6.3

再び 3 次収束の例

参考文献

例 5 例 1 と例 2 の数列を 8 倍精度で計算すると表 7[1] C. Brezinski and R. Zaglia, Extrapolation
methods theory and practice, North-Holland,のようになる。同じ数字 n が d 個続くとき、nd と書くことにする [3, p.271]。例えば、2.00045 を 2.03 45

1991.と書く。本例では、f ′ (2) = 1 なので、f (xν ) ; xν −2[2] 二宮市三編、数値計算のわざ、共立出版、2006となるため、関数値 f (xν ) は省略した。

x2 , x3 , x4 によるニュートン法の計算収束次数は、[3] J.F. Traub, Iterative methods for the solution
2.0044、ハレー法は 3.00043 である。
of equations, Chelsea, 1964
表 7 よりハレー法では(おさだなおき/東京女子大学)|xν+1 − 2| ≤ 3.1|xν+1 − 2|3 ν = 0, 1, 2, 34