Bernoulli Numbers and Harmonic Series
Kevin Brown
i.e.,
1 1 1 1
H(n) = 1 + --- + --- + --- + ... + ---
2 3 4 n
It's well known that
H(n) -> log(n) + gamma as n -> inf
where gamma = 0.57721... is Euler's constant and log is the natural
logarithm. Interestingly, it appears that the discrepancy between
H(n) and log(n)+gamma can be computed very directly from the Bernoulli
numbers. This relation not only allows us to determine H(n) very
precisely, it can also be used to compute the value of gamma from a
known value of H(n).
For any integer n>1 it appears that
inf B_k
H(n) = log(n) + gamma - SUM ----- (1)
k=1 k n^k
where B_k denotes the kth Bernoulli number, i.e.,
B_0 = 1 B_1 = -1/2 B_2 = 1/6 B_3 = 0 B_4 = -1/30
and so on. The summation _appears_ to rapidly converge, so that just
ten non-zero terms are needed to give
H(256) = 6.12434496281728...
Conversely, given (for example) the fact that H(8) = 2283/840, we can
easily compute
2283 inf B_k
gamma = ---- - 3 log(2) + SUM ----- (2)
840 k=1 k n^k
Again taking just the first ten non-zero terms of the summation, this
formula gives
gamma = 0.577215664901533....
which is correct to 15 decimal places. Of course, the value of log(2)
is easily computed using the formula
/ 1 + 1/3 \ / 1 1 1 \
log(2) = log( --------- ) = 2 ( --- + ----- + ----- + ... )
\ 1 - 1/3 / \ 1*3 3*3^3 5*3^5 /
Anyway, more digits of gamma can be determined simply by using a
larger value of n (e.g., n=16) with the same number of terms in the
Bernoulli summation. However...
When I said the summation _appears_ to converge very rapidly, I was
referring to the fact that the terms initially become extremely small
as k increases. This is because the denominator k n^k increases
rapidly whereas the numerator B_k actually gets smaller...for a while.
But in the long run the values of the Bernoulli numbers increase at a
super-exponential rate, as can be seen from the inequality
2(2k)!
| B_(2k) | > ----------- (3)
(2 PI)^(2k)
Therefore, in the long run the numerators of the terms in the
summation rise faster than the denominators, so the summations in
equations (1) and (2) actually diverge. This strikes me as a very
interesting example of asymptotic divergent series. Has anyone ever
studied this to see how much precision can be achieved with this
series?
Kevin Brown
> Let H(n) denote the sum of the first n terms of the harmonic
>
> inf B_k
> H(n) = log(n) + gamma - SUM ----- (1)
> k=1 k n^k
>
Thanks to Kurt Foster, David Einstein, Herman Rubin, and Terry Moore
for pointing out that the above divergent series is given directly
by the Euler-Maclaurin summation formula. I hadn't immediately
recognized this as an EM summation because of the round-about path
by which I arrived at it. For any positive integers b and N you can
construct a Nth order linear recurring sequence approximating the
values of H(b^k), k=0,1,2,... For example, setting h(k) = H(2^k)
we have the (approximate) 3rd order recurrence
2h(k) - 5h(k-1) + 4h(k-2) - h(k-3) = 0
Taking the exact initial values h(0)=1, h(1)=3/2, and h(2)=25/12, the
closed-form solution of this recurrence gives the (approximate) values
2 2 1
h(k) = --- k + --- + -----
3 3 3*2^k
Of course, this is exact only for the three initial values (k=0,1,2),
so let's call this function h3(k). Similarly the fourth order
"harmonic recurrence" gives the function
218 110 17 64
h4(k) = --- k + --- + ------ - ----------
315 189 35*2^k 945*2^(2k)
and so on. Obviously the first j values of hj(k) are the exact values
h(k) of the harmonic series. Going to higher and higher order
formulas it quickly becomes apparent that the coefficient of k is
approaching ln(2), and the constant coefficient is approaching gamma.
For example, the constant coefficient of h7(k) is
330229884852951965355239887485346
--------------------------------- = 0.5772156...
572108310987828197293898319807675
Also, the coefficients of the subsequent terms approach B_n/n, so it's
natural to jump to the conclusion of equation (1)...until realizing
that the summation diverges.