tricky definite integration

[hugh steers]

Suppose you wanted to evaluate the following definite integral:

      /1

     |     1

I=   | cos(-) dx

     |     x

    /0

the idea is to try to evaluate it on a calculator. so far, no calculator i've tried will

handle the input as stated, so i'm looking for a manipulation or substitution

that will break down the problem or transform it into a do-able calculation

that will yield over 5 significant figures. a bit like some of the hp

manuals treated some difficault cases.

But first: the answer!

We actually know the symbolic answer:

                                          /x

          1                              | Sin(t)

I=  x cos(-) + Si(x) + c     where Si(x)=| ------ dt

          x                              |   t

                                        /0

therefore I = cos(1) + Si(1) - Si(x)         and Si(x)     = pi/2

                               lim x->inf        lim ->inf

       so I = -0.084410950559573886889

Now suppose we didn't know the above and wanted to evaluate the answer on a

calculator and i don't mean using the breakdown of `I' in terms of Si(x) and

calculating Si(1) on a calculator because this is cheating.

let's try some substitions:



plan A:                /inf

                      |    cos(u)

 u = 1/x,   gives I = |   ------ du

                      |      2

                     /1     u

this isn't any good either because of the `inf'. you can try numbers like 1000

for inf and get some of the answer. for example, on a 71b i (eventually) get

integral(1,1000,1e-5,cos(ivar)/ivar^2) = -0.0844101253838 which is correct to

the 5 figures i asked of it (ie 1e-5), but i didn't prove the remainder of the

integral was small enough to ignore and getting more than 5 this way would

be very tedious.



plan B:                           /pi/2

                                 |     cos(tan(y))

 u = tan(y), gives         I =   |    ----------- dy

                                 |          2

                                /pi/4    sin(y)



now we get finite limits, but still the function oscillates so wildly,

there's insufficient calculator convergence. for instance,

integral(pi/4,pi/2,1e-3,cos(tan(ivar))/sin(ivar)^2) = -0.084561887...

after some time and this is still quite vague.

so does anyone have a transform that gets fast enough convergence?

[Marcus von Cube, Germany]

Hugh,

you'll find this thread interesting:
Torture Integral. Look for message #5.

This beast is even worse:

Integral between 0 and 1 of f(x) = Cos(Ln(x)/x)/x

The trick is to integrate along a complex path where the trigonometric functions stop oscillating.

Marcus

[GE]

Plan A is not so bad as you think, since we can find an upper bound for the remainer R=integral(cos(u)/u^2,1000,inf) : abs(R)<integral(1/u^2,1000,inf)=1/1000.

So you get your 4 exact digits.

I wonder if you can take u = cos(1/x), since arccos is relatively well behaved. I have not tried.

Thanks for this mind-streching post.

[Valentin Albillo]

Hi, Hugh:

The HP-71B can compute your integral as is, for instance let's get and display four decimal places:

     FIX 4 @ INTEGRAL(0,1,1E-4,COS(1/IVAR))

        -0.0844

but of course it's rather slow. The optimum way to compute it, exact analytic methods apart, is to do the integration in the complex plane, where your trigonometric function becomes just the real part of some complex exponential, and using a path of integration that mostly avoids the oscillations. A suitable change of variable before performing the complex integration might help as well, if the upper or lower limits are infinite. This is the standard "trick" and it does work.

The integral then becomes much easier and can be computed to 12-digit accuracy in a matter of minutes (or seconds under Emu71).

Best regards from V.

[Arnaud Amiel]

I am loking to refresh my memory on complex plane integration. Quite a few years ago, it was covered in my university maths course. Unfortunately it was on Friday afternoons so only remember it is very useful for tricky integral and that there are a few pitfalls in the selection of integration path.

Any pointers to go and relearn this would be helpful.

Thanks

Arnaud

[hugh steers]

thanks for your complex contour suggestion,

i decided to try out the contour approach. i was hoping to find a really

cunning path that might make the problem really easy. but i can't see one

(anyone?). so here's my attempt:

working now in the complex plane, the problem i have is to start from (0,0)

which is a pole. so what i'm going to do is not start on the pole and instead

start on (eps,0) [eps = epsilon].

my path from (eps,0) to (1,0) has three parts:

part A: quarter circle anticlockwise from (eps, 0) to (0, eps).

part B: line along y, from (0,eps) to (0,1).

part C: clockwise quarter cicle from (0,1) to (1,0).

so my answer, I, is Re(A + B + C), providing i can make eps->0 without

problems (Re = real part).

now, i can say right away that Re(B) = 0, since it is entirely on the

imaginary line. which leaves me with two problems lim eps->0 Re(A) and Re(C).

this is where the calculator can only partly help me. i have to do part A

and maybe it (the 71b in this case) can do part C.

                1

let  f(z) = cos(-)

                z

PART A:

                /pi/2

               |                                                            it

       A =     | f(z(t)) z'(t) dt   and z(t) = eps(cos(t)+ i sin(t)) = eps e

               |

              /0


                   /pi/2

                  |        1       i t

         = i eps  | cos(--------) e    dt

                  |          i t

                 /0     eps e

now, the bit i want is Re(A) whilst eps->0. however this is where i cheat

because i can't let eps->0 before integrating and i don't know how to

integrate it.

now it turns out that Re(A) = -eps cos(1/eps) - Si(1/eps), where Si(z) is the

sine integral as before. knowing this, it's clear than Re(A) lim eps->0

 = 0 - si(inf)

 = -pi/2

but that was because i knew the original indefinite integral all along.

maybe i can put a small eps into the 71b and integrate numerically and guess

the answer is pi/2. is there any way to deal with the limit before

integrating??  so, by guessing or otherwise i have Re(A) = -pi/2

PART C:

this part i have to manipulate the problem into a form i can feed to the 71b.

as in part A, but with no eps, i have:


              /pi/2

             |      -(i t)   i t

    C =  -i  | cos(e      ) e    dt   [note: minus comes from other direction]

             |

            /0



the first idea was to feed this directly to the 71b as a complex numerical

integration, like this:

integral(0,pi/2,1e-5,cos(exp((0,-1)*ivar))*exp((0,1)*ivar))

but it doesn't like it. its a pity, i was hoping it could perform its

summation with complex numbers. is there a way??

so, i have to do more more work. i can expand C into real and complex

parts. if i expand the integrand and take the imaginary part, it will form the

real part of the answer when put with the -i outside. so after simplification,

i get:

       /pi/2

      |

Re(C)=| (sin(t) cos(cos(t)) cosh(sin(t)) + cos(t) sin(cos(t)) sinh(sin(t))) dt

      |

     /0

and into the 71b (with the additional -pi/2 from part A) as this:

10 DEF FNA

20 S=SIN(IVAR) @ C=COS(IVAR)

30 FNA=S*COS(C)COSH(S)+C*SIN(C)*SINH(S)

40 END DEF

50 DISP INTEGRAL(0,PI/2,1E-8,FNA)-PI/2

which finally gives the very good (correct) answer:

-0.08441095058.

[Tizedes Csaba [Hungary]]

Hello,

I played with the Taylor-approximation of

1/X

function on my 48SX. Some results:

1/X = 1+SUM(I=1,INFINITY,(1-X)^I) = (1-(1-X)^(N+1))/X (where N is big enough, but not too ;) )

I know, this is not an elegant way, but I spend some time to play with it. I want only a fastest result. I wrote few versions of this method, finally the "PRG1" was the integrator with 'COS(INV(X))' intengrand, and "PRG2" was the integrator with 'COS((1-(1-X)^(N+1))/X)'. The big fault is the numerical error when the N is too big!

--------------------------------------------------------

N=24  RESULTS                RUNNING TIMES

      PRG1       PRG2        PRG1         PRG2

FIX1  -0.07725   -0.19917       1.2sec     1.9sec

FIX2  -0.083396  -0.084487     37.6sec    14.7sec

FIX3  -0.084375  -0.084487   1193  sec(!) 14.1sec(!)

--------------------------------------------------------

N=49  RESULTS                RUNNING TIMES

      PRG1       PRG2        PRG1         PRG2

FIX1 (AS ABOVE)  -0.06555    (AS ABOVE)    1.9sec

FIX2 (AS ABOVE)  -0.085261   (AS ABOVE)   28.7sec

FIX3 (AS ABOVE)  -0.085261   (AS ABOVE)   28.7sec

--------------------------------------------------------

I tried some other way too: changing the boundaries of the interval from [0...1] to [1...H], where H->0, tried to integrate 1+COS(INV(X)), then substract 1 from the result (with this method no change sing on the interval), and so on...

Csaba