Handheld Calculator Helps Solving Kahan's Integral Symbolically



#5

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_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/(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.


#6

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

#7

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.


#8

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 605 11-21-2013, 08:25 PM
Last Post: Chris Pem10
  HP-PRIME CAS SOLVING fabrice48 8 618 10-19-2013, 01:21 PM
Last Post: Han
  HP Prime Solving Nonlinear System of Equations for Complex Results Helge Gabert 11 906 09-30-2013, 03:44 AM
Last Post: From Hong Kong
  Integral Richard Berler 32 1,693 09-22-2013, 08:59 PM
Last Post: Les Koller
  HP charger 82002c for handheld calculators Dave Burson 2 275 05-19-2013, 12:51 PM
Last Post: Dave Burson
  [WP34s et al.] Solving the TVM equation for the interest rate Dieter 24 1,386 12-01-2012, 05:53 AM
Last Post: Paul Dale
  Solving with the built-in equations of the 35s Palmer O. Hanson, Jr. 4 341 10-17-2012, 03:17 PM
Last Post: Dieter
  Another tough integral Gerson W. Barbosa 59 1,984 05-20-2012, 05:59 PM
Last Post: Gerson W. Barbosa
  can't calculate integral on hp 41 Hans Holzach 12 701 05-16-2012, 02:59 PM
Last Post: Hans Holzach
  Let's see if HP helps out Matt Agajanian 13 590 04-27-2012, 12:54 PM
Last Post: bill platt

Forum Jump: