The Museum of HP Calculators

HP Forum Archive 15

[ Return to Index | Top of Index ]

tricky definite integration
Message #1 Posted by hugh steers on 1 Feb 2006, 10:16 p.m.

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?
Re: tricky definite integration
Message #2 Posted by Marcus von Cube, Germany on 2 Feb 2006, 3:17 a.m.,
in response to message #1 by hugh steers


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.


Re: tricky definite integration
Message #3 Posted by GE on 2 Feb 2006, 3:45 a.m.,
in response to message #1 by hugh steers

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.

Re: tricky definite integration
Message #4 Posted by Valentin Albillo on 2 Feb 2006, 7:05 a.m.,
in response to message #1 by hugh steers

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))


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.

Re: tricky definite integration - help wanted to relearn
Message #5 Posted by Arnaud Amiel on 3 Feb 2006, 3:37 a.m.,
in response to message #4 by Valentin Albillo

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.



Re: tricky definite integration
Message #6 Posted by hugh steers on 3 Feb 2006, 12:13 p.m.,
in response to message #4 by Valentin Albillo

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


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:


which finally gives the very good (correct) answer:


Re: tricky definite integration
Message #7 Posted by Tizedes Csaba [Hungary] on 6 Feb 2006, 3:08 a.m.,
in response to message #1 by hugh steers


I played with the Taylor-approximation of

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...


[ Return to Index | Top of Index ]

Go back to the main exhibit hall