数学公式（３） 調和級数の公式

１）　　　　　　Σ1/n　（自然数の逆数のベキ和）の続き

前ページの公式（２）では、Σ(1/n) にはたどり着けませんでした。簡単な式なのですが、すべての項目となると見つけるのが大変です。事実、これを発見したのは公式（２）を見つけたあとのことでした。

というわけで、Σ(1/n)　については、もう少しお待ちいただいて、その前提となる別の公式を先に説明させていただきます。Σ1/n については、このページの（３）で説明してあります。


２）　　　　　　Σ(1/n^r)

Σ1/n^r = 1+1/2＾ｒ+1/3^r+1/4^r+1/5^r.... 　と並ぶ単純な級数です。（ゼータ関数と呼ばれているとのことです。そんなことも知りませんでした。） 　
　　(TT)　注　　Σ1/n^r　はゼータ関数ではないそうですよ。間違った説明をして申し訳ありませんでした。しかし、納得いかないところがあるので、訂正はしません。ただ、一般的には　ζ(r) = Σ_[n=1,n→∞]1/n^r 　がゼータ関数だとのことですので、ご注意下さい。「大した違いはない」とも思えるのですが、たしかに細部の議論になると違いが出てきますね。数学はネチコイですね。　 　
　これは、r=1 のとき、Σ(1/n) ですから、無限に発散しますが、r>1 の時、ある定数に収束します。ここではおもに　r>1 の時を想定していますが、ここで発見する式がΣ(1/n) に大いに関連することは容易に予測できます。 　
　
まず、すでに判っていることを確認します。
Σ1/n^2 = 1+1/2^2+1/3^2+1/4^2.... ≒ 1.644934067=(pi)^2/6

定数はpiで表せるのですね。不思議な現象です。
それはそれとして、この定数を　ｋとします。すると、Σ1/n^2 = k*f(n) と表記できる　f(n) が存在するはずです。このf(n)を求めてみましょう。

f(n) = 1+a/n+b/n^2.... として、電卓で　a を求めると、 a=0.60745.. と答えが出ました。これは　1/1.644934067 ですから、f(n)　の設定は変更したほうがよさそうです。
それゆえ、Σ(1/n^2)=k+f(n) と問題を変更します。
Σ1/n^2 = k-1/n+1/2/n^2-1/6/n^3+0+1/30/n^5+0+... となります。

これを一般化すると、
Σ1/n^r = k(r)-1/(r-1)/n^(r-1) + 1/2/n^r - r/12/n^(r+1) + 0 + r*(r+1)*(r+2)/6!/n^(r+3) + 0 + ..... となります。

ここで 0 が交互に現れるところに注目します。また、Σ(1/n^2) とは、Σ(n^r) の変形であり、　r<0　の場合ですから、Σ(n^r) の結果との類似性が生じることは充分に考えられます。ですから、Σ(1/n^r) は su(r) を使って表現できる可能性があることはすぐに思いつきます。しかし、符号が異なるところがあるのはやっかいです。

Σ1/n^ｒ=k（ｒ）-1/(r-1)!/n^(r-1)*{su(0)*(r-2)!-su(1)*(r-1)!/n+su(2)*r!/n^2+su(4)*(r+2)!/n^4+su(6)*(r+4)!/n^6....} となります。

別の表現も可能です。
Σ1/n^ｒ = k（ｒ） + 1/(r-1)/n^(r-1) * { -1 + su(1)*(r-1)/n - su(2)*(r-1)*r/n^2 - su(4)*(r-1)*r*(r+1)*(r+2)/n^4 - su(6)*(r-1)*r*...*(r+4)/n^6 .... } となります。

いずれにせよ、　n^1 の項目符号が全体と一致しません。これは　su(r) の持つ、特殊な性質から来る現象ですので、ここではこのまま受け入れておくことにします。符号が不自然な点は注意が必要です。



３）　　　　　　再び Σ(1/n) 　（調和級数の結論）　

Σ1/n^r でsu(r) が使われたことは、大きな示唆を与えてくれます。Σ(n^r) で、 r<0 の時が２）の式になるのですから、 r=-1 という特殊なケースでも、何らかの形で su(r) が関係している可能性はあります。

そこで、今までのΣ1/n = log(n)/log(e) + 0.577215664... + 1/2n - 1/12n^2+1/120n^4 - 1/280n^6 + 1/504n^8... という推定を見直して、su(r) で表現してみることにします。

Σ1/n^r の時と同様に、符号については混乱していますが、su(r) で置き換えられる可能性はあります。
つまり、
Σ1/n = log(n)/log(e) + 0.577215664... + 1/2n - su(2)/n^2 - 0 - su(4)*3!/n^4 - 0 - su(6)*5!/n^6 - 0 - su(8)*7!/n^8...

su(r) の表は別ページにありますが、参考までにここにも書いておきましょう。
su(0)=1
su(1)=1/2
su(2)=1/12
su(3)=0 以下、r が奇数の時は =0 となります。
su(4)=-1/720
su(6)=1/(6!*42)
su(8)=-1/(8!*30)
su(10)=5/(11!*6)
su(12)=-691/15!
su(14)=7/(14!*6)
su(16)=-3617/(17!*30)

su(r) を使わないで表現すると、
Σ1/n = log(n)/log(e) + 0.577215664... + 1/2n - 1/12n^2 + 1/120n^4 - 1/252n^6 + 1/250n^8 - 5/660n^10 +... となります。　（ベルヌーイ数表示については別ページを参照してください。）

これが正しいかどうかの検証は、厳密ではありませんが、数値を入れて計算してみれば、だいたいの所は判ります。
たとえば、n=2 として、上式のΣ(1/n) に直接入れると、　1+1/2=1.5 となります。一方で、公式の方に n=2 を代入しますと、極限値を計算しなければなりませんが、f(n)=1.499999723 となりますので、正しそうです。
n=10 の時、2.928968252 ですが、公式では、 2.928968262 あたりを振動します。正しいことにしておきましょう。

＜追加＞
以上の結果は推定に基づいていたので、以下のようなやり方で上記の式が正しいことを証明してみました。
Σ(1/n) + 1/(n+1) = Σ(1/(n+1)) 　という等式を利用して、この両辺に上記の式を当てはめてみたところ、左辺＝右辺という結果を得ることが出来ました。これは意外と簡単でした。

少し詳しく説明すると、左辺Σ(1/n)には上記の式を代入し、1/(n+1) には　1/n-1/n^2+1/n^3-1/n^4+ . . . .を代入します。右辺 Σ(1/(n+1)) には、log(n+1)/log(e) + 0l577215664 ... + 1/2(n+1) - 1/12(n+1)^2 + 1/120(n+1)^4 - ..... を代入します。log(n+1) = log(n) + log(e) * (1/n - 1/2n^2 + 1/3n^3 - 1/4n^4 .... ) ですから、両辺からlog(n)を消すことが出来、また、0.577...も共通ですから消すことが出来ます。

残りを整理すると、左辺は　1/2n -1/12n^2 +1/120n^4 -1/252n^6 +... +1/n-1/n^2+1/n^3-1/n^4+ .... となります。右辺は、 1/n-1/2n^2+1/3n^3- .... +1/2 * {1/n-1/n^2+1/n^3-....} -1/12 * {1/n^2-2/n^3-3/n^4+.... } + 1/120 * {1/n^4-4/n^5+5*4/2n^6 -.... } - ..... となります。これを 1/n , 1/n^2, 1/n^3 と順番に整理してゆくと、すべての項目で左辺＝右辺を導くことが出来ます。「すべて」と言っても、「計算した範囲では」と言うことですが、公式の証明としてはこれで充分です。

なお、「岩波・数学公式�U、p7」にΣ(1/n) の解説が載っていました。それによると、Σ1/n = g+logn+1/2n + 1/2 - a2/n(n+1) - a3/n(n+1)(n+2) . . . ということで、a2=1/12, a3=1/12, a4=19/120, a5=9/20, となっていました。これには「近似値である」との但し書きが付いていましたが、それにしても随分大雑把な近似値です。計算してみたところ、1/2n - 1/12/n^2 + 1/120/n^4 までしか特定されていませんでした。上記の式とどちらが優れているかは一目瞭然です。かなり自信が湧いてきました。

＜追加　その２＞
ベルヌーイ数は振動発散するので数式での厳密な一致は難しいのですが、多桁計算ソフトを使って、以下のような形で、かなりの精度の一致を確認することが出来ました。

n=10 の時、上記の式を使って、左辺は普通に加算した結果 sahen、右辺は式にもとづいてベルヌーイ数１０まで計算した結果 uhen を算出してみました。すると、sahen/uhen = 1.000000000..... と０が１４個並ぶほどの精度で一致していました。ベルヌーイ数を増やしてみると、B(58)で０が２８個並び、B(66)で２８個並んでいるので、このあたりが精度の限界のようです。
n=20 の時、同様に計算すると、ベルヌーイ数１0で、sahen/uhen = 1.0000000..... と０が１０個並びました。ベルヌーイ数を増やしてみると、B(130)あたりが精度の限界のようで、０が５５個並びました。
ベルヌーイ数を増やしても精度は限界がありますが、ｎが増えればいくらでも精度は上がります。たとえば、n=100で、B(130)までの精度で計算したところ、０が１４７個並んでいました。

＜追加　その３＞
「調和級数」 (harmonic series) という名前がどうして付けられたか、気になっていたのですが、「数学おもしろ講座」のHPでその説明を見つけました。「調和数列」というのがあって、それは各項の逆数が等差数列になっているものだそうです。とすると、１，２，３，・・・だけでなく、１５，５，３，１５／７，５／３．．．．なども調和数列になります。これらが調和数列と呼ばれる由来は、音階がこのような数列になっているからだそうです。ドの音を出す弦の長さを１とすると、レは８／９、ミは６４／８１、ファは３／４、ソは２／３、ラは１６／２７、シは１２８／２４３、ドは１／２となっているそうです。全音を２，半音を１として数えると、ドからソまでが７，ソから上のドまでが７ですから、これをならべると、1, 2/3, 1/2, 2/5, 1/3 .... となります。この逆数は等差数列になっています。というわけで、和音 (harmony) から作られた数列ということで、調和数列 （harmonic progression）と命名されたとのことです。とても面白いですね。


４）　　　　Σn^r　　（　r が奇数、および実数の場合）

Σn^r は公式（１）で明らかにしたつもりですが、ｒが偶数でないときにはやや問題が生じます。そのあたりを解決しておきたいと思います。

たとえば、r=2 の時、公式（１）に代入して、Σn^2 = 1/3*n^3 (1+3/2/n+3*2/12/n^2) を得ます。これを整理すると、n*(n+1)*(2n+1)/6 となり、目的の式を得ることが出来ます。しかし、r=3にすると、Σn^3 = 1/4*n^4 ( 1 + 4/2/n + 4*3/12/n^2 - 4*3*2*1/6!/n^4 ) となり、最後の項目が余計になります。これを整理すると-1/5!という綺麗な数字になるのですが、これでは式が成り立ちません。

また、r=1.5 だと、Σn^1.5 = 1/2.5*n^2.5* ( 1 + 2.5/2/n + 2.5*1.5/12/n^2 - 2.5*1.5*0.5*(-0.5)/720/n^4 +....) のような無限級数になり、計算不能となります。しかも、問題なのは、不能だけでなく、どこか間違っているようです。なぜなら、電卓で計算すると、Σn^1.5 = 1/2.5*n^2.5 (1 + 1.25/n + 0.3125/n^2 - 0.0637/n^2.5 +....) というように、n^2.5 の項目が登場するからです。これは公式（１）の中にはありません。

それゆえ、rが偶数の時に消える項目が存在していることは間違いないと判断できます。この項目がどうなっているかを調べてみると、
Σn^r = n^(r+1)/(r+1) * { 1+(r+1)/2/n + (r+1)r/12/n^2 - (r+1)r(r-1)(r-2)/720/n^4 +...... } - r!su(r+1) と確定できるように思えてきました。しかし、su(r)はrが実数のときはまだ定義されていません。

そこで、厳密に考えるために、以前提示した公式を再定義して、su(r)ではなく、se(r) = 2/(2pi)^r* { 1 + 1/2^r + 1/3^r + 1/4^4.... } と定義しておきます。こうすると符号を無視することが出来、rが正の実数の時にもse(r)は確定することになります。しかし、それと上記の r!su(r+1) は別のものです。そもそも r!su(r+1) は今の段階では推定に過ぎません。

少し複雑になっていますが、su(1)=1/2, su(2)=1/12, su(3)=0, su(4)=-1/720... であるのに対し、se(1)=∞, se(2)=1/12, se(3)=0.00967..., se(4)=1/720,... となります。

そこで、se(r) と r!su(r+1) の関係を調べるために、Σn^r = 1/(r+1)*n^(r+1)* { 1+(r+1)/2/n +..... }　-r!se(r+1)f(r)　となるf(r)を定義して、このグラフを作ってみました。つまり、
f(r) = 1/r!se(r+1)* { -Σ(n^r) + 1/(r+1)*n^(r+1)* { 1+(r+1)/2/n + (r+1)r/12/n^2-..... }
という式になります。右辺は振動して確定しにくいのですが、r!はこのページの５）を見ると、実数の時でも計算可能となります。その結果、nが大きいと間違いなくsin曲線になり、nが小さいときでもかなりsin曲線に近いことが判りました。たとえば、f(1)=1.01218となり、f(0.5)=0.7924となりました。f(2)=0,f(3)=-1,f(4)=0,f(5)=1と変化します。

ゆえに、この曲線をsin曲線と見なすと、su(r) = sin(pi*r/2-pi/2)*se(r) という結果を得ました。（ r!se(r+1) である点にご注意ください。）これは、su(r) = -sin(pi*(r+1)/2)*se(r) とも表記することが出来ます。いずれにせよ、rが奇数、正の実数の時も成り立つと言うことです。

この式は、計算上は問題ないのですが、理屈の上では確定というにはもう一歩のところがあります。これについては、これからさらに検討してゆきたいと思います。

再度書いておくと、
Σn^r = n^(r+1)/(r+1) * { 1 + (r+1)/2n + (r+1)r/12n^2 - (r+1)r(r-1)(r-2)/720n^4 + ...... } -r!su(r+1)
ということです。



５）　　　n! 、階乗の計算、階乗の求め方　（　n が実数の場合）

n! は通常正の整数に対して定義されています。つまり、 n!=n*(n-1)*(n-2)*(n-3)*....3*2*1 となります。
しかし、公式（２）を用いると、n が-1より大きな実数の時にも計算可能となります。（n<-1の時も計算できることが判りました。n=-2,-3 などの負整数のときは無限になりますが、それ以外では収束します。）

n!=(n/e)^n*√(2*pi*n)*e^(1/12n-1/360n^3+1/1260n^5....)
つまり、
n!=(n/e)^n*√(2*pi*n)*e^｛su(2)/n+su(4)*2!/n^3+su(6)*4!/n^5+su(8)*6!/n^7....｝
です。

この式の su(2)/n+su(4)*2!/n^3+su(6)*4!/n^5+su(8)*6!/n^7.... 部分は振動し、nが小さいとそれだけ確定しにくくなります。

そこで、(n-a)!=n!/n/(n-1)/(n-2)/..../(n-a+1) を利用することにします。

たとえば、0.25! を求めてみましょう。公式（２）に代入すると、振動部分が拡大して計算できません。そこで、n-a=0.25と置いて、aは適当に都合の良い大きさの数にします。たとえばa=3としてみましょう。すると、n-3=0.25ですから、n=3.25 となります。

一方、(n-3)!=n!/n/(n-1)/(n-2) ですから、(n-3)!=(n/e)^n*√(2*pi*n)*e^｛su(2)/n+su(4)*2!/n^3+su(6)*4!/n^5....｝/n/(n-1)/(n-2) ですから、n=3.25として右辺を計算します。すると、振動部分がある程度確定しますので、それにより、0.25!=0.9064026024...を得ます。しかし、これでは、小数点５桁以下は信用できませんので、さらに厳密な答えを知りたいときは、a=10とか、a=20として計算すればよいのです。すると、計算は面倒になりますが、0.25!=0.9064024774...を得ます。

別の例をやってみましょう。n=-0.6とします。マイナスの数は、そのままではルートがありますので、実数内での計算が出来なくなります。しかし、a=3 として、n-3=-0.6と置きます。すると、n=2.4となり、これを(n-3)!の式に代入すると、実数内での計算となります。(-0.6)!=2.218162014... もっと精度を上げたいときは、aを大きくすればよいだけです。n=3.4にしてみると、(-0.6)!=2.218159793... となります。