 The Museum of HP Calculators

HP Forum Archive 16

 So your HP can INTEGRATE ...Message #1 Posted by Valentin Albillo on 3 July 2006, 6:58 a.m. Hi, all: There's recently been some threads related to numerical integration and the respective effectiveness and robustness of such methods as Simpson's Rule, Trapezoidal Rule, Gaussian quadrature, and Romberg's method as implemented by HP in the HP-11C, HP42S, HP-41C's Advantage ROM and HP-71B's Math ROM, to name a few models and plug-ins. As a mention was made to a certain difficult-to-numerically-compute integral I posted a few months ago, and lest you think that HP's Romberg implementation is so robust that either just by itself or with a little bit of help on the user's part, it can cope with most anything thrown at it, go take your best HP hardware and software and boldly set out to try and numerically compute these six little integrals o'mine. By the way, no great precision is required to prove the point, just getting as few as 4 correct digits for each will do. ``` /1 I1 = | cos(x)*ln(x) .dx /0 /2*Pi I2 = | ln(1+x)*sin(10x) .dx /0 /2 I3 = | xPi/4*sin(Pi/(8-4x)) .dx /0 /Inf I4 = | sin(x)*sin(x2) .dx /0 /1 I5 = | sin(1/x)/x .dx /0 /1 I6 = | cos(ln(x)/x)/x .dx /0 ``` I'll give numerical answers to all of them in a few days, see if you can deliver before I do, and get some fun in the process. At the very least, they'll be useful to have at hand in order to test new integrators or integrating algorithms you might eventually get involved with. Best regards from V.

 Re: So your HP can INTEGRATE ...Message #2 Posted by John Smitherman on 3 July 2006, 8:39 p.m.,in response to message #1 by Valentin Albillo Valentin, the first three solved quickly on an HP33s, 49G+ and TI83+SE yielding the following: 1. -0.9461 2. -0.1976 3. 1.0108 The last three would not solve on the TI83+SE resulting in a tolerance exceeded message. I tried various large numbers for infinity for problem number four. I do not believe that the last three would solve on the 33s within my lifetime and so I interrupted the execution after many minutes. I tried various large numbers for infinity for problem number four. I was not sure how to handle infinity on the 49G+ for problem number four. I tried using the infinity symbol as well as several large numbers and got very different results. Any clues as I am relatively new to the 49G+? The 49G+ gave the following results for number five and six as follows: 5. -1.6117 6. 0.1691 I am hoping for partial credit teacher ;-) Regards, John

 Re: So your HP can INTEGRATE ...Message #3 Posted by Valentin Albillo on 4 July 2006, 4:15 a.m.,in response to message #2 by John Smitherman Hi, John: John posted: ```"1. -0.9461 2. -0.1976 3. 1.0108" ``` Right, right, and not-so-right. "Any clues as I am relatively new to the 49G+?" Never used one, can't say. "The 49G+ gave the following results for number five and six as follows: ``` 5. -1.6117 6. 0.1691" ``` Dead wrong, dead wrong. "I am hoping for partial credit teacher ;-)" You've got 4.75/10 so far. You can do better ... ! :-) Thanks for your interest and results and Best regards from V. Edited: 4 July 2006, 4:51 a.m.

 Re: So your HP can INTEGRATE ...Message #4 Posted by Karl Schneider on 4 July 2006, 3:47 a.m.,in response to message #1 by Valentin Albillo Quote: ...and lest you think that HP's Romberg implementation is so robust that either just by itself or with a little bit of help on the user's part, it can cope with most anything thrown at it... Thou dost smite me... :-) Well, just for the record, I meant the term "robust" to describe the practical soundness of the approach, not necessarily that the mathematical algorithms employed were the most capable, most sophisticated, or the quickest possible. Perhaps "rigorousness" would have been a better term than "robustness" I have a few answers, and hope to obtain the rest: I1: f(x) = cos(x)*ln(x)dx from 0 to 1 To enhance accuracy and speed, I substituted u = ln x, which yielded sin(x)*ln(x) - integ [sin(x)/x dx] from 0 to 1. Limit analysis on [sin(x)/x]*[x*ln(x)] using L'Hopital's Rule showed that sin(x)*ln(x) = 0 for both x = 0 and x = 1, and integrating the sinc function on an HP-48G with FIX 6 yielded I1 = -0.94608307068 Integrating the original function directly with FIX 6 yielded I1 = -0.94608289406 with more computational time. I2: f(x) = ln(1+x)*sin(10x) from 0 to 2*pi. This integrated directly and quickly on an HP-42S to I2 = -0.197627 Too easy. Did I miss something? I5: f(x) = sin(1/x)/x from 0 to 1 f(x) = sin(1/x)/x has high frequency and exploding magnitudes near x = 0. f(x) benefits from the change of variable u = 1/x Leading to integ [f(u) = sin(u)/u du] from 1 to infinity. A second change of variable w = u2 gives an integrand with lower frequency and more attenuation as w approaches infinity: f(w) = sin(sqrt(w))/(2w) from w = 1 to infinity Taking it in subdivisions up to w = 107 and FIX 6, I5 = 0.62479588 -- KS Edited: 4 July 2006, 3:56 a.m.

 Re: So your HP can INTEGRATE ...Message #5 Posted by Valentin Albillo on 4 July 2006, 5:45 a.m.,in response to message #4 by Karl Schneider Hi, Karl: Karl posted: "Thou dost smite me... :-)" Nay, but this man hath cause enow to smite ... :-) "Well, just for the record, I meant the term "robust" to describe [...] Perhaps "rigorousness" would have been a better term than "robustness"" Details, details ... "I1 = -0.94608307068" I only asked for 4 correct digits, so Ok. "I2 = -0.197627. Too easy. Did I miss something?"" No, it's just that I wanted to include an easy one so people would be able to get at least one correct result, even if doing it with Simpson's rule on an HP-25. Avoid frustrated readers. Good marketing technique. Ah, and yes, your result, rounded to 4 digits, is correct as well. "I5 = 0.62479588" This rounds to 0.6248, which is wrong by one ulp. You can do better. So far, a 4.9/10 for you. You *can* do better ! :-) Thanks a lot for your interest and detailed explanations, and Best regards from V.

 (Edited) So your HP can INTEGRATE ...Message #6 Posted by Karl Schneider on 6 July 2006, 2:41 a.m.,in response to message #5 by Valentin Albillo One more result, and some discussion: I3: ```f(x) = xpi/4* sin(pi/(8-4x)); x_lo = 0; x_hi = 2 ``` The integrand is positive from x = 0 to x = 7/4; the integral of that subinterval came out to 1.103119 at FIX 6. The remainder from x = 7/4 to x_hi produces high-frequency oscillations, and took quite a while to run on an HP-49G; the result at FIX 6 was -0.091867. (NOTE: This value, as well as the resulting sum, were corrected from the "-0.0918967" that appeared originally, due to clumsy fingers. The reader may refer to VA's response below.) The sum of these partial results was 1.101252 (which appears notto meet the "4 decimal digits" standard, if Rodger Rosenbaum's answer, obtained using more andvanced techniques, is correct). I4: ```f(x) = sin(x)*sin(x2); x_lo = 0; x_hi = infinity ``` This integrand is bounded between -1.00 and +1.00, but does not attenuate (converge) to zero as x approaches infinity. In fact, the integrand will probably take any and all values over that range. The upper limit of the variable of integration is unbounded. Thus, I am not yet convinced that the integral converges or can be proved to be bounded within a narrow specified range for all values of x beyond some value. Here's something interesting: Perform one integration by parts of integral (f(x).dx). The result is ```-cos(x)*sin(x2) + 2*integral (x*cos(x)*cos(x2).dx) ``` As x approaches infinity, the first term is bounded between -1.00 and +1.00, but does not approach a limit; the integrand is unbounded. I6: ```f(x) = cos(ln(x)/x)/x); x_lo = 0; x_hi = 1 ``` I remember this as one of the "torture" integrals discussed during December 2005, where you posted an approach based on complex numbers. Conventional approaches appear to be hopeless (and indeed they were, as described in the thread from that month.) I see a substitution that might make it converge a bit faster, but nothing to "tame" the integrand as x approaches zero. -- KS Edited: 7 July 2006, 1:11 a.m. after one or more responses were posted

 Re: (More) So your HP can INTEGRATE ...Message #7 Posted by Valentin Albillo on 6 July 2006, 4:13 a.m.,in response to message #6 by Karl Schneider Hi again, Karl: Karl posted: "I3: [...] the integral of that subinterval came out to 1.103119 at FIX 6. [...] -0.0918967. [...] The sum of these partial results was 1.101252" Strange, very strange ... On all my RPN calculators and my HP-71B (not to mention my SHARPs), the sum of your partial results comes out as: ` 1.103119 - 0.0918967 = 1.0112223` and not 1.101252, as your 49G gets ... Must be RPL at work again ... next time I suggest you try an RPN model, they might not include a built-in "Gas Compresibility Z Factor Function" (ZFACTOR), but sure they do add up easily and correctly ;-). Note to readers: if the above paragraph doesn't make sense, it's because Karl cheated and edited his post after I pointed out the blatant error :-) "I4: [...] I am not yet convinced that the integral converges [...]" Trust me: it does. Unlike some moron who posted an insolvable 'challenge' a while ago, mine are always solvable, save typographical errors of course. "I6: [...] I'm reticent to try it, even on an HP-48, HP-49, or HP-71B/Math ROM." Back in that thread I demonstrated how you could easily compute it to full accuracy and fairly quickly on an HP-71B (the real thing, not the 'light' version). If you run out of ideas, you can have a look at my solution back then and score an easy point :-) Thanks again for your unfailing interest in my 'challenges' and your always interesting (and didactic!) contributions. You can be 100% sure that your inputs are much appreciated by a substantial number of lurkers that never post anything but enjoy reading our ramblings. I'll give the solutions tomorrow, so you've got 24 additional hours to try and correctly solve all remaining integrals. Best regards from V.

 Re: (More) So your HP can INTEGRATE ...Message #9 Posted by Valentin Albillo on 7 July 2006, 4:49 a.m.,in response to message #8 by Karl Schneider Hi, Karl: Karl posted: "On I5: [...]You said, 'Trust me: it does' This brings to mind a statement [...]: 'Trust, but verify' [...] I did find the contradiction a bit amusing..." It is amusing to me as well. Seriously, Karl, I've never posted an insolvable challenge, not even on my April 1st S&SMC. And the integral is typographically simple enough that I would have detected any typo at once. So you can rest assured that it does indeed converge and it has a finite, small value. "It certainly seems that the integral would be finite and of small magnitude. However, I'd like to see the convergence or tight bounding proved. [...] However, there's apparently a closed-form solution, according to Rodger" Yes, there's a generalized closed-form solution to a family of integrals of which this is but one of the simplest particular cases. As for proving convergence, I'll give you a few hints for you to have a start at proving it: Using elementary trigonometric identities, express the product as a sum of two cosines. Using a trivial, implied variable substitution, see that both cosines are equivalent regarding convergence, so that convergence of one of them implies that of the other, thus it suffices to prove convergence for one of them. Using a very simple variable substitution, get rid of both the square and the linear term inside the cosine (very similar to that used to solve a polynomial quadratic equation), and you'll find yourself with an integral which is a cosine divided by the square root of a linear polynomial. Integrating by subintervals between the zeros of the cosine, you'll see that you get an alternating series where the general term tends to zero (like, say, the well-known logarithmic series: 1-1/2+1/3-1/4+...). This fact (i.e, alternating series and general term tending to zero) is enough to prove convergence. Also, using well-known results, you can shorten it up a bit, because once you get both cosines as above, and see that they're essentially the same, you can do the quadratic-equation-like substitution to render them both as particular cases of improper Fresnel integrals, i.e, integral between 0 and infinite of sin(x2) or cos(x2), both of which have very well-known, closed form values. "I accept the compliment, and certainly assume that your use of "didactic" (which you've done before) refers to definition 1 from www.dictionary.com 'Intended to instruct' and not definition 3: :-) 'Inclined to teach or moralize excessively' Certainly. The English term 'didactic' translates in my native Spanish as 'didactico', both of them with a Latin origin, and in Spanish your definition 1 is the only meaning, the Spanish term has no meaning equivalent to your definition 3 above. For that matter, I wasn't aware of that particular meaning in English, one always learns something new each and every day. "Sometimes, an enlightened transformation of the problem is required in order to achieve success." Sometimes even that won't save you. As an example, see if you can compute this simple integral, how many digits can you confidently say are correct in your result, and further, if you do recognize the exact, closed form: ``` / Inf | FRAC(x) I7 = | ------- .dx | x2 / 1 ``` where FRAC(x) is, of course, the fractional part of x, i.e., FRAC(3.14) = 0.14. I think that no "enlightened" transformation will save you now. But if you do succeed, then by all means: Give the man a cigar ! :-) Best regards from V.

 Re: (More) So your HP can INTEGRATE ...Message #10 Posted by Kiyoshi Akima on 7 July 2006, 5:23 p.m.,in response to message #9 by Valentin Albillo I7 ~ 0.4228 +/- 0.00005 This may not be the best way, but here goes... ``` /Inf I7 = | FRAC(x)/x2 dx /1 /2 /3 = | (x-1)/x2 dx + | (x-2)/x2 dx + ... /1 /2 /2 /2 /3 /3 = | x/x2 dx - | 1/x2 dx + | x/x2 dx - | 2/x2 dx + ... /1 /1 /2 /2 ``` Summing the x/x2 terms gives ln(Inf)-ln(1)=ln(Inf). Deferring the infinite limit for now, we go on to the other terms. ``` /k+1 | k/x2 dx = 1/(k+1) /k ``` The sum of these also diverge. Now for the upper limit. Since FRAC(x)<=1, the integrand is bounded by 1/x2. Actually, since FRAC(x) averages 0.5 and the larger parts get divided by larger quotients, it's bounded by 0.5/x2. (This also gives us 0.5 as a first, rough approximation.) For whatever accuracy epsilon we desire, we merely need to select U such that ``` /Inf | 0.5/x2 dx < epsilon /U |Inf 0.5/x | < epsilon |U 0.5/U < epsilon 0.5/epsilon < U ``` For four decimal digits, since we already have a rough estimate between 0.1 and 1.0, we select epsilon=0.00005 so that U>10000. Thus: ``` 10001 I7 ~ ln(10002) - SIGMA 1/(i+1) i=1 ``` which produces the value given at the beginning. So we've evaluated a bounded integral by splitting it into two unbounded expressions, evaluating each of them separately, and adding them together. Since U is an inverse linear function of the error tolerance, to get one more decimal digit we need to sum ten times as many terms (and adjust the argument of the ln()). Does that qualify as an "enlightened" transformation? I suppose I should add that the closed form for the result is 1-gamma, where gamma is Euler's constant. Edited: 7 July 2006, 6:23 p.m.

 Re: (More) So your HP can INTEGRATE ...Message #11 Posted by Rodger Rosenbaum on 7 July 2006, 7:00 p.m.,in response to message #10 by Kiyoshi Akima OK, now how about this one: ``` /Inf I8 = | F(x) dx /0 ``` where F(x) is defined as 0 if x is irrational and 1 if x is rational.

 Re: (More) So your HP can INTEGRATE ...Message #12 Posted by Valentin Albillo on 7 July 2006, 10:10 p.m.,in response to message #10 by Kiyoshi Akima Hi, Kiyoshi Akima: Kiyoshi Akima posted: "Does that qualify as an "enlightened" transformation?" Not quite, but your concoction's pretty funny, if not exactly rigorous. The closed-form value, 1-Gamma, is actually correct. Nevertheless, this integral isn't particularly difficult at all to numerically evaluate, it's just weird so to say, and the fact that it does have such a nice closed-form value is frankly unexpected. This simple HP-71B code computes its numerical value to 5 correct digits: ``` 10 S=0 @ FOR I=1 TO 6 @ S=S+INTEGRAL(10^(I-1),10^I,1E-12,FP(IVAR)/IVAR/IVAR) 20 NEXT I @ FIX 5 @ DISP S >RUN 0.42278 ``` so no sweat with this one. Best regards from V.

 #3Message #13 Posted by Crawl on 5 July 2006, 10:54 a.m.,in response to message #1 by Valentin Albillo So far, I've only had time to do number 3 (and verify that 1 and 2 give the results already given, without extra prompting). For #3, I first integrated the integral as given from 0 to 1. I don't remember my rational for that, but it shouldn't hurt. The result was 0.312942299.. For the remainder, I did the change of variable u = 1/(8-4x) This gives the integral ``` /-\ inf | | (2-1/(4u))^(pi/4) * sin(pi*u)/(4*u^2) du \-/ 1/4 ``` You can see now why the integral converges: It's due to the 1/u^2 term. We can put a rigorous limit on this integral. The first term is always less than or equal to 2^(pi/4). The next term is always (in absolute magnitude) less than 1. So, let's suppose, instead of integrating from 1/4 to infinity, we went from 1/4 to a. Then there error would be less than ``` /-\ inf | 2^(pi/4) | du/(4*u^2) \-/ a ``` = 2^(pi/4)/(4a) Let's say we want this error to be less than 0.00005. Then a would have to be greater than 8617.8396. For the heck of it, I rounded to 10,000. The previous integral, between 1/4 and 10,000, is 0.698310297... So, the final answer is 1.01125

 #5 Message #14 Posted by Crawl on 7 July 2006, 12:13 a.m.,in response to message #13 by Crawl The first step is to substitute u = 1/x, giving ``` /-\ Inf | sin u / u du | \_/ 1 ``` There's an "easy" way of doing this one if you happen to know that that integral from 0 to infinity is pi/2. Then you could just integrate that function from 0 to 1 (giving 0.946083) (which isn't very difficult to get with a calculator), and subtract it from pi/2 (giving 0.624713). Perhaps another way to do it would be by identifying that integral as being similar to the sum 1 - 1/3 + 1/5 - 1/7 + 1/9 - ... = pi/4 and using the Euler Maclaurin formula. We'll use ``` /-\ Inf | sin u / u du == | \_/ 1 /-\ a*pi/2 | sin u / u du + | \_/ 1 /-\ Inf | sin (pi * u / 2) / u du . | \_/ a ``` The first integral we'll do with straightforward numeric techniques. The second equals ```--- Inf \ / sin(pi*x/2) / x --- x=a + 1 ``` + (1/2) sin(pi*a/2)/a - (1/12)*(sin(a*pi/2)/a^2-a*pi/2*cos(a*pi/2))/a^2 We could add more correction factors to make the error smaller, but they become tedious to calculate. The infinite sum will be calculated not by actually calculating THAT sum, but rather by subtracting the missing terms from pi/4. As is, we'll use a value of 86 for "a". The first integral is 0.632115. The sum is -0.005813 and the next two correction terms are: 0 -(1/12)*(43*pi)/(86)^2 = -.00152208.. The error is guaranteed to be less than that last value. In fact, it's even better; all of the first 4 digits are correct. Putting those together gives 0.62477 Edited: 7 July 2006, 12:15 a.m.

 Re: So your HP can INTEGRATE ...Message #15 Posted by Kiyoshi Akima on 5 July 2006, 12:58 p.m.,in response to message #1 by Valentin Albillo Here are my initial results to four decimal digits. I won't explain how I got them yet in order to avoid leading other people down the wrong path :-) Everything was done on an HP 48GX and a couple of sheets of paper. I1 ~ -0.9641 I2 ~ -0.1976 I3 ~ 1.011 I4 ~ ???? I5 ~ 0.6247 I6 ~ ???? Now to go to work on the other two (and to verify the four results)...

 Re: So your HP can INTEGRATE ...Message #16 Posted by Valentin Albillo on 6 July 2006, 4:21 a.m.,in response to message #15 by Kiyoshi Akima Hi, Kiyoshi Akima: Your result #1 surely has a typo. Thanks for your interest and keep on trying, 24 hours still left. Best regards from V.

 Re: So your HP can INTEGRATE ...Message #17 Posted by Kiyoshi Akima on 6 July 2006, 1:45 p.m.,in response to message #16 by Valentin Albillo Right. A transposition error. I1 ~ 0.9461 and not 0.9641 as originally given. I can type on a calculator. It's these QWERTY keyboard that give me trouble :-)

 Re: So your HP can INTEGRATE ...Message #18 Posted by Kiyoshi Akima on 6 July 2006, 3:10 p.m.,in response to message #17 by Kiyoshi Akima To save space I won't repeat the methods used by other people, but only describe what I did differently. I1 ~ -0.9461 to four digits A singularity at the left endopoint. This can be done by almost any open method that avoids the singularity. The built-in integrator on the HP 48 in FIX 5 uses 511 function evaluations. A good adaptive routine may do better. But we can do even better than that. By using the series expansion for cos(x), the integral becomes ``` /1 x2 x4 x6 | (1 - -- + -- - -- + ... ) ln(x) dx /0 2! 4! 6! ``` Using the elementary integral ``` /1 x(i+1) 1 |1 1 | xi ln(x) dx = ------- (ln(x) - ---) | = - ------- /0 i+1 i+1 |0 (i+1)2 ``` the integral is replaced by the series ``` 1 1 1 1 -1 + ----- - ----- + ----- - ----- + ... 32 2! 52 4! 72 6! 92 8! ``` the first four terms of which sum to the value given above with an error less than the magnitude of the fifth next term ~ 3E-7. I2 ~ -0.1976 to four digits No singularities but is highly oscillatory, best handled by methods that don't sample the interval uniformly. The HP 48 integrator in SCI 5 gives the above result with 127 function evaluations. However, this is tailor-made for Filon's method for ``` /b | g(x) sin(kx) dx /a ``` which gives the same result with only 16 evaluations of g(x).

 Re: So your HP can INTEGRATE ...Message #19 Posted by Rodger Rosenbaum on 6 July 2006, 11:31 a.m.,in response to message #1 by Valentin Albillo Numbers 1, 2, 4 and 5 have closed form solutions, found in Gradshteyn and Ryzhik and other large tables of integrals. Number 6 is, of course, the one you discussed a while ago. Number 3 is therefore the most interesting, requiring a contour integration for the best result. ```No. 1 -.946083 No. 2 -.19763 No. 3 1.01123909053 ``` Number 3 was done on my HP71. I had to make the parabolic contour be further from the real axis to get the last two digits right with a tolerance of 1E-9 in the INTEGRAL command. It took 915 seconds on my sped-up HP71. ```No. 4 .4916997 No. 5 .62471 No. 6 .323367+, the same result obtained in your posting of some months ago. ```

 Difficult Integrals: My Original Solutions and CommentsMessage #20 Posted by Valentin Albillo on 11 July 2006, 5:18 p.m.,in response to message #1 by Valentin Albillo Hi all: Thanks for all your interesting inputs, whether actual results or comments, to my question about the actual efectiveness of HP calcs' built-in methods of integration when applied to these not-so-easy definite integrals. I'll give now correct numerical results for all six integrals, and will discuss what happens when you try to compute them blindly, and what can you do when this naive approach doesn't deliver. All given code is for the HP-71B, which uses much the same Romberg's quadrature method as former (HP-34C) and later models, so your results should be about the same no matter what model you do use. ``` / Inf I1 = | sin(x)*sin(x^2) .dx = 0.491695777984 / 0 ``` This is a very difficult integral to compute numerically to any decent accuracy, because the integrand does not tend to zero as the limit of integration tends to infinity but oscillates instead with increasing frequency, so at first sight it might be questionable whether it actually converges or not (which is does). The naive approach, using 10^6 for infinity, takes a looong time, and blatantly fails: ``` INTEGRAL(0,1E6,1E-12,SIN(IVAR)*SIN(IVAR*IVAR)) -.356719258 ``` Insisting on blindly going on by using a plain 4096-point Gaussian quadrature applied to the interval (0, 10^6), subdivided in 16384 subintervals does succeed in producing four correct digits, 0.4917, but it takes a long time as well and requires a fast PC, so it seems plain that some math thinking will be needed if we're going to go much farther (or even that far) while using a much slower handheld. Fortunately, there's a nifty closed form for a generalized family of similar integrals, this being but one of the simplest cases, namely: ``` / A / A I1 = (Cos(A)* | Cos(x)/Sqrt(x).dx + Sin(A)* | Sin(x)/Sqrt(x).dx )/2 / 0 / 0 ``` where A=1/4 and both non-elementary integrals are particular cases of Fresnel functions. Using the above expression, the integral can be easily and accurately computed right from the keyboard (no programming needed) as follows: ``` >A=1/4 @ P=1E-12 >(COS(A)*INTEGRAL(0,A,P,COS(IVAR)/SQR(IVAR))+SIN(A)*INTEGRAL(0,A,P,SIN(IVAR)/SQR(IVAR)))/2 .491695777984 ``` ``` / 1 I2 = | cos(ln(x)/x)/x .dx = 0.323367431678 / 0 ``` This one also takes quite long, if using the naive approach, and the initial result is not exactly promising: ``` >INTEGRAL(0,1,1E-12,COS(LN(IVAR)/IVAR)/IVAR) .169114784125 ``` Also, this time even a ponderous 8192-point Gaussian formula applied to 4096 subintervals does utterly fail, producing the completely wrong value 1.7656... In this case, the easiest way to numerically compute this integral is by integrating in the complex plane, as discussed previously in my post of the reference: ``` "Re: Integration Times "Old" 33s vs "New" 33s Message #5 Posted by Valentin Albillo on 13 Dec 2005, 4:28 a.m., in response to message #1 by John Smitherman" ``` Self-quoting from that post: You just need to remember that ``` e^ix = cos(x) + i*sin(x) ``` and the original integral becomes: ``` Integral = Integral(0, 1, RealPart(e^(i*log(x)/x)/x)) = Integral(0, 1, RealPart(x^(i/x-1))) ``` which can readily be integrated by considering it a line integral in the complex plane along some integration path. A simple parabolic path will do: ``` z = t + t*(1-t)*i ``` which takes us far from the oscillations of the real line. Using this path, your faithful HP-71B can compute the integral to full 12-digit precision in a few minutes, as follows: ``` 10 DESTROY ALL @ COMPLEX Z @ SFLAG -1 @ DISP INTEGRAL(0,1,1E-9,FNY(IVAR)) 20 DEF FNY(T) @ Z=(T,T-T*T) @ FNY=REPT(Z^((0,1)/Z-1)*(1,1-2*T)) >RUN .323367431678 ``` which, as stated, is the correct result rounded to the 12 decimal digits shown. This technique of computing difficult real-line integrals by taking them to the complex plane is actually very useful in many situations, and is discussed at length (with an example) in the HP-15C Advanced Functions Handbook. ``` / 1 I3 = | sin(1/x)/x .dx = 0.624713256428 / 0 ``` Once again, the naive approach takes quite long, and miserably fails: ``` >INTEGRAL(0,1,1E-12,SIN(1/IVAR)/IVAR) -1.6116828454 ``` Also, the 8192-point Gaussian formula, applied to 1024 subintervals, fares no better, producing utter garbage: -2.342289... The way to go is by using 1/x = t, which, together with the well-known elementary result for the integral between 0 and Infinity of Sin(t)/t .dt, results in the following expression: ``` / Inf / Inf / 1 / 1 I3 = | sin(t)/t .dt = | sin(t)/t .dt - | sin(t)/t .dt = PI/2 - | sin(t)/t .dt / 1 / 0 / 0 / 0 ``` which can be computed very quickly, to provide 12-digit accuracy: ``` >PI/2-INTEGRAL(0,1,1E-12,SIN(IVAR)/IVAR) .624713256433 ``` ``` / 1 I4 = | cos(x)*ln(x) .dx = -0.946083070367 / 0 ``` For one, the naive approach does work, even if it takes quite long: ``` >INTEGRAL(0,1,1E-12,COS(IVAR)*LN(IVAR)) -.946083070191 ``` Using just one subinterval gets us 10 more-or-less correct digits. We can get full 12-digit accuracy by subdividing further. ``` / 2*Pi I5 = | ln(1+x)*sin(10*x) .dx = -0.197626807719 / 0 ``` The easiest of the pack, the naive approach works too, and very quickly: ``` >INTEGRAL(0,2*PI,1E-12,LN(1+IVAR)*SIN(10*IVAR)) -.197626807716 ``` and produces nearly 12-digit accuracy. ``` / 2 I6 = | x^(Pi/4)*sin(Pi/4/(2-x) .dx = 1.01123909053 / 0 ``` The naive approach also works, but with reduced accuracy and taking very long: ``` >INTEGRAL(0,2,1E-12,IVAR^(PI/4)*SIN(PI/4/(2-IVAR))) 1.0112556592 ``` In search of increased accuracy, we can make it last twice as long by using two subintervals: ``` >P=PI/4 @ E=1E-12 >INTEGRAL(0,1,E,IVAR^P*SIN(P/(2-IVAR)))+INTEGRAL(1,2,E,IVAR^P*SIN(P/(2-IVAR))) 1.01119507067 ``` or even ten times longer: ``` >S=0 @ FOR X=0 TO 1.9 STEP .1 @ S=S+INTEGRAL(X,X+.1,E,IVAR^P*SIN(P/(2-IVAR))) @ NEXT X @ S 1.01122722759 ``` But actually, this is yet another integral which is much easier to compute in the complex plane: ``` 10 DESTROY ALL @ COMPLEX Z @ SFLAG -1 @ DISP INTEGRAL(0,2,1E-10,FNY(IVAR)) 20 DEF FNY(T) @ Z=(T,T+T-T*T) @ FNY=REPT((0,-1)*EXP(PI/4*(LOG(Z)+(0,1)/(2-Z)))*(1,2-2*T)) >RUN 1.01123909053 ``` Best regards from V.

 Re: Difficult Integrals: My Original Solutions and CommentsMessage #21 Posted by Walter B on 11 July 2006, 6:12 p.m.,in response to message #20 by Valentin Albillo Hi Valentin, are you teaching math? Just wonder... I faintly remember solving integrals in the complex plane way back at the university. Never needed it thereafter anymore - carefully chose my jobs this way ;) My congratulations to your work!

 Re: Difficult Integrals: My Original Solutions and CommentsMessage #22 Posted by Valentin Albillo on 12 July 2006, 4:24 a.m.,in response to message #21 by Walter B Hi, Walter: Walter posted: "Hi Valentin, are you teaching math? Just wonder..." Wish I would ! ... :-( ... But no, my career has been mainly centered in systems engineering and things computational. Nevertheless I was fond of mathematics long before I knew what a slide rule was, much less a calculator or computer. Since I was 8 years old or so I was already doing math in my spare time and liking it a lot. This has continued throughout my adult life, and it's my life-long passion. When I was 17, I actually wanted to go for the equivalent of a PhD in Math, but was forced to chose a different career instead, which I've always regretted. As for teaching math, I've certainly imparted some classes here and there to individuals, with surprising success as everyone enjoyed them a lot and said individuals went on to get high grades in math where previously they had miserably failed. I credit that success to my loving the subject matters and most importantly, being capable of teaching it in such a way that even the most boring subject seemed enjoyable to the student. Just for instance, I would chose interesting, unusual examples that had particularly likeable or unexpected results, then discussed why they were so, thus enhancing the subject matter's understanding, etc. "My congratulations to your work! " Thank you very much, your kind comment is much appreciated and I'm truly glad you've enjoyed my humble efforts. Best regards from V.

 Re: Difficult Integrals: My Original Solutions and CommentsMessage #23 Posted by Karl Schneider on 12 July 2006, 2:07 a.m.,in response to message #20 by Valentin Albillo Hi, Valentin -- Thank you for the numerical answers and informative discussion -- and particularly for the HP-71B code. I'll take some time to try your solutions. -- KS

 Thanks a lot, Karl ! [No Text]Message #24 Posted by Valentin Albillo on 12 July 2006, 4:27 a.m.,in response to message #23 by Karl Schneider Best regards from V. Go back to the main exhibit hall