HP48 integration error?



#10

When I evaluate the integral of exp(-x^2) from 0 to 1000 [by using the equation writer, then pressing "Evaluate"], the HP48 gives an answer of zero. Anyone know why this happens? When I evaluate the integral from 0 to 100, the HP48 gives an answer of 0.8862, which is the right answer.


#11

Quote:
When I evaluate the integral of exp(-x^2) from 0 to 1000 [by using the equation writer, then pressing "Evaluate"], the HP48 gives an answer of zero. Anyone know why this happens? When I evaluate the integral from 0 to 100, the HP48 gives an answer of 0.8862, which is the right answer.

The exact answer is 0.5 * sqrt(pi). The reason for the erroneous answer of zero for a large upper limit of integration is sparse sampling. Too many values of the integrand are extremely small when the samples are stretched out to x = 1000, such that the calculated integral doesn't change very much when additional widely-spaced samples are taken. There's absolutely no point in setting the upper limit greater than 34, because the value of the integrand is less than 1E-499 from x = 33.89675 onward. In fact x = 6 would suffice as an upper limit.

This issue might not be covered in the HP-48 manuals, but a very similar problem is covered in the HP-15C Advanced Functions Handbook.

-- KS

Edited: 16 Mar 2007, 12:21 a.m.

#12

I echo Karl's answer.

There is a lot of talk in the 15C manuals about how the algorithm samples the "interesting" part of the curve, and if there is lots of uninteresting tail stuff, these samples predominate. exp(-x^2) goes to zero for big x, so lots of zeros get the weight.

Improper integrals are often handled best by a transformation of variables. In the case of the bell curve, integrating 0.5/Sqrt[-Log[t]] from 0 to 1 is equivalent to integrating your problem from zero to infinity--still an improper integral with poles at each end, but a finite interval, which is nicer. I leave it to you to determine the variable transform that leads from one to the other, because I am a lazy devil and it is late.

HTH,

Les

Edited: 16 Mar 2007, 4:03 a.m.


#13

Actually, I think in this case cutting off the tail as Karl suggested is faster and easier.

Where the tranform could be useful is if you want to find the area of the tail itself, from x to infinity, as in the computation of the complementary error function. But if you want erfc(x), SQ 2 * 1 SWAP UTPC after putting your x on the stack is a lot faster. As a matter of fact, the following little code snippet of mine handles all real arguments:

<< DUP SQ 2 * 1 SWAP UTPC IF SWAP 0 < THEN 2 - NEG END >>

HTH,

Les

Edited: 16 Mar 2007, 7:11 p.m.

#14

Les, Karl,

Thank you for your response, and your suggestions on evaluating the integral appropriately.

Still, what seems strange is that when evaluating the integral from 0 to 35, it takes about 5 seconds for the HP48G to yield the (correct, to specified decimal places) answer.

However, when evaluating the integral from 0 to 1000, the HP48 in a split second returns a zero answer. I would think it would take longer to evaluate the integral from 0 to 1000, compared to the time is took to evaluate the integral from 0 to 35, which is not the case.

Any thoughts on this?

In any case, when using the HP integrate function, it is clear the user must make sure the answers make sense.

This issue on evaluating integrals is not discussed in the HP48G Users Manual. Is there anyway you can copy the relevant excerpts from the HP15C manual into a response to this message?

Any help you could provide would be appreciated.

Thanks - kelly


#15

Yup, I have an idea.

The integration algorithm makes an estimation at a few points, then doubles the number of points and makes another estimation, refines this by something called Richardson extrapolation, etc., and continues until two or three subsequent estimations are equal to each according to the tolerance implied by the setting for FIX, SCI, or ENG.

When estimating this integral, the value of the integrand is < 1e-499 for the huge majority of the range 0 to 1000. Basically, this is zero as far as the calculator is concerned. The value of integrand only becomes nonzero, as far as the calculator is concerned, somewhere between x = 33 and x = 34. That part of the range of integration is stuck way on the left end, and it is unlikely that in a few sequential estimates the algorithm even samples much from that end. This means that you are going to get two or three consecutive estimates that are effectively zero, since the integral estimates are basically weighted averages of the integrand at the sample points. The weighted average of a bunch of zeros is always zero. This happens a a few times in a row (I think at least three) and tells the calculator convergence has been reached. That happens really quickly, right at the beginning, so false convergence is reached quickly.

The stuff in the 15C manuals is not conveniently copied in this forum, as there is pages of stuff. But, in a general sense, if you want to basic idea of what the 48G does in integration, google Romberg Method or Romberg algorithm. The math is not very hard, even though the notation can be confusing--subscripts to keep track of and that sort of thing. If you can appreciate graphically this sampling issue, you may grasp why the HP48G goes to 0 when you set up the integral with the 0-1000 limits.

Les

Edited: 17 Mar 2007, 8:32 a.m.


#16

Actually, trying plotting exp(-x^2) directly on the calculator. Set the vertical range to -.1..1 and the horizontal to 0..1000.

When you execute DRAW you won't see much--mostly blank plot with maybe some extra pixels up against the y-axis. Even with the range set to 0..35 the visible part of the curve is scrunched up to the left, but the integrator will still sample some points over there. If you play around with the upper limit and replot, you will see why one can chop off the tail and still get a good estimate of the integral--there is just not much area under that tail once you get past x = 4, even. Choosing x = 35 as an upper limit is sort of overkill.

Plotting a function before trying to integrate it always makes sense. And do try to learn about the Romberg method. I am of the school that no one who uses numerical integration should treat them just as black box routines.

Les

Edited: 17 Mar 2007, 8:44 a.m.


#17

Quote:
And do try to learn about the Romberg method

Just a warning--most theoretical discussions of the Romberg method seem to talk about it as extrapolating the trapezoid rule. In HP calculators since the 34C, the Romberg integrator actually use the midpoint rule (after certain transformations of variable), so that the endpoints of the interval of integration are never directly sampled.

For example, if integrating on the interval -1..1, the first iteration will sample at 0, the next at +1/2 and -1/2, the next at +1/4, +3/4 and the negative of these, the next at +1/8, +3/8, +5/8, +7/8 and the negative of these, etc.

On the other hand, the trapezoid rule would sample the endpoints 1 and -1 first, then the 0 between them, then the series of points given above, BUT the estimate integral would be trapezoid estimates based on the integrand values at the end of each subinterval, not midpoint estimates based on the middle point of each subinterval.

It really is a lot easier to look at pictures or draw this stuff out.

Les

#18

Les,

Your explanation below (shown in quotes) makes sense. Thank you. Guess I had to know why the HP was yielding the zero answer, even though the Calculus was telling me it was not zero.

Cheers - Kelly


"The integration algorithm makes an estimation at a few points, then doubles the number of points and makes another estimation, refines this by something called Richardson extrapolation, etc., and continues until two or three subsequent estimations are equal to each according to the tolerance implied by the setting for FIX, SCI, or ENG.

When estimating this integral, the value of the integrand is < 1e-499 for the huge majority of the range 0 to 1000. Basically, this is zero as far as the calculator is concerned. The value of integrand only becomes nonzero, as far as the calculator is concerned, somewhere between x = 33 and x = 34. That part of the range of integration is stuck way on the left end, and it is unlikely that in a few sequential estimates the algorithm even samples much from that end. This means that you are going to get two or three consecutive estimates that are effectively zero, since the integral estimates are basically weighted averages of the integrand at the sample points. The weighted average of a bunch of zeros is always zero. This happens a a few times in a row (I think at least three) and tells the calculator convergence has been reached. That happens really quickly, right at the beginning, so false convergence is reached quickly."


Possibly Related Threads...
Thread Author Replies Views Last Post
  Integration question and "RPN" mode comment Craig Thomas 16 1,296 12-05-2013, 02:32 AM
Last Post: Nick_S
  MS advert shows spreadsheet with obvious error BruceH 3 440 11-14-2013, 09:50 AM
Last Post: Bill (Smithville, NJ)
  HP Prime: Rounding error in determinant Stephan Matthys 3 347 10-25-2013, 09:29 PM
Last Post: Walter B
  Prime Error or Mine? toml_12953 12 849 10-22-2013, 10:35 AM
Last Post: toml_12953
  WP34s integration trapped in infinite loop Bernd Grubert 25 1,413 10-17-2013, 08:50 AM
Last Post: Dieter
  Explaination on How to Reset Caculator in Users guie: error Harold A Climer 5 475 10-15-2013, 02:11 AM
Last Post: cyrille de Brébisson
  Repair of HP-34C - Error 0 and Error 9 Jeff Kearns 3 340 10-11-2013, 12:29 PM
Last Post: Randy
  HP Prime integration Richard Berler 1 261 10-06-2013, 10:52 PM
Last Post: Helge Gabert
  Do You Think a Listing Error Was Made? Jim Johnson 13 819 09-04-2013, 09:23 AM
Last Post: Eddie W. Shore
  HP Prime: Error reports : Here Patrice 111 4,056 07-24-2013, 05:52 PM
Last Post: Thomas Klemm

Forum Jump: