お話：数値解析 第 5 回
収束の加速法 (前編)長田直樹

1

はじめにC1

収束する数列について何らかの規則性が事前に分0.9かっているときは、規則性を活用することにより、少0.8ない項数の数列から極限値を推定することができる

ADB0.7ことがある。このような方法を収束の加速法あるい0.6は補外法という。加速法を適用すると、少ない項数0.5で極限値の近似値が得られるばかりでなく、加速に0.4より得られる値は、加速しない場合に比べ丸め誤差0.3の影響は少ない。

0.2

今回取り上げる加速法は、エイトケン ∆2 法とリO0.1

チャードソン補外である。これら 2 つの方法は、多0くの加速法の基礎になっている。

−0.6−0.4−0.20x0.20.40.6ν

図 1: 内接正 2 角形の周長の計算

2関孝和の円周率の計算

4CAD に三平方の定理を用いると、(1)(3) より

関孝和 (せきたかかず) がヤコブ・ベルヌーイとほAC =ぼ同じ時期にベルヌーイ数を発見したことを前々回

√

CD

(4)

(2008 年 7 月号) に話した。今回は関の円周率の計算 が得られる。
から話を始める。
関は CD を勾 (こう)、AD を股 (こ)、AC を弦、sν
直径 1 の円に内接する正 2ν 角形の周の長さを sν を周と呼び、ν = 2, 3, . . . , 17 に対し (2)(3)(4) によとする。関は以下に述べる方法で s15 , s16 , s17 まで計 り計算した値を与えた。なお、和算では直角三角形算し、それらの値を用いて円周の長さを推定した。 の直角を挟む短い辺を勾、長い辺を股、斜辺を弦と
図 1 は O を中心、半径 AO = 12 の円弧である。 いう。
∠BOA = π/2ν−2 とすると sν = 2ν AC である。 関の計算した sν の値を見易いように 5 桁毎に空4AOD に三平方の定理を用いると白をつけて示す。イタリックの数字は円周率の真値
1CD2 − CD + AB2 = 04と異なる部分、下線の数字は周の真値と異なる部分(1)である。四捨五入により切り捨てた場合が強、切りが得られる。CD < 12 に注意して (1) を CD につい上げた場合が弱である。とくに末位が 0 を切り捨てて解くと、た場合は微強、末位が 9 を切り上げた場合は微弱でCD =1−√

1 − AB22ある。

(2)となる。D は AB の中点だから

1AD = AB2(3)s2= 2 .82842 71247 46190 0976 微強s3= 3.06146 74589 20718 1738 強中略1s15数列 {sν } の差分を ∆sν = sν+1 − sν 、2 階差分を3.14159 26487 76985 6 708 弱=s16=

3.14159 26523 86591 3 571 強s17=3.14159 26532 88992 7 759 弱

t15= s16 +から

(s16 − s15 )(s17 − s16 )(s16 − s15 ) − (s17 − s16 )3.14159 26535 89793 2476=∆ sν = ∆sν+1 − ∆sν により定義すると、
2(sν+1 − sν )(sν+2 − sν+1 )(sν+1 − sν ) − (sν+2 − sν+1 )(sν+1 − sν )2= sν −sν+2 − 2sν+1 + sν(∆sν )2(7)= sν −∆2 sνtν= sν+1 +(5)を計算しと書ける。数列 {sν } を数列 {tν } に対応させる変換

3.14159 26535 9 微弱

は、エイトケン ∆2 法あるいはエイトケン加速法とを定周とした [6, 三四八頁]。小数第 9 位以下 3589呼ばれている。

を切り上げて 359 微弱としているので、円周率を名称は統計学者の A.C. エイトケン [1] が、1926 年13 桁 (整数部分 1 桁を含む) 書き示したことになに代数方程式の最大根を求める過程で使ったことにる。実際には、t15 は 17 桁正確に求まっている。 由来する。エイトケンの再発見は、関の発見から二{sν } ν=2,3,...,17 は小数第 20 位を丸めて 20 桁表示してあるが、s15 , s16 , s17 の小数 17 位以下は桁落ち (近百数十年後のことになる。∆2 法は (7) の分母から来ている。

い数同士の引き算により有効桁数が少なくなる現象)などのため正しくない。そのため t15 は小数第 16 位3までが正しい数字である。

加速法

上記の計算は、関の没後 1712 年に弟子達により出加速法は、収束する数列 {sν } を同じ極限により速版された括要算法で公表されている。ベルヌーイ数く収束する数列 {tν } への変換としてとらえることと同様に関は (5) の理由を説明していないが、増約ができる。

術と呼ばれていた無限等比級数の和の公式によった

s に収束する数列 {sν } に対し、ものと考えられている [6, pp.183-184]。

limtν − sν→∞ sσ(ν) − ssν+2 − sν+1sν+3 − sν+2== ··· = rsν+1 − sνsν+2 − sν+10 < |r| < 1=0(8)を満たす数列 {tν } を得る方法を加速法という。ここで、σ(ν) は tν を求めるために用いる最大の添字である。例えば、エイトケン ∆2 法の場合 tν はを仮定するとsν , sν+1 , sν+2 から計算するので σ(ν) = ν+2 となる。
sν+k+1 − sν+k = (sν+1 − sν )rk ,k = 1, 2, . . .

が成立する。k = 1, 2, . . . , l について和を取ると

sν+l+1 − sν+1 = (sν+1 − sν )l∑4

エイトケン ∆2 法については次の定理が基本的でrkある。

k=1だから、l → ∞ として増約術を用いるとs = sν+1 +(sν+1 − sν )r1−r

エイトケン ∆2 法定理 1 数列 {sν } が ν → ∞ のときsν = s + c1 λν1 + c2 λν2 + o(λν2 )(6)(9)と漸近表示されるものとする。ここで、s は未知の極となる。ここで s は数列 {sν } の極限である。(6) の限値、c1 , c2 , λ1 , λ2 は未知の定数で 1 > |λ1 | > |λ2 | >右辺の r に (sν+2 − sν+1 )/(sν+1 − sν ) を代入した0 である。このとき、()2λ1 − λ2tν = s + c2
λν2 + o(λν2 )λ1 − 1(sν+1 − sν )(sν+2 − sν+1 )tν = sν+1 +(sν+1 − sν ) − (sν+2 − sν+1 )を満たす。

は関が定周を求めた式 (5) になっている。

2(10)証明 (9) より

5∆sν = c1 λν1 (λ1 − 1) + c2 λν2 (λ2 − 1) + o(λν2 )関の弟子の建部賢弘 (たけべかたひろ) は、1722年に著した綴術算径 (てつじゅつさんけい) においが成り立つので(∆sν )2∆2 sν

建部賢弘の円周率の計算て、円周率の求め方を解説し、円周率を 42 桁与えた。最後の 1 桁を除き正確である。同書には、師匠= c21 λ2ν
(λ1 − 1)2[1( )ν
λ2 − 12c2 λ2× 1+c 1 λ1
λ1 − 1(( )ν )]λ2+oλ1であった関の方法論との違いをのべた興味深い記述もある。建部は円に内接する正 2ν 角形の周長 sν の2 乗 σν = s2ν を計算したと述べているが、説明の通りに計算を行っても 38 桁しか正しい値は得られない。そのため、sν を用いて円周率 41 桁を得たと考= c1 λν1 (λ1 − 1)2[( )ν ()2c2 λ2λ2 − 1× 1+c1 λ1
λ1 − 1(( )ν )]λ2+oλ1えられている。詳細は [2] を見よ。
建部は sν の差分の比 (sν+2 − sν+1 )/(sν+1 − sν )が 1/4 に近づくのを見て取り、(ν)

s1= sν+1 +がいえる。これらより、{sν } が公比 r の無限等比級数の部分和になると仮定
すると、sν = a(1 + r + · · · + rν−1 ) の極限は増約術により a/(1 − r) である。sν , sν+1 , r を用いて表すと

asν+1 − rsνsν+1 − sν== sν+1 + 1
1−r1−rr −1¤となる。
(ν)

建部はさらに s1

系 1 数列 {sν } が定理 1 の条件を満たすとき、エイの差分の比が 1/16 に近づくことから、トケン ∆2 法は収束を加速する。

(ν+1)

(ν)

s2

証明 定理 1 よりlimtν − sν→∞ sν+2 − s=0s1+(ν)− s1
16 − 1を計算した。そして、s0 = 2 とし、ν = 2, 3, . . . , 10に対し(ν)

s0∞∑π(−1)j π 2j+1 −2j νsν = 2ν sin ν = π +(2 )2(2j + 1)!
j=11
16= sν(ν−k+1)(ν−k)

sk(ν−k+1)= sk−1+k = 1, . . . , ν − 1

より、λ1 = 1/4, λ2 = 1/16, c2 = π 5 /5! とおくと((ν+1)= s1(1)¤例 1 関の方法は定理 1 により正当化される。

π5tν − π ;5!

(11)を計算した。(11) 式は次のように説明できる。数列(∆sν )2∆2 sν
[( )ν
()
c2 λ2
λ2 − 1λ2 − 1= c1 λν1 1 +2−c1 λ1
λ1 − 1
λ1 − 1(( )ν )]λ2+oλ1が導けるので、漸近公式 (10) が得られる。

sν+1 − sν4−1)ν+1sk−1(ν−k)− sk−1(12)2k2 −1(1)を計算し、s9 を円周率とした。[2, pp.6-34](12) によって収束を加速する方法は、リチャード

ソン補外あるいはリチャードソン加速法と呼ばれてが導ける。関の結果の理論上の誤差は ν = 15 より、 いる。これを k = 1 の場合に最初に適用したのは、土星の衛星を発見した天文学者で物理学者の C. ホイ

π 5 −1616; 1.382 × 10−195!

ヘンスである。1654 年に円周率の計算において (11)を適用した。リチャードソン補外の名称の由来は、数である。

値予報のパイオニアである L.F. リチャードソン [4]3定理 2 関数 F (h) が漸近展開 (13) を満たすときが次節で述べるような提案を行ったことによる。建部はリチャードソンよりも約 200 年早くリチャードFk (h)ソン補外を用いた。

∞∑∼ s+( k)∏ 22i−2j − 1cj

i=1j=k+1

6

リチャードソン補外

6.1

連続量の補外∞∑(ν)

Tk(h → +0)

cj h2j ,h2j (15)

証明 k に関する数学的帰納法で証明できる。¤= Fk (h/2ν ) とおくと、(ν)

T0

h > 0 で定義された関数 F (h) が漸近展開
F (h) ∼ s +22i − 1h= F( ν )2(ν−k+1)

(ν−k)

(13)

(ν−k+1)Tk= Tk−1+Tk−1

j=1(ν−k)− Tk−12k2 −1k = 1, . . . , νを満たしているとする。ここで、s は未知の極限値、(ν)と書ける。2 次元配列 {Tk } は

c1 , c2 , . . . は未知の定数である。(13) のような漸近展開を持つ関数は、h > 0 を刻み幅とする数値積分公(0)T0&

式や数値微分、常微分方程式の初期値問題にしばしば登場する。

(1)
T0→&

T1

T0(2)→

T1(3)
T0&→(2)
T1

h > 0 を適当に取る。
F (h) ∼ s +∞∑cj h2jj=1∞∑hF( ) ∼ s +cj 2−2j h2j2j=1(0)&(1)→(0)T2&→&(1)T2→(0)

T3のように計算する。この表は T -表と呼ばれる。

から c1 の項を消去することを考える。

数列の補外

6.2hF ( ) − 2−2 F (h)2∞∑cj (2−2j − 2−2 )h2j
∼ (1 − 2−2 )s +数列 {sν } が漸近展開sν ∼ s +j=2∞∑cj λνj ,(1 > |λ1 | > |λ2 | > . . . > 0)

j=1よりF1 (h)=(16)を満たしているとする。ここで、c1 , c2 , . . . は未知の定数で、λ1 , λ2 , . . . , は既知の定数とする。数列 {sν }F ( h2 ) − 2−2 F (h)1 − 2−2F ( h2 ) − F (h)h= F( ) +2
22 − 1に対するリチャードソン補外は(ν)

T0

とおくと

F1 (h) ∼ s +∞∑− 1 2jh22 − 1(ν−k)

Tk2−2jcj

j=22= sν(14)(ν−k+1)= Tk−1+)λk ( (ν−k+1)(ν−k)
Tk−1− Tk−1
1 − λkk = 1, . . . , νとなる。リチャードソンは、F (h) から F1 (h) を作るである。

過程を h2 -補外と名づけた。

定理 3 漸近展開 (16) を満たす数列 {sν } にリチャー

この過程を繰り返し

ドソン補外を適用すると

F0 (h) = F (h)

(ν−k)

Tk

Fk−1 ( h2 ) − Fk−1 (h)hFk (h) = Fk−1 ( ) +2
22k − 1

∼ s+∞
∑
j=k+1

を F (h) に対するリチャードソン補外という。

4cj( k)∏ λj − λi
1 − λii=1λjν−k証明 定理 2 と同様である。¤

数列 {sν } が (9) を満たしているとき、∆2 法 tν−2の誤差は定理 1 より c2 ((λ1 − λ2 )/(λ1 − 1))2 λ2ν−2系 2 λj = 2−2j (j = 1, 2, . . .) の場合、ν → ∞ とす(ν−1)

でリチャードソン補外 T1

ると(ν−k)Tk= s+(−1)k ck+1 2(k+1)(k−2ν) +o(2(k+1)(k−2ν) )例 3 sν = 2ν sin(π/2ν ) にエイトケン ∆2 法とリチ

例 2 sν を直径 1 の円に内接する正 2ν 角形の周の長ャードソン補外を適用した際の誤差は表 1 のようにさとする。建部の方法は、(ν−1)

なる。tν−2 − s ; −4(T1∞∑πsν = 2 sin( ν ) = π +cj λνj2j=1
ν2j+1cj = (−1)jπ,(2j + 1)!の誤差は定理 3 より

c2 ((λ2 − λ1 )/(1 − λ1 ))λν−1である。λ1 = 1/4, λ2 =2(ν−1)
1/16 のとき tν−2 − s ; −4(T1− s) となる。

− s), (4 ≤ ν ≤ 8) となっている。計算は gcc 4.1.2 の倍精度 (10 進 16 桁弱) で行った。

λj = 2−2j

表 1: sν = 2ν sin(π/2ν ) を加速(ν−1)にリチャードソン補外を適用したものである。理論νsν − πtν−2 − πT1

上の誤差は系 2 より4−2.0 × 10−2

6.4 × 10−4−1.6 × 10−4567−5.0 × 10−3−1.3 × 10−3−3.2 × 10−43.9 × 10−5
2.4 × 10−6
1.5 × 10−7−9.7 × 10−6−6.1 × 10−7−3.8 × 10−88−7.9 × 10−5

9.5 × 10−9−2.4 × 10−98

ロンベルク積分法(1)

T9− π ; (−1)9 c10 2−110 ; −4.154 × 10−43である。

7−π∆2 法とリチャードソン補外
数列 {sν } がsν ; s + c1 λν

0 < |λ| < 1

関数 f (x) は閉区間 [a, b] で C ∞ 級とする。区間∫b[a, b] を n 等分した際 I = a f (x)dx に対する複合台(17)を満たしていると仮定する。sν に対するリチャード形公式は h = (b−a)/n, xj = a+jh(j = 0, 1, . . . , n),

n−1
∑
1
1
Tn = h  f (a) +f (xj ) + f (b)22

ソン補外の最初のステップは(ν)

T1= sν+1 +λ(sν+1 − sν )1−λ(18)j=1である。一方、である。

sν+2 − sν+1;λsν+1 − sν

オイラー・マクローリンの公式によりTn ∼ I +より、(18) の λ を (sν+2 − sν+1 )/(sν+1 − sν ) で置き∞∑
B2j (j=1換えると、(2j)!

)
f (2j−1) (b) − f (2j−1) (a) h2jが成立する。したがって、
sν+2 − sν+1(sν+1 − sν )T1 , T2 , T4 , T8 , T16 , . . .
sν+1 − sνsν+1 +sν+2 − sν+11−にリチャードソン補外を適用することができる。
sν+1 − sνν = 0, 1, . . . に対し(sν+2 − sν+1 )(sν+1 − sν )= sν+1 +(19)(sν+1 − sν ) − (sν+2 − sν+1 )(ν)
T0 = T2ν
(ν−k+1)(ν−k)となりエイトケン ∆2 法が得られる。リチャードソ
Tk−1− Tk−1(ν−k)(ν−k+1)
T=T+(20)kk−1ン補外は収束率 λ が未知の場合は利用できないが、22k − 1
k = 1, . . . , ν
その場合でもエイトケン ∆2 法は利用できる。

5(ν)

W. ロンベルク (1909-2003) に生前インタビューを行った藤野清次九大教授は、ロンベルクがロンベルにより {Tk } を計算する数値積分公式をロンベルク
積分法 [5] という。
(0)

(0)

− Tν−1 | < ²

ク積分法を思いついた際のエピソードを紹介してい(0)が成立するとき計算を停止し、Tν を定積分の近似る [7, p.75]。ある積分を複合シンプソン公式で計算値に採用するしようとしたが、当時 (1950 年代前半) 使えるコン与えた許容誤差 ² > 0 に対し、|Tν

ピュータの記憶容量は 30 語しかなく、精度が足りな定理 4 (ロンベルク積分法の誤差)かった。次に、ガウス型公式で計算しようとすると

関数 f (x) は閉区間 [a, b] で C ∞ 級とする。(20) は)( k∞∏ 22i − 22j ( b − a )2j
∑
(ν−k)
Tk
∼I+cj
22i − 12νj=k+1

係数を記憶させるだけで記憶容量が一杯になった。
そこで、複合台形公式を加速して精度を高める方法を思いついた。

i=1

ロンベルクの発見のきっかけにもなっているようを満たす。ここで)B2j ( (2j−1)
cj =f(b) − f (2j−1) (a)(2j)!に、リチャードソン補外は記憶容量の使用を少なく(ν−k)

することができる。そのためには、Tkの値を(ν−k)
Tk−1 に上書きするようにプログラミングする。

である。

ロンベルク積分法は数値積分の標準的な算法であ証明 定理 2 より得られる。¤

る。例 4 で用いた C 言語のプログラムは系 3 ν → ∞ のとき、次の漸近公式が成立する。

http://www.cis.twcu.ac.jp/~osada/rikei/
においておく。

(ν−k)

Tk= I + (−1)k ck+1 2(k+1)(k−2ν) (b − a)2k+2
例4+o(2(k+1)(k−2ν) )

参考文献

∫ 1

[1] A.C.Aitken, On Bernoulli’s numerical solutionof algebraic equations, Proc. Roy. Soc. Edin-xe dx = 1.71828 18284 59045
0burgh Ser A 46(1926), 289–305.にロンベルク積分法を適用した結果を表 2 に示す。

[2] 小川・平野、数学の歴史、朝倉書店、2003∫1
表 2: 0 ex dx にロンベルク積分法を適用(ν)[3] 長田直樹、数値微分積分法、現代数学社、1987(ν−1)νT00121.85914 09142 29523
1.75393 10924 64825
1.72722 19045 57517

1.71886 11518 76593
1.71831 88419 21747

Roy. Soc. London Ser. A 226(1927), 299-349.

31.72051 85921 64302

1.71828 41546 99897[5] W. Romberg, Vereinfachte numerishe Integra-

T1

(ν−2)

νT3

2

1.71828 26879 24757

3

1.71828 18422 18440(ν)

tion, Kgl. Norsk. Forhandinger, 32(1955), 30-36

(ν−3)

T2

T0

[4] L.F. Richardoson, The diferred approach to thelimit, Part I – Single Lattice, Philos. Trans.

1.71828 18287 94530(ν−1)

は複合台形公式、T1[6] 平山・下平・広瀬共編、関孝和全集、大阪教育図書、1974[7] 藤野清次、数値計算の基礎、サイエンス社、1998

は複合シンプソン公(0)の誤差は 3.35 ×10−10 であるので、9.5 桁正確な値が得られている。
(0)
T3 の系 3 による理論上の誤差は式である。8 分割 (ν = 3) で T3

B8 −122 (e − 1) ; 3.47 × 10−108!
である。

6(おさだなおき/東京女子大学)