数式処理ソフト DERIVE(デライブ) de ドライブ

54.数値積分(3)(Euler-Maclaurin展開とラプラス変換、Wijngaarden変換)

Euler-Maclaurin展開のラプラス変換法による導出の試み

「数式処理ソフト DERIVE(デライブ)の第53回で『Euler-Maclaurin(オイラー-マクローリン)展開』を森口先生の『計算数学夜話』(森口繁一:1978年4月30日:日本評論社)に基づき、いわば、発見法的に導き、その二、三の応用例を紹介したのじゃな。
 なお、同書では、導出の次の章で、より数学的な方法による証明も行われておる」

「『発見法的』か~。なるほど。確かに、ちょっと、だまされたような気もするわね」

「そこで、Euler-Maclaurin(オイラー-マクローリン)展開のオーソドックな証明は、同書をご覧いただくことにして、ここでは、ラプラス変換法による導出を試みてみよう。
とは言っても、それほど難しいところはないので、安心して欲しい。
まずは、関数 f(t)とそのラプラス変換をF(s)とする。
なお、ここでは、ラプラス変換の一般的なテキストにならって、変数をtとしておく。
最初に関数の和、f(t+a)+f(t+a+h)+f(t+a+2h)+・・+f(t+a+(n-1)h) を考えよう。
ここで、tは、最終的には、ゼロと置くのじゃが、tの関数にしないとラプラス変換が有効に使えないからだ。
また、t<0のとき、f(t)=0、従って、a<bは、いずれも正とし、h=(b-a)/nと置く。
ともちゃん、この関数和のラプラス変換を求めてご覧」
「う~ん。
近藤先生の本(『演算子法』:近藤次郎:1970年6月30日:培風館)にならうと、cを、c>0の定数として、
f(t)⊃F(s)のとき、f(t+c)⊃∫(0~∞)exp(-st)f(t+c)dt、u=t+cとおくと、t=u-cであるけど、
元の積分範囲は、t=0~∞だから、うっかりすると、uの範囲は、c~∞かと思うと、u=0~cの積分はどうなるか。
実は、この辺、ちょっと違うのよね」
「そう、ここは、うっかりしてしまう。
t<0では、f(t)=0としているので、uの0~cまでは、f(u)は、ゼロとなるのじゃ。
なんとなれば、仮定では、t<0では、f(t)=0なのだから、u=0では、t=-c、u=cで、ようやく、t=0だからじゃ。
従って、f(t+c)⊃∫(0~∞)exp(-s(u-c))f(u)du=exp(cs)∫(0~∞)exp(-su)f(u)du=exp(cs)×F(s) となる」
「よって、関数和⊃Σ(k=0~n-1)exp((a+kh)s)F(s)=exp(as)(1-exp(nhs))/(1-exp(hs))×F(s)
=(exp(as)-exp(bs))/(1-exp(hs))×F(s)となるわ。ただし、a+nhs=bsを使う」
「そう。ここで、ベルヌイ数の母関数、x/(exp(x)-1)の展開が、Σ(k=0~∞)B(k)xk/k! を思い出す必要がある」
「ていうか、絶対、思い出すよ。xをhsと考えるんだけど、
Σ(k=0~n-1)f(a+kh)⊃(exp(as)-exp(bs))×(1/hs)(Σ(k=0~∞)B(k)(hs)k/k! )F(s)
ところが、(1/hs)Σ(k=0~∞)B(k)(hs)k/k!=1/hs+Σ(k=1~∞)B(k)(hs)k-1/k! であるので、
Σ(k=0~n-1)f(a+kh)⊃(exp((bs)-exp(as))F(s)/hs+(exp(bs)-exp(as))(Σ(k=1~∞)B(k)(hs)k-1/k! )F(s) 」
「これを逆変換するために、(exp(bs)-exp(as))F(s)/hs をまずは、考えよう。
F(s)/s⊂∫(τ=0~t)f(τ)dτであることを使うと、(exp(bs)-exp(as))F(s)/hs ⊂(1/h)(∫(τ=0~t+b)f(τ)dτ-∫(τ=0~t+a)f(τ)dτ)
=(1/h)∫(τ=t+a~t+b)f(τ)dτ となる。これは、t→0で、(1/h)∫(τ=a~b)f(τ)dτとなる。
次は、Σの部分じゃ」
「Σのところは、(exp(bs)-exp(as))(Σ(k=1~∞)B(k)(hs)k-1/k! )F(s)で、まず、第1項 k=1の部分は、B(1)=-1/2なので、
-1/2(f(t+b)-f(t+a))となり、これも、t→0で、-1/2(f(b)-f(a))となる。
第2項以下は、Σ(k=2~∞)B(k)h(k-1)/k!×(f(k-1)(t+b)-f(k-1)(t+a))なので、Σ(k=2~∞)B(k)h(k-1)/k!×(f(k-1)(b)-f(k-1)(a))となる。
従って、f(t+a)+f(t+a+h)+f(t+a+2h)+・・+f(t+a+(n-1)h)→f(a)+f(a+h)+f(a+2h)+・・+f(a+(n-1)h)は、
(1/h)∫(τ=a~b)f(τ)dτ-1/2(f(b)-f(a))+Σ(k=2~∞)B(k)h(k-1)/k!×(f(k-1)(b)-f(k-1)(a))となって、
前回の結果と一致している」
「そうじゃな。定積分としては、∫(τ=a~b)f(τ)dτ=h×Σf(a+kh)+h/2(f(b)-f(a))-Σ(k=2~∞)B(k)hk/k!×(f(k-1)(b)-f(k-1)(a))
と移項してh倍した方が見やすい。
なお、Σ(k=2~∞)B(k)hk/k!×(f(k-1)(b)-f(k-1)(a))は、実際は、kが2以上の偶数でのみB(k)は、ゼロでないので、前回の記事のように
Σ(j=1~∞)B(2j)h2j/(2j)!×(f(2j-1)(b)-f(2j-1)(a))とした方が照合しやすいじゃろう。
ラプラス変換では、t>=0のみを考えているんじゃが、f(x)のxが負の範囲で定義されていても、有限範囲で、ゼロとなる関数であれば、
原点を移動すれば、同じ議論が成り立つじゃろう」
「x→-∞まで、f(x)がゼロでないときは、どうするのかな?」
「ラプラス変換の代わりに『フーリエ変換』を使えば、同種の議論が成り立つのではないかと思っているがの。
そのうち、調べて見てご覧」
「付け加えると、マクローリン(Maclaurin)さん(1698年~1746年)人は、どこの人かと思って調べて見たの。
岩波数学辞典では、18世紀のスコットランドの数学者で、あのテイラー級数のTaylerさんと同時代の人と紹介されていたわ。
マクローリンさんは、あのオイラーさん、Leonhard Euler(1707年~1783年)とも、ほぼ、同時代なの。
また、Web(『科学史の小窓』有賀暢迪氏のHP。http://www.ariga-kagakushi.info/)で更に調べると、コリン・マクローリン、Colin Maclaurinは、ニュートンの後継者としてみなされるとのことだけど、46才という若さで亡くなったことがイギリスの科学界において惜しまれたそうよ」

Euler-Maclaurin展開によるζ関数の計算

「Euler-Maclaurin展開は、Euler変換のように交代級数である必要は無いので、たとえば、ツェータ関数、ζ(x)=Σ(k=0~∞)(1/kx)、x>1 のような関数に適用可能ね。
これまでも計算したけど、ζ(1)は、発散するけど、x>1である、ζ(2)=π2/6≒1.644934066などは、収束するのね」

「そうじゃったな。オイラー・マクローリン展開を使ってζ関数の値を計算してみようか。
g(x):=1/x^2として、h=(b-a)/n、ベルヌイ数B(m)をm項まで計算するとして、f(a,b,n,m)を下記のように定義する。
f(a,b,n,m):=∫(g(x), x, a, b) - h(a, b, n)/2(g(b) - g(a)) + IF(m = 0, 0, ∑(h(a, b, n)^(2j)B(2j)/(2j)!(LIM(∂(g(x), x, 2j - 1), x, b) - LIM(∂(g(x), x, 2j - 1), x, a)), j, 1, m))
ζ(2)=π2/6≒1.6449340668482264364・・の場合は、a=1、b=n+1、h=1、n=nとして、
f(1,n+1,n,0)=3/2 - (2n + 3)/(2(n + 1)^2)→3/2=1.5
f(1,n+1,n,1)=5/3 - (6n^2 + 15n + 10)/(6(n + 1)^3)→5/3≒1.6666666666666666666
f(1,n+1,n,2)=49/30 - (30n^4 + 135n^3 + 230n^2 + 175n + 49)/(30(n + 1)^5)≒1.6333333333333333333
f(1,n+1,n,3)=58/35 - (210n^6 + 1365n^5 + 3710n^4 + 5390n^3 + 4403n^2 + 1911n + 348)/(210(n + 1)^7)≒1.6571428571428571428
f(1,n+1,n,4)=341/210 - (210n^8 + 1785n^7 + 6650n^6 + 14175n^5 + 18893n^4 + 16107n^3 + 8573n^2 + 2607n + 341)/(210(n + 1)^9)≒1.6238095238095238095
のように、初項から計算すると、小数点以下1桁程度までがやっとで、真値までは、遠く届かない」
「そうなのよね。
でも、森口先生の『計算数学夜話』に出てくるように、初項から、適当な項数までは、正直に計算して、その後の項を近似すれば、かなりの精度で計算できるんだ。
ここでは、初項から、20項までを正直に計算して、Σ(k=1~20)1/k^2≒1.5961632439130233166 を得る。
一方、この値にf(21,n+21,n,m)を加えて、極限値を取ると、ζ(2)の近似値として、
m=0では、 1.6449160783801434980 :小数点以下4桁まで。
m=1では、 1.6449340749967795704 :小数点以下7桁まで。
m=2では、 1.6449340668350486697 :小数点以下9桁まで。
m=3では、 1.6449340668482681885 :小数点以下13桁まで。
m=4では、 1.6449340668482262217 :小数点以下15桁まで。
と俄然、正しい値が出てくるのよ。ここが不思議~」
「一方、正直に計算してみると、
10万項まで、 1.6449240668982262698
20万項まで、 1.6449290668607264156
30万項まで、 1.6449307335204486525
40万項まで、 1.6449315668513514338
50万項まで、 1.6449320668502264351
60万項まで、 1.6449324001829486579
70万項まで、 1.6449326382778182727
80万項まで、 1.6449328168490076861
90万項まで、 1.6449329557377326090
100万項まで 1.6449330668487264363となって、90万項と100万項の間は、ずいぶんと精度に差が出てくる。
これは、DERIVEでは、nが大きいとき、何らかの近似計算を使っているのじゃなかろうかと思うがの」

Wijngaarden(ヴァインハーデン)変換

「この人は、なんと読むのかな?」

Aad van Wijngaarden氏(1916年11月2日~1987年2月7日)は、オランダの数学者で「ヴァインハーデン」と書き表されるようじゃな。
コンピュータ科学や計算数学で多くの業績を残したそうだ。残念ながら、同氏に関する日本語のページは、ほとんどないようだが。
BASICによる高速ラプラス変換』(細野敏夫:共立出版:1984年4月)のP44にその変換方法と例が載っている。
この変換は、Euler変換が交代級数(正負の項が交代して現れる級数)にのみ適用されるのを補おうという目的で考えられた変換のようじゃ。
よって、すべての項が正(正確には定符号)であるような級数に適用可能じゃ。その処方じゃ。
まずは、S=a(1)+a(2)+・・において、a(k)は、すべて正とする。そして、この級数は、『絶対収束』とする。
なお、『絶対収束』とは、Σ|a(k)|が収束する場合である。このときは、項の順序を変えてもその和は変わらない。
ここで、a(k)から次の規則で導かれる、v(k)を考える。
なお、上記書籍に合わせて、k=1から数える。
v(k)=20a(20k)+21a(21k)+22a(22k)+23a(23k)+・・、k=1~
そして、S'=v(1)-v(2)+v(3)-v(4)+-・・を作ると、これは、交代級数なので、Euler変換で高精度の近似計算が可能というものじゃ」
「じゃ、実際に計算してみるよ。
まずは、a(k)を1/k^2と定義する。k=1、2,・・
次に、v(m,p)を∑(2^(k - 1)a(2^(k - 1)m), k, 1, p)と定義する。ここで、pは、v(k)=a(k)+2^1a(2^1k)+・・+2^(p-1)a(2^(p-1))のときの正整数。
そして、f(n,p)=IF(n < 1, 0, ∑((-1)^(k - 1)v(k, p), k, 1, n)) と定義する。このf(n,p)がS'を定義式により、第n項まで計算したものよ。
以下では、p=40として、f(n,p)をn=10~100まで10刻みで計算してみた。
10: 1.6359243512204824050
20: 1.6425587566580120444
30: 1.6438599517567130501
40: 1.6443246820993519954
50: 1.6445420636505620617
60: 1.6446609174136442383
70: 1.6447329000713462196
80: 1.6447797696666976614
90: 1.6448119816294313153
100:1.6448350667467604507
となって、小数点以下、3桁程度は、真値(ζ(2)=π2/6≒1.6449340668482264364・・)に近づいている」
「なるほど。では、今度は、pを変えて、n=100項まで計算してみよう。
10: 1.6432287825033850753
20: 1.6448334981112985373
30: 1.6448350652163843903
40: 1.6448350667467604507
50: 1.6448350667482564194
60: 1.6448350667482564194
70: 1.6448350667482564194
80: 1.6448350667482564194
90: 1.6448350667482564194
100:1.6448350667482564194
p=40で十分のようじゃな。
確かに、p=40ということは、k=2^39、a(k)≒1.818989403×10^(-12))という小さな値じゃからな」
「では、いよいよ、Euler変換をS'に適用してみるよ。
階差関数 g(n,m,p)を、IF(m = 0, v(n, p), g(n + 1, m - 1, p) - g(n, m - 1, p))、と再帰的に定義するの。
これは、第41回の「積分(1)(単積分、数値積分とEuler変換)」に出てきてるよ。
みなさん、覚えてますか?
また、S'≡S(n,p,c)を、f(n, p) + ∑((-1)^(m - 1)g(n + 1, m - 1, p)/2^m, m, 1, c)、として定義する。
ここで、階差関数g(n,m,p)は、Δmv(n)であり、v(n)の第m階差を表す。
このΔも階差の記号:Δv(k)=v(k+1)-v(k) だったわ。
ただし、vは、k=1~k=2pまででの和として、以下の数値計算では、p=40に固定している。
また、S(n,p,c)では、f(n,p)で、正直にn項計算する。(nは偶数とする)
そして、第2項のΣは、n+1項からのオイラー変換による和を行う。オイラー変換の階差は、c個とする。
以下では、c=10として、n=0~20までを2刻みで計算してみた。
n=0の場合は、初項からEuler変換(第10階差まで)を施していることになる。
0: 1.6444257276110604086
2: 1.6449308074731404738
4: 1.6449339402754930964
6: 1.6449340563000305612
8: 1.6449340654630464128
10:1.6449340665999047493
12:1.6449340667916606853
14:1.6449340668321526774
16:1.6449340668423125830
18:1.6449340668452361203
20:1.6449340668461767780 :小数点以下11桁まで。
真値(ζ(2)=π2/6≒1.6449340668482264364・・)と比較すると、初項から20項目まで、vの和を計算してから、Euler変換を施すと、小数点以下、11桁程度まで、一致する」
「う~ん。微妙じゃな。前節のEuker-Maclaurin展開の方法じゃと、初項から20項まで、こちらはa(k)じゃが、を正直に計算してから、近似すると、小数点以下、15桁まで、出ておるからな。
試しに、s(20,80,16)を計算すると、1.6449340668482263918 となって、小数点以下、15桁まで一致する。
これは、初項から第20項までは、上と同様に正直に計算して、21項目以降を、p=80として、階差もc=16階差まで、計算したものじゃ。
BASICによる高速ラプラス変換』に記載されているように、精度と階差数は、だいたい比例するようなので、c=16で、小数点以下15桁まで求まったのは、偶然ではないようだ」
「上の計算では、公式をそのまま、関数として、表現しているんだけど、以前に注意があったように、∑((-1)^(m - 1)g(n + 1, m - 1, p)/2^m, m, 1, c)は、展開して、すべて、v(k)で表しておくと、計算が速い。
実際、a(k)=1/kxとして、n=20、p=80、c=16として、S(20,80,16)の計算をしておいてから、xを2から10まで変化させて代入して計算すると、
ζ(2)≒ 1.6449340668461870104
ζ(3)≒ 1.2020569031583689331
ζ(4)≒ 1.0823232337101141681
ζ(5)≒ 1.0369277551424178108
ζ(6)≒ 1.0173430619843914074
ζ(7)≒ 1.0083492773810166239
ζ(8)≒ 1.0040773561979302870
ζ(9)≒ 1.0020083928260262379
ζ(10)≒1.0009945751269068449 は、すぐに求められる」
「そうじゃのう。重要な点じゃな。
また、今回のζ(x)については、v(k)の定義式から、S'を計算すると、x>1のとき、
ζ(x)=(1/(1-2(1-x))(1-1/2x+1/3x-1/4x+・・)であることが、同書のP45に示されている。
上式とEuler変換を組み合わせると、同書にあるように、いっそう速い計算が行えるであろう」

最終更新日 2011/7/31