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

1

補題 1 ν → ∞ のとき、漸近展開

はじめにsν =これまで 3 回にわたり、円周率およびライプニッ∞ν∑1∑ B2k1
1π2− + 2−∼2i6ν2νν 2k+1i=1
k=1

ツ級数の加速にまつわる収束の加速法の話をしてきた。前編 (2008 年 8 月号) では、線型収束の加速にが成立する。ここで、B2k はベルヌイ数である。

有効なエイトケン ∆2 法とリチャードソン補外につ証明 ν < µ を任意の自然数とする。区間 [ν, ν + µ]

いて、中編と中編の 2(9-10 月号) では、交代級数のにおいて f (x) = 1/x2 にオイラー・マクローリンの加速に有効なオイラー変換、e 変換 (実装は ² 算法)、 公式を適用する。µ 分割した複合台形公式による値
レヴィン t 変換を取り上げた。中編の 2(10 月号) でを Tµ とする。

は、一般化リチャードソン補外とその実装である E算法についても話した。
今回は平方数の逆数の和11−22ν2(ν + µ)2(2k)!
f (2k−1) (x) = − 2k+1x[]ν+µ∫ ν+µ1
11f (x)dx = −= −x ννν+µνTµ = sν+µ − sν +∑∞2i=1 1/i を巡る加速法の話から始める。

2平方数の逆数の和より、

sν+µ − sν ∼自然数の平方の逆数の和∞∑1i2i=1=1
11+ 2 + 2 + ···1223k=1µ → ∞ とすると、(1) より求める結果が得られる。
¤
偶数番目のベルヌイ数を書き出しておく。

を求めることは、17 世紀後半から 18 世紀半ばまでヨーロッパの数学界において大きな問題 (バーゼル問題) であった。この問題に解答を与えたのは L. オイラーである。彼は 1735 年サンクトペテルスブル

表 1: 偶数番目のベルヌイ数

グのアカデミーに提出した『逆数の級数の和についkて』[4] においてB2k∞∑
1i=1i2=π261
111−− 2+νν + µ 2ν2(ν + µ)2()∞∑
11B2k −+ 2k+1+(ν + µ)2k+1
νk(1)
B2k1
166691−
2730

21−
30
7763
14283617−510

41−30943867
798

556610174611−330と解決した。
バーゼル問題の答えを数値的に検証するためオイ3ラーは、オイラー・マクローリンの公式 (本誌 2008∑ν年 7 月号) を用いて sν = i=1 1/i2 の漸近展開を求オイラーの数値計算
補題 1 より∑ B2kπ211= sν + − 2 ++O6
ν 2νν 2k+1mめた。1755 年『微分学教程』[5, p.361] において、次の補題を与えている。

k=11(1ν 2m+3)(2)

として sν の極限値 (= π /6) を計算した。オイラーを満たすとき、対数収束 (2008 年 5 月号) という。級∑ν
数 (1) の部分和 sν = i=1 1/i2 のなす数列は、補題の計算は表 2 のようである。

1 よりが成立する。オイラーは (2) において ν = 10, m = 82

表 2: オイラーの計算 (ν = 10)sν
1/ν

1.54976 77311 66540 690
0.1−1/2ν 2−0.005
1.64476 77311 66540 690
0.00016 66666 66666 666

B2 /ν3lim

ν→∞なので π 2 /6 に対数収束する。
対数収束は一般に収束が遅い。sν =B14 /ν
B16 /ν 17∑ν2i=1 1/i は1000(103 ) 項までから 3 桁、100 万 (106 ) 項までから1.64493 43978 33207 356
中略151
1
π2+ O()−ν+1(ν + 1)2
6 = lim=111ν→∞π2−+O()sν −νν26sν+1 −6 桁しか正しい値は得られない。収束の速さはライプニッツ級数と同程度である。対数収束する数列の極限値を得るには、一般には何らかの加速法が必要

1.64493 40668 48225 335
0.00000 00000 00001 166−0.00000 00000 00000 071になる。
対数収束する数列 {sν } の主な加速法は、大きく 4

1.64493 40668 48226 430つに分けられる。

小数点以下第 19 位を切り捨て (一部は四捨五入し)1. 補正項 (誤差項) の情報を用いて加速する。ているため、小数点以下第 18 位に丸め誤差がある。

2. 部分列 (たとえば {sν }ν=1,2,4,8,16,... ) を加速す

オイラーは計算の最後に「(ベルヌーイ数は) 発散する。

るにもかかわらず、それ (1.644934066848226430) は真の和 (π 2 /6) を与えることは明らかである。」[8] と3. 漸近展開の情報を用いて加速する。

述べている。

4. 数列 {sν } に数列変換を適用する。

オイラーは漸近展開の 1/ν 17 までで打ち切っているが、B18 /ν

19= 0.00000 00000 00000 005 まで用いると、より正確な 1.64493 40668 48226 435 が得られる。真値は 1.64493 40668 48226 43647 である。 5オイラーは漸近展開の 1/ν

17までで打ち切った理補正は由を書いてないが、1.

∑∞補正項の利用真値 = 近似値 + 補正k, k = 2, . . . , 16 の値を 17 桁与えてあ により定義する。(2) の右辺第 2 項以下は補正項でる [5, p.365] ので、17 桁必要だったが余裕を持っ ある。誤差 = 近似値 − 真値 なので補正は誤差の符て 19 桁計算した。
号を変えたものである。
i=1 1/i対数収束級数に対してではないが、交代級数に対2. B18 = 43867/798, B20 = −174611/330 と絶対する補正項を用いた加速は、オイラーより 300 年以値が大きくなる。

上前にインドのマーダヴァによるライプニッツ級数ことなどが考えられる。の和を計算する公式 (本連載 10 月号) が起源である。
オイラー・マクローリンの公式により漸近展開の4係数を求めることができれば、補正項による加速計対数収束の加速算ができる。この方法は、3 節で述べたオイラーのバーゼル問題の数値計算 (表 2) に始まる。

s に収束する数列 {sν } がlimν→∞sν+1 − s=1sν − s26部分列の加速sν =∑ν

2i=1 1/i の漸近列は補題 1 より、

ν −1 , ν −2 , ν −3 , ν −5 , ν −7 , ν −9 , . . .

実用上現れる対数収束数列の多くはsν ∼ s +c1c2c3+ 2 + 3 + ···ννν

がとれるので、(3){gj (ν + l) =あるいはsν ∼ s +c1c2c3+ 4 + 6 + ···2νν
νの型の漸近展開を持つ。

(7)

し、最良の値は小数第 15 位まで示した。
∑ν2
表 3:
i=1 1/i に E 算法を適用(5)なので、(4) の場合である。なお、(5) は漸近展開の係数に未知の π を含んでいるので補正による加速には使えない。
数列 {sν } がsν ∼ s + c1 λν1 + c2 λν2 + c3 λν3 + · · · ,j≥4

算し小数第 12 位を四捨五入し表示している。ただの円に内接する正 ν 角形の周長はππ3π5=π−+− +···ν
3!ν 25!ν 4

j = 1, 2, 3(ν + l)3−2j(4) とおくと E 算法が適用できる。結果は表 3 のようである。以下の表 3-5 では倍精度 (10 進 16 桁弱) で計

バーゼル問題の級数は (3) の場合である。直径 1ν sin(ν + l)−j(6)

(ν)

(ν−1)

(ν−2)

νE0

E1

E2

12
3
4
5
10

1.00000000000
1.25000000000
1.36111111111
1.42361111111
1.46361111111
1.54976773117

1.50000000000
1.58333333333
1.61111111111
1.62361111111
1.63976773117

1.62500000000
1.63888888889
1.64236111111
1.64470600277

νE3(ν−3)(ν−4)

E4(ν−5)E5

開 [12] と呼ばれる。幾何学的漸近展開をもつ数列は、

4 1.64351851852
5 1.64467592593 1.64483855434
10 1.64493105421 1.64493397472 1.64493406098

エイトケン ∆2 法、リチャードソン補外あるいは ² 算

ν法などで加速できる。

10 1.64493406611 1.64493406667 1.64493406677の形の漸近展開を持つとき、(6) は幾何学的漸近展(3) あるいは (4) を満たす数列の部分列 {s2ν } は(ν−6)

E6(ν−7)E7(ν−8)E8(ν−9)νE9

10 1.644934066785491

s2ν ∼ s + c1 (2−1 )ν + c2 (2−2 )ν + c3 (2−3 )ν + · · ·(1)s2ν ∼ s + c1 (4−1 )ν + c2 (4−2 )ν + c3 (4−3 )ν + · · ·10 項までの部分和から E9 = 1.644934066785491が得られ、π 2 /6 と 10.42 桁一致している。

なる漸近展開を持つ。円周率の計算で関孝和がエイ
トケン ∆2 法を用い、建部賢弘がリチャードソン補外を用いそれぞれ加速に成功した理由は、部分列 {s2ν }数列変換の利用

8が幾何学的漸近展開を持つことによる。

対数収束数列の加速に適した数列変換は各種開発

複合台形公式の誤差は分割数の逆数のオーダーでされている。ここでは、レヴィン u 変換とウィンのρ 算法を取り上げる。

あるので、分割数を 2 倍 2 倍と増やしていくと幾何学的収束になる。ロンベルク積分法はこのことに基づきリチャードソン補外を適用したものである。

8.17レヴィン u 変換

レヴィン変換は、E 算法において漸近展開の利用

gj (ν + l) = Rν+l (ν + l)1−j対数収束数列の漸近列と係数が両方とも既知のときは、5 節の補正による加速が適用できる。(漸近列としたものである (2008 年 11 月号)。ここで、Rν =については 2008 年 6 月号を見よ。) 漸近列が既知でν∆sν−1 = ν(sν − sν−1 ) とおくのが u 変換、あれば、係数が未知であっても前回話した E 算法にRν =より加速することができる。

3∆sν ∆sν−1(sν+1 − sν )(sν − sν−1 )=∆2 sν−1sν+1 − 2sν + sν−1とおくのが v 変換、Rν = ∆sν = sν+1 − sν とおくある。数列 {sν } に対する ρ 算法はのが t 変換である。
∑ν2sν =i=1 1/i にレヴィン u 変換を適用する。

ρ−1 = 0(ν)(ν)

Rν = 1/ν 2 なので、j = 1, 2 . . . に対し

ρ0gj (ν + l) = (ν + l)(ν + l)−2 (ν + l)1−j = (ν + l)−jρk = ρk−2 += sν(ν)

k(ν+1)(ν+1)(ν)ρk−1 − ρk−1,(8)k = 1, 2, . . . ,である。(j = 1, 2, 3 の場合は、(7) に一致するので(ν)

Tk , k = 1, 2, 3 は表 3 に一致する。なお、このよう∑νなことが生じるのは sν = i=1 1/i2 のみである。)
表 4:

∑ν

ようになる。

2i=1 1/i にレヴィン u 変換を適用
(ν−1)

(ν)

νT0

T1

123
4
5
10

1.00000000000
1.25000000000
1.36111111111
1.42361111111
1.46361111111
1.54976773117

1.50000000000
1.58333333333
1.61111111111
1.62361111111
1.63976773117

ν

T3

T4

4 1.64351851852
5 1.64467592593
10 1.64493105421

1.64496527778
1.64493499097

(ν−3)

(ν−4)

(ν−6)

ν

により定義される。
∑νρ 算法を sν = i=1 1/i2 に適用した結果は表 5 の

T6

(ν−7)

T7(ν−2)T2

表 5:

1.62500000000
1.63888888889
1.64236111111
1.64470600277
(ν−5)

T5

2
i=1 1/i に ρ 算法を適用

(ν)

(ν−2)

ν

ρ0

12
34
510

1.00000000000
1.25000000000
1.36111111111
1.42361111111
1.46361111111
1.54976773117

νρ6(ν−4)

ρ2

(ν−6)

ρ4

1.65000000000
1.64682539683
1.64583333333 1.64489489489
1.64503088906 1.64493357289
(ν−8)

ρ8

10 1.64493407610 1.644934066277279

1.64493417081

(2)

(ν−8)

ρ8 = 1.644934066277279 は π 2 /6 と 9.46 桁一致しており、レヴィン u 変換と同程度である。この例∑νに見られるように、ρ 算法は i=1 1/i2 などに効率よく働く加速法であるが、あまり利用されてない。

T810 1.64493405358 1.644934057048 1.64493406373
(ν−9)

ν∑ν

T9

10 1.644934066239456

最大の理由は発表以来 30 数年にわたり収束定理がな= 1.644934066239456
かったことによると思われる。あるいは規則性が分が得られ、π 2 /6 と 9.43 桁一致している。
からなかったということかもしれない。ρ 算法の収(1)

10 項までの部分和から T9束定理は筆者により与えられた。

8.2

ρ 算法

定理 1 (長田 [6]) 漸近展開()
c1c2sν ∼ s + ν −1 c0 ++ 2 + ···νν

D.A. スミスと W.F. フォードが加速法の比較を行ったことは前回 (2008 年 11 月号) 話した。彼らはを持つ数列 {sν } に ρ 算法を適用したとき、8 種類の対数収束級数を取り上げた。そのうち 6 個の級数に対してはウィンの ρ 算法が最良の結果を与えたが、残り 2 つに対しては ρ 算法は加速に失敗しρ2k = s + O((ν + k)−1−2k ),(ν)

た。8 種すべてに対してはレヴィン u 変換が最良でが成立する。

あった。

証明 略

ν→∞¤

この結果についてスミスとフォードは、「ρ 算法は定理 2 (長田 [7]) 漸近展開()
c1c2sν ∼ s + ν θ c0 ++ 2 + ··· ,νν

8 つの問題のうち 6 つの問題では、u 変換や θ 算法よりよい実行結果である。(2 つで失敗したという) 信頼性のなさと相殺してもわずかな利点である。」とRe θ < 0 (9)を持つ数列 {sν } に ρ 算法を適用したとき、した。

(ν)ρ 算法は P. ウィン [13] が ² 算法と同時期に提案しρ2k − s=0
ν→∞ sν+2k − slim

た収束の加速法である。² 算法に類似の菱形算法で4となるための必要十分条件は、θ が負の整数のとき漸近展開した際の、漸近列と係数の両方のデータをである。

用いたオイラーの計算 (補正) が最も精度が良く、漸近列のデータを用いた E 算法がそれに続く。s1 , . . . , s10¤ の値のみを用いた ρ 算法、レヴィン u 変換は 1 桁程スミスとフォードの問題で ρ 算法が失敗した 2 題度劣る。といっても、スミスとフォードの比較テスの問題はいずれも (9) の θ が整数ではなかった。収トに取り上げられた他の 3 つの算法に比べれば、格束定理が得られた後も ρ 算法の利用は増えてないよ段に優れている。
うである。
対数収束数列に対しては万能の加速法は存在しな
http://www.cis.twcu.ac.jp/~osada/rikei/ いので、与えられた数列の漸近展開を考え、適切なに ρ 算法の C 言語によるプログラムをおいておく。
加速法を選ぶのが良い。対数収束数列には、漸近列
証明 略に (log ν)p /ν q , (p, q ∈ N) などの項を持つものもある。これらに対してはレヴィン u 変換や ρ 算法は無加速不可能性

9力であるが、レヴィン変換を拡張した変換 GREP で2加速できることがある。詳細は [9] を見よ。

線型収束する数列はエイトケン ∆ 法で加速することができる。対数収束する数列を１つの数列変換で加速することはできないだろうか。これについて11は、J.P. ドゥラエと B. ジェルマンボンによって否定的に解決されている。これらのことは１冊のテキ加速法を終えるにあたって前編、中編、中編の２、後編と 4 回にわたった加スト [2] にまとめられている。

速法の話はこれで終える。数値計算において加速法
J.P. ドゥラエは数列の集合に残存性という概念を を選択するポイントを挙げておく。
導入し、残存性をもつ数列の集合に属す数列全体を1. 反復が 1 次収束より高次数のときは、加速法は１つのアルゴリズムで加速することはできないとい必要ない。

うことを証明した [1]。

2. 線型収束数列 {sν } が漸近展開

そのことから次の定理が導かれる。
定理 3 (J.P. ドゥラエ・B. ジェルマンボン [3]) 対数sν = s +収束する数列全体の集合 LOG は１つのアルゴリズ

λj 6= 1(j = 1, 2, . . .) (10)を満たすとき、証明 日本語で書かれたものに [11] がある。

¤(a) c1 , c2 , . . . ; λ1 , λ2 , . . . が既知の場合、補正による加速。
(b) λ1 , λ2 , . . . が既知の場合、リチャードソン

対数収束数列の加速 — まとめ補外による加速。

今回はバーゼル問題の級数を巡る加速法の話をし(c) λ1 , λ2 , . . . が未知の場合、² 算法による加速。

た。10 項までの部分和からそれぞれの加速法で π 2 /6

3. 線型収束であるが漸近展開が未知の場合は、エイトケン ∆2 法、レヴィン u 変換による加速。と何桁一致するかをまとめておく。補正による方法
以外の計算は倍精度 (10 進 16 桁弱) によっている。
表 6: 級数cj λνj ,j=1ムで加速することはできない。

10∞∑

4. 交代級数の場合は、エイトケン ∆2 法、レヴィ∑ν2i=1 1/i に対する加速法の比較

ン u, t 変換による加速。

方法

データ一致桁数

オイラーの計算 (補正)表217.41

漸近列による E 算法

表3sν ∼ s +を満たすとき、

レヴィン u 変換

表410.429.43

ρ 算法

表59.465. 対数収束数列 {sν } が漸近展開∞∑cj gj (ν),

j=15limν→∞g1 (ν + 1)
=1
g1 (ν)

(a) cj と gj (ν) がともに既知の場合、補正によ [9] A. Sidi, Practical Extrapolation methods, Cambridge, 2003.
る加速。
(b) gj (ν) が既知の場合、一般化リチャードソ [10] D.A. Smith and W.F. Ford, Acceleration of
ン補外 (E 算法) による加速。
linear and logarithmic convergence, SIAM J.
−j
Numer. Anal. 16(1979), 223–240.
• gj (ν) = ν のときは ρ 算法、レヴィン u 変換による加速。

[11] 杉原正顯・室田一雄、数値計算法の数理、岩波書店、1994• gj (ν) = ν −mj , m ∈ N のときは、部分列 {s2ν } を加速。

[12] G. Walz, Asymptotics and Extrapolation,
Akademie Verlag, 1966(c) gj (ν) が未知の場合は難しい。

[13] P. Wynn, On a Procrustean technique for the

参考文献

numerical transformation of slowly convergent
sqeuences and series, Proc. Camb. Phil. Soc.,[1] J.P. Delahaye, Algorithmes pour suites non con-52(1956), 663–671.

vergentes, Numer. Math. 34(1980), 333–347.
[2] J.P. Delahaye,Springer, 1988

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

Transformation,[3] J.P. Delahaye and B. Germain-Bonne, The setof logarithmically convergent sequences cannot
be acceleratrd, SIAM J. Numer. Anal. 19(1982),840–844.
[4] L. Euler, De summis serierum reciprocarum (逆数の級数の和について), Commentarii academiae

scientiarum Petropolitanae 7(1740), 123-134
http://www.math.dartmouth.edu/~euler/pages/E41.html[5] L. Euler, Institutiones calculi diﬀerentialis cumeius usu in analysi ﬁnitorum ac doctrina serierum(微分学教程), 1755
http://www.math.dartmouth.edu/~euler/
pages/E212.html[6] N. Osada, A convergence acceleration method
for some logarithmically convergent sequences,
SIAM J. Numer. Anal 27(1990), 178–189.
[7] N. Osada, An acceleration theorem for the ρalgorithm, Numer. Math. 73(1996), 521–531.
[8] D.J. Pengelley, Dances between continuous anddiscrete: Euler’s summation formula, in R.E.
Bradley et al. Euler at 300, MAA, 2007.

6