Handheld Calculator Helps Solving Kahan's Integral Symbolically Message #1 Posted by Gerson W. Barbosa on 25 Feb 2007, 1:12 p.m.
The Kahan's integral has appeared in some recent threads, like this one, started by Les Wright:
http://www.hpmuseum.org/cgisys/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, HewlettPackard Jounal, Aug/1980, pages 2332.
This is the famous integral:
/1

I =  (sqrt(x)/(x1)  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.738E3 then decreases more smoothly to 0 as x approaches 1. This was troublesome for the HP Solver on the HP34C, then on the HP15C and even on the HP32SII, 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 HP34 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 HP50G 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 EulerMascheroni 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/EulerMascheroni_constant
That 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/(1x))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/(1x))dx
/0 /0
Replacing this in Kahan's integral:
/1

I =  (sqrt(x)/(x1) + 1/(1x))dx  E
/0
This should be easy to integrate, but let's the HP49G/G+/50G do it:
'\v/X/(X1)+1/(1X)' INTVX > '2*\v/X+(LN(\v/X+1)+LN(ABS(\v/X1)))LN(ABS(X1))'
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 HP49G/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 HP50G, we get
0.03648997398
On the HP200LX 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 HP33S:
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.
