Inverse Of Harmonic Numbers
Leroy Quet
H(n) = H_n = 1 + 1/2 + 1/3 +...+ 1/n.
So, for our purposes, H(x) = sum{k=1 to oo} [1/k - 1/(k+x)] =
c + d ln(Gamma(x+1)) /dx, where c is Euler's constant, and Gamma(x+1)
is the Gamma function.
How could we calculate the inverse of H(x) (the J(x) such that J(H(x))
= x)?
The inverse of H(x) is multivalued.
But for each real x, there is only one value > -1. This is the branch
I'm interested in.
Thanks,
Leroy Quet
David W. Cantrell
As noted in recent articles of mine about the gamma function and its
inverse, a nice, simple approximation for the gamma function is
Sqrt(2*pi)*((x-1/2)/e)^(x-1/2). From this, it is easy to get a
corresponding approximation for the digamma function: psi(x) is
approximately Ln(x-1/2). We can then get a nice approximation of the
desired branch of the inverse of H(x) = g + psi(x+1). [This is equivalent
to Leroy's definition of H. I prefer to use g to denote the Euler gamma
constant in ASCII.]
InvH(x) is approximately e^(x-g) - 1/2.
How good is this approximation? Well, of course, it is not good for
small x. But as x increases without bound, the error (not just the
relative error) approaches 0. As examples,
e^(H(25)-g) - 1/2 = 25.0016 and
e^(H(1000)-g) - 1/2 = 1000.000042
I suppose, Leroy, that you wanted a precise inverse. Offhand, I can't
help with that; but better approximations could be devised if you're
interested.
Regards,
David W. Cantrell
--
-------------------- http://NewsReader.Com/ --------------------
Usenet Newsgroup Service
David W. Cantrell
> qqq...@mindspring.com (Leroy Quet) wrote:
> > Let H(x) be the continuous form of the harmonic numbers, where
> > H(n) = H_n = 1 + 1/2 + 1/3 +...+ 1/n.
> >
> > So, for our purposes, H(x) = sum{k=1 to oo} [1/k - 1/(k+x)] =
> > c + d ln(Gamma(x+1)) /dx, where c is Euler's constant, and Gamma(x+1)
> > is the Gamma function.
> >
> > How could we calculate the inverse of H(x) (the J(x) such that J(H(x))
> > = x)?
> > The inverse of H(x) is multivalued.
> > But for each real x, there is only one value > -1. This is the branch
> > I'm interested in.
>
> As noted in recent articles of mine about the gamma function and its
> inverse, a nice, simple approximation for the gamma function is
> Sqrt(2*pi)*((x-1/2)/e)^(x-1/2). From this, it is easy to get a
> corresponding approximation for the digamma function: psi(x) is
> approximately Ln(x-1/2). We can then get a nice approximation of the
> desired branch of the inverse of H(x) = g + psi(x+1). [This is equivalent
> to Leroy's definition of H. I prefer to use g to denote the Euler gamma
> constant in ASCII.]
>
> InvH(x) is approximately e^(x-g) - 1/2.
>
> How good is this approximation? Well, of course, it is not good for
> small x.
In particular, e^(1-g) - 1/2 = 1.026 .
Perhaps it's worth noting that, if one is interested in inverting just
the _discrete_ H(n), then we can get a precise result by merely rounding
to the nearest integer. In other words, for all positive integers n,
Round( e^(H(n)-g) ) = n
> But as x increases without bound, the error (not just the
> relative error) approaches 0. As examples,
>
> e^(H(25)-g) - 1/2 = 25.0016 and
>
> e^(H(1000)-g) - 1/2 = 1000.000042
>
> I suppose, Leroy, that you wanted a precise inverse. Offhand, I can't
> help with that; but better approximations could be devised if you're
> interested.
In fact, there's a better inverse which can be expressed using the
principal branch of the Lambert W function: InvH(x) is approximately
1/2*( 1/Sqrt(-3*W(-1/(12*e^(2*(x-g))))) - 1 ).
Evaluating this improved approximate inverse for x = H(1), H(25), and
H(1000), we obtain 0.9981, 24.99999956, and 999.9999999999927, resp.
Robert Israel
David W. Cantrell <DWCan...@sigmaxi.org> wrote:
>I suppose, Leroy, that you wanted a precise inverse. Offhand, I can't
>help with that; but better approximations could be devised if you're
>interested.
A couple of iterations of Newton's Method with David's approximation as
starting point ought to work very well. In Maple notation, the
Newton iteration for solving H(x)=y would be x -> x - (H(x)-y)/Psi(1,x+1).
e.g. for y = 1, H(1) = 1, David's approximation is 1.026205112,
one Newton iteration gives 0.9997859102, a second gives 0.9999999860.
Robert Israel isr...@math.ubc.ca
Department of Mathematics http://www.math.ubc.ca/~israel
University of British Columbia
Vancouver, BC, Canada V6T 1Z2
David W. Cantrell
> In article <20011111161607.086$K...@newsreader.com>,
> David W. Cantrell <DWCan...@sigmaxi.org> wrote:
>
> >I suppose, Leroy, that you wanted a precise inverse. Offhand, I can't
> >help with that;
Perhaps I can help with that now.
Two useful series are given toward the end of this article.
> >but better approximations could be devised if you're interested.
>
> A couple of iterations of Newton's Method with David's approximation as
> starting point ought to work very well. In Maple notation, the
> Newton iteration for solving H(x)=y would be x -> x - (H(x)-y)/Psi(1,
> x+1). e.g. for y = 1, H(1) = 1, David's approximation is 1.026205112,
> one Newton iteration gives 0.9997859102, a second gives 0.9999999860.
This will work well only if y is large enough. Remember that I said that
e^(x - g) - 1/2 does not approximate InvH(x) well if x is small. But I
suppose that you, Robert, were considering just the case when x is
sufficiently large. Nonetheless, it seems that Leroy wanted to invert
H(x) for x > -1, and so for completeness, I should be more thorough.
Also, in any event, I should have made clear (but did not do so) in my
first two articles here that:
(1) The approximation above is very poor for negative x, obviously
approaching -1/2 as x decreases without bound, rather than -1 as does
InvH(x).
(2) My other approximate inverse, in terms of the Lambert W function,
which is a much better approximation of InvH(x) for large x, is not even
defined for x < g + (1 - Ln(12))/2 or about -0.163 .
When x << 0, it is a bad idea to use my simple approximate inverse
above to obtain a starting value for Newton's method. But in this case,
a far simpler approximation, just -1 - 1/x, provides a reasonable
starting value. (Better starting values can be obtained using more terms
of a series, given below.)
As a modest example, suppose that x = -pi/2 - Ln(8), and note then that
InvH(x) = -3/4 precisely. Using Newton's method: With x_0 = -1 - 1/x, we
obtain successively -0.726, -0.7521, -0.750017, -0.7500000011, ...
whereas with x_0 = e^(x - g) - 1/2, we obtain successively -0.485, -0.982,
-0.965, -0.935, ..., x_8 = -0.750010, ... Had a substantially smaller x
value been chosen for the example, the contrast between use of the
different starting values would have more extreme.
Two useful series for InvH(x):
(1) Expanding at -oo gives InvH(x) =
-1 - 1/x + pi^2/(6*x^3) + Zeta(3)/x^4 - (2*pi^4)/(45*x^5)
+ (6*Zeta(5) - 5*pi^2*Zeta(3))/(6*x^6) + ...
(2) Expanding at 0 gives InvH(x) =
6/pi^2*x + 216*Zeta(3)/Pi^6*x^2
+ 72*(-pi^6 + 1080*Zeta(3)^2)/(5*pi^10)*x^3
- 2592*(pi^6*Zeta(3) - 540*Zeta(3)^3 - 3*pi^4*Zeta(5))/pi^14*x^4 + ...
In summary, we can approximate InvH(x) well (and simply) using
(a) e^(x - g) - 1/2 when x >> 0,
(b) several terms of the above Maclaurin series when x is close to 0, and
(c) several terms of the other series above, in powers of 1/x, when x << 0.