 Handheld Calculator Helps Solving Kahan's Integral Symbolically - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: Handheld Calculator Helps Solving Kahan's Integral Symbolically (/thread-109056.html) Handheld Calculator Helps Solving Kahan's Integral Symbolically - Gerson W. Barbosa - 02-25-2007 The Kahan's integral has appeared in some recent threads, like this one, started by Les Wright: 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_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): /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. Re: Handheld Calculator Helps Solving Kahan's Integral Symbolically - Karl Schneider - 02-25-2007 Hi, Gerson -- That's very good work. I, too, had attempted to derive the exact answer, but didn't know of an "insightful" way to integrate 1/ln(x). I'd certainly seen the Euler-Mascheroni Constant before, but didn't know its applications. It should be an MoHPC article, or part thereof. -- KS Re: Handheld Calculator Helps Solving Kahan's Integral Symbolically - Valentin Albillo - 02-26-2007 Hi, Gerson: First of all, excellent work ! And a truly didactic presentation, I'd say. Just a comment: Gerson posted: "E appears in various integrals. One of them is the classical form for Euler's constant [...] /1 | E = | (1/ln(x) + 1/(1-x))dx /0 " Personally, I find this other integral form simpler and more aesthetically pleasing: /1 | E = | -ln(-ln(x)) .dx /0 where you have the curiosity of an iterated logarithm (the log of a log) and further, not only two logs but also two minus signs (and nothing else!), one of them applied to the argument of the outer log, which at first glance can trick you into thinking about logs of negative arguments, thus complex numbers, which isn't the case of course. Best regards from V. Re: Handheld Calculator Helps Solving Kahan's Integral Symbolically - Gerson W. Barbosa - 02-26-2007 Hello Valentin, Quote: Personally, I find this other integral form simpler and more aesthetically pleasing: /1 | E = | -ln(-ln(x)) .dx /0 This is really an elegant form, and easier to remember. On the 48G/GX and the 71B, I think, it can be evaluated numerically to give the first nine digits of E (For plain ASCII text y might be a better notation). Special functions and constants are interesting as they can provide closed form solutions to problems previously only being able to be solved numerically, as the HP-67 Math Pac I example in your excellent HP-71B Math ROM Baker's Dozen (Vols. 1 & 2) articles. Perhaps the next version of the HP-50G CAS should include Lambert's W-function. Best regards, Gerson.