この記事のコードをJupyter notebook形式で見るにはこちら。
せっかく作ったオイラーマクローリン総和公式パッケージを使って、具体的な応用を一つ。同じオイラーが考えたオイラーマスケローニ定数を数値的に求めてみます。（Maximaには%gammaという定数として実装されていますから、テスト的な意味合いしかありません。逆に答え合わせは簡単です）。 ちなみにこの定数の定義式は、
$$\notag\lim_{M\rightarrow\infty}\sum_{n=1}^{M}{\frac{1}{n}}-\log M$$
です。この調和級数から対数を引いたものが収束するのですが、それが求めるオイラーマスケローニ定数です。
早速、すでにダウンロードしてあるeuler-maclaurin-sumパッケージを読み込みます。
(%i1) asdf_load_source("euler-maclaurin-sum");
$$ \tag{%o1} \#<LOAD-SOURCE-OP > $$
オイラーマクローリン総和公式を表示してみます。
(%i2) ems;
$$ \tag{%o2} \sum_{n=N}^{M}{f\left(n\right)}=\sum_{k=1}^{K-1}{\frac{\left(-1\right)^{k+1}\,B_{k+1}\,\left(\left.\frac{d^{k}}{d\,x^{k}}\,f\left(x\right)\right|_{x=M}-\left.\frac{d^{k}}{d\,x^{k}}\,f\left(x\right)\right|_{x=N}\right)}{\left(k+1\right)!}}+\frac{\left(-1\right)^{K+1}\,\int_{N}^{M}{\overline{B}_{K}\left(x\right)\,\left(\frac{d^{K}}{d\,x^{K}}\,f\left(x\right)\right)\;dx}}{K!}+\int_{N}^{M}{f\left(x\right)\;dx}+\frac{f\left(N\right)}{2}+\frac{f\left(M\right)}{2} $$
Mは1以上として、
(%i3) assume(M>1);
$$ \tag{%o3} \left[ M>1 \right] $$
emsにf(x)=1/x, N=1, K=1を代入してみます。
(%i4) ems,f(x):=1/x,N:1,K:1;
$$ \tag{%o4} \sum_{n=1}^{M}{\frac{1}{n}}=\int_{1}^{M}{\frac{d}{d\,x}\,\frac{1}{x}\,\overline{B}_{1}\left(x\right)\;dx}+\int_{1}^{M}{\frac{1}{x}\;dx}+\frac{1}{2\,M}+\frac{1}{2} $$
(%i5) %,nouns;
$$ \tag{%o5} \sum_{n=1}^{M}{\frac{1}{n}}=-\int_{1}^{M}{\frac{\overline{B}_{1}\left(x\right)}{x^2}\;dx}+\log M+\frac{1}{2\,M}+\frac{1}{2} $$
左辺は調和級数の和で、右辺にはちょうどlog(M)があります。これを移行すれば定義式に近づきます。
(%i6) EuGamma:%-log(M);
$$ \tag{%o6} \sum_{n=1}^{M}{\frac{1}{n}}-\log M=-\int_{1}^{M}{\frac{\overline{B}_{1}\left(x\right)}{x^2}\;dx}+\frac{1}{2\,M}+\frac{1}{2} $$
右辺の最初の積分はM→∞の時広義積分になりますが、グラフで様子を確認します。区間[1,5]での分子の1次の周期的ベルヌーイ多項式、分母の1/x^2、両者をかけた被積分関数をプロットしてみます。
(%i7) set_plot_option([svg_file, "maxplot.svg"])$
(%i8) plot2d([bernpoly(x-floor(x),1),bernpoly(x-floor(x),1)/x^2,1/x^2],[x,1,5])$
1次の周期的ベルヌーイ多項式は鋸の歯のグラフです。これに薄緑の2次の双曲線をかけると、赤い、鋸の歯がだんだん潰れるグラフになります。収束しそうですね。
この積分を変形するとこんな感じになります。
(%i9) apply1(first(rhs(EuGamma)), emsRemInt);
$$ \tag{%o9} -\sum_{k=1}^{M-1}{\left(\log \left(k+1\right)-\log k+\frac{1}{2\,k+2}+\frac{k}{k+1}-\frac{1}{2\,k}-1\right)} $$
100項くらいまでの和で広義積分の近似とすると、積分の値は、
(%i10) rem:%,M:100,nouns,numer;
$$ \tag{%o10} 0.0772073316515306 $$
となります。
それ以外の項はM→∞を計算し、上記の近似値と足し合わせます。
(%i11) limit(rest(rhs(EuGamma)),M,inf)+rem,numer;
$$ \tag{%o11} 0.5772073316515306 $$
さて、答え合わせをしてみます。%gammaを数値で表示させます。
(%i12) %gamma,numer;
$$ \tag{%o12} 0.5772156649015329 $$
うーむ、４桁しか合っていません。元々収束遅そうだし、適当にやったので仕方ないのですが、、、。
というわけで加速的なやり方に変更です。積分の部分が十分に小さくなるように、N=10000より先の部分をこの公式で求め、N=9999までは普通に足算することにします。また被積分関数の分母の次数はK+1ですからK=6と大きく取ります。
(%i13) ems,f(x):=1/x,N:10000,K:6,nouns;
$$ \tag{*} \verb|Is |M-9999\verb| positive, negative or zero?| $$
pos;
$$ \tag{%o13} \sum_{n=10000}^{M}{\frac{1}{n}}=-\int_{10000}^{M}{\frac{\overline{B}_{6}\left(x\right)}{x^7}\;dx}+\log M+\frac{1}{2\,M}+\frac{\frac{1}{100000000}-\frac{1}{M^2}}{12}-\frac{\frac{3}{5000000000000000}-\frac{6}{M^4}}{720}+\frac{\frac{3}{25000000000000000000000}-\frac{120}{M^6}}{30240}-\log 10000+\frac{1}{20000} $$
右辺のlog(M)を左辺に移行してオイラーマスケローニ定数の定義に合わせます。
(%i14) EuGamma:%-log(M);
$$ \tag{%o14} \sum_{n=10000}^{M}{\frac{1}{n}}-\log M=-\int_{10000}^{M}{\frac{\overline{B}_{6}\left(x\right)}{x^7}\;dx}+\frac{1}{2\,M}+\frac{\frac{1}{100000000}-\frac{1}{M^2}}{12}-\frac{\frac{3}{5000000000000000}-\frac{6}{M^4}}{720}+\frac{\frac{3}{25000000000000000000000}-\frac{120}{M^6}}{30240}-\log 10000+\frac{1}{20000} $$
積分の部分の評価をしてみます。分子は6次の周期的ベルヌーイ多項式なので、その最大最小はベルヌーイ多項式の区間[0,1]での最大最小と一致します。
(%i15) diff(bernpoly(x,6),x),factor;
$$ \tag{%o15} \left(x-1\right)\,x\,\left(2\,x-1\right)\,\left(3\,x^2-3\,x-1\right) $$
最後の2次式は[0,1]で0になりません。従ってこの区間で極値は0,1/2,1でとることが分かります。そこでの値を計算すると、
(%i16) [bernpoly(0,6),bernpoly(1/2,6),bernpoly(1,6)];
$$ \tag{%o16} \left[ \frac{1}{42} , -\frac{31}{1344} , \frac{1}{42} \right] $$
というわけで絶対値の最大は1/42です。そこで周期的ベルヌーイ多項式を定数1/42で置き換えて広義積分を計算します。
(%i17) [err:integrate(%[1]/x^7,x,10000,inf),ev(err,numer)];
$$ \tag{%o17} \left[ \frac{1}{252000000000000000000000000} , 3.968253968253969 \times 10^{-27} \right] $$
という訳で、この積分の寄与は3.9683*10^(-27)で抑えられることがわかります。無視しても26桁は求められる訳です。
多倍長浮動小数の桁数を30桁に設定して、1から9999までの調和級数の和と、(%o14)の右辺の、積分以外の項をM→∞として足し合わせます。
(%i18) fpprec:30;
$$ \tag{%o18} 30 $$
(%i19) bfloat(sum(1/n,n,1,9999)+limit(rest(rhs(EuGamma)),M,inf));
$$ \tag{%o19} 5.77215664901532860606512090082B-1 $$
では答え合わせです。こちらも多倍長浮動小数で求めます。
(%i20) bfloat(%gamma);
$$ \tag{%o20} 5.77215664901532860606512090082B-1 $$
30桁全て一致しました！
想定よりも高精度に計算できているのは、広義積分の評価で6次の周期的ベルヌーイ多項式を1/42で抑えてから積分していますが、ここの部分が雑な評価だったということですね。
