The Kahan's integral has appeared in some recent threads, like this one, started by Les Wright:
http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv016.cgi?read=104550
Perhaps it should be more appropriately named Kahan/Woodyard Integral, as did Warren Anderson in a recent thread, since it was presented to Prof. Kahan for analysis by his electrical engineering colleague, Prof. J. R. Woodyard, as reported in Kahan's article
Handheld Calculator Evaluates Integrals, Hewlett-Packard Jounal, Aug/1980, pages 23-32.
This is the famous integral:
/1
|
I = | (sqrt(x)/(x-1) - 1/ln(x))dx
/0
The function is steep near 0: it is 0 when x=0, rises suddenly to 0.1177 when x=4.738E-3 then decreases more smoothly to 0 as x approaches 1. This was troublesome for the HP Solver on the HP-34C, then on the HP-15C and even on the HP-32SII, at FIX 5 display mode. Those calculators give I1 = 0.03662 +/- 0.00001, instead of the correct result at this accuracy, 0.03649, an "error too small to be obvious and too large to ignore", in Kahan's words. However, at FIX 7 the problem disappears, but the HP-34 takes 23 minutes to display the right answer. Prof. Kahan then presents a technique to tame what he calls "wild integrals", which is detailed in his article. Even the newest HP-50G takes 37 minutes and 30 seconds to compute the integral in its original form to maximum accuracy. It places the result between 0.0364899739827 and 0.0364899739835, the lower limit being slightly above the exact answer as we will verify later.
Kahan's integran can be solved exactly with help of the EulerGamma or Euler-Mascheroni constant. According to some references, its the third most important constant in Mathematics, after pi and e. It's represented by the greek letter gamma. Because the lack of this letter here, we'll represent it as simply E, for short.
It is defined as:
E = lim (1 + 1/2 + ... + 1/n - ln n)
n->inf
E = 0.577215664901533, to fifteen places. More about E at
http://en.wikipedia.org/wiki/Euler-Mascheroni_constantThat said, let's try to find an exact result for Kahan's integral. As mentioned in the references, E appears in various integrals. One of them is the classical form for Euler's constant, as shown in this paper (Formula 9, page 4):
http://arxiv.org/PS_cache/math/pdf/0306/0306008.pdf
/1
|
E = | (1/ln(x) + 1/(1-x))dx
/0
Hence we can compute the second term in Kahan's integral, the one we cannot normally compute symbolically:
/1 /1
| |
| (1/ln(x))dx = E - | (1/(1-x))dx
/0 /0
Replacing this in Kahan's integral:
/1
|
I = | (sqrt(x)/(x-1) + 1/(1-x))dx - E
/0
This should be easy to integrate, but let's the HP-49G/G+/50G do it:
'\v/X/(X-1)+1/(1-X)' INTVX -> '2*\v/X+(-LN(\v/X+1)+LN(ABS(\v/X-1)))-LN(ABS(X-1))'
where \v/ is the square root symbol.
| |x=1
| |
I = | 2 sqrt(x) - ln(|sqrt(x) + 1|) + ln(|sqrt(x) - 1|) - ln(|x - 1|) | - E
| |x=0
We should take care not to cancel the last two terms when replacing x by the integration limit 1. Instead, we should replace them by ln[|(sqrt(x) - 1)/(x - 1)|) and find that the limit when x tends to 1 is 1/2, by applying l'Hopital's rule, for instance.Finally, we obtain:
I = 2 - ln(2) + ln(1/2) - E
or
I = 2 - ln(4) - E
We are now ready to compute the numeric result of Kahan's integral. There is no EulerGamma in the HP-49G/G+/50G keyboard. But there is a special function under MTH menu, digamma function, Psi. In this reference, we notice that
E = - Psi(1) (formula 16)
Thus,
I = Psi(1) + 2 - ln(4)
On the HP-50G, we get
0.03648997398
On the HP-200LX we can compute
I = 2 - ln(4) - 0.577215664901533 = 0.036489973978576
An approximation good to 12 digits for E is
E ~ 1/sqrt(3) - 1/sqrt(55192773)We can check it on the HP-33S:
2 - 1/sqrt(3) + 1/sqrt(55192773) - LN(4) = 0.03648997398
Gerson.
P.S.: I am no expert in these topics to fully understand the references I have mentioned. My goal here was just to provide an exact answer to Kahan's integral.
Edited: 25 Feb 2007, 2:12 p.m.