Handheld Calculator Helps Solving Kahan's Integral Symbolically « Next Oldest | Next Newest »

 ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 02-25-2007, 01:12 PM 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. ▼ Karl Schneider Posting Freak Posts: 1,792 Threads: 62 Joined: Jan 2005 02-25-2007, 03:48 PM 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 Valentin Albillo Posting Freak Posts: 1,755 Threads: 112 Joined: Jan 2005 02-26-2007, 05:24 AM 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. ▼ Gerson W. Barbosa Posting Freak Posts: 2,761 Threads: 100 Joined: Jul 2005 02-26-2007, 11:28 AM 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.

 Possibly Related Threads... Thread Author Replies Views Last Post [HP Prime] Tips for Solving Differential Equations More Efficiently Chris Pem10 8 2,295 11-21-2013, 08:25 PM Last Post: Chris Pem10 HP-PRIME CAS SOLVING fabrice48 8 2,279 10-19-2013, 01:21 PM Last Post: Han HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 3,352 09-30-2013, 03:44 AM Last Post: From Hong Kong Integral Richard Berler 32 5,911 09-22-2013, 08:59 PM Last Post: Les Koller HP charger 82002c for handheld calculators Dave Burson 2 1,227 05-19-2013, 12:51 PM Last Post: Dave Burson [WP34s et al.] Solving the TVM equation for the interest rate Dieter 24 4,725 12-01-2012, 05:53 AM Last Post: Paul Dale Solving with the built-in equations of the 35s Palmer O. Hanson, Jr. 4 1,152 10-17-2012, 03:17 PM Last Post: Dieter Another tough integral Gerson W. Barbosa 59 9,336 05-20-2012, 05:59 PM Last Post: Gerson W. Barbosa can't calculate integral on hp 41 Hans Holzach 12 2,825 05-16-2012, 02:59 PM Last Post: Hans Holzach Let's see if HP helps out Matt Agajanian 13 2,957 04-27-2012, 12:54 PM Last Post: bill platt

Forum Jump: 