▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
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
HP11C, HP42S, HP41C's Advantage ROM and HP71B's Math ROM,
to name a few models and plugins.
As a mention was made to a certain difficulttonumericallycompute
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
I_{1} =  cos(x)*ln(x) .dx
/0
/2*Pi
I_{2} =  ln(1+x)*sin(10x) .dx
/0
/2
I_{3} =  x^{Pi/4}*sin(Pi/(84x)) .dx
/0
/Inf
I_{4} =  sin(x)*sin(x^{2}) .dx
/0
/1
I_{5} =  sin(1/x)/x .dx
/0
/1
I_{6} =  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.
▼
Posts: 177
Threads: 12
Joined: Jan 1970
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
▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, John:
John posted:
"1. 0.9461
2. 0.1976
3. 1.0108"
Right, right, and notsoright.
"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"
"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.
Posts: 1,792
Threads: 62
Joined: Jan 2005
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:
I _{1}: 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 HP48G with FIX 6 yielded
I_{1} = 0.94608307068
Integrating the original function directly with FIX 6 yielded
I_{1} = 0.94608289406
with more computational time.
I_{2}: f(x) = ln(1+x)*sin(10x) from 0 to 2*pi.
This integrated directly and quickly on an HP42S to
I_{2} = 0.197627
Too easy. Did I miss something?
I_{5}: 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 = u^{2} 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 = 10^{7} and FIX 6,
I_{5} = 0.62479588
 KS
Edited: 4 July 2006, 3:56 a.m.
▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
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""
"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 HP25. 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.
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
One more result, and some discussion:
I_{3}:
f(x) = x^{pi/4}* sin(pi/(84x)); 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 highfrequency oscillations, and took quite a while to run on an HP49G; 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).
I_{4}:
f(x) = sin(x)*sin(x^{2}); 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(x^{2}) + 2*integral (x*cos(x)*cos(x^{2}).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.
I_{6}:
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
▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
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 HP71B (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 builtin "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 HP48, HP49, or HP71B/Math ROM."
Back in that thread I demonstrated how you could easily compute it to full accuracy and fairly quickly on an HP71B (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.
▼
Posts: 1,792
Threads: 62
Joined: Jan 2005
Hello, Valentin 
On I_{3}:
there were two blatant errors. Fatigue led to several careless mistakes, and I corrected them. However, my answer appears to be incorrect by one ulp (unit in the last place?). I suppose that FIX 7 would have done the trick, but how long should one wait for an answer, and still not know for sure whether it is correct?
On I_{5}:
I said,
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.
You said,
Trust me: it does.
This brings to mind a statement regarding armslimits treaty compliance, made in the 1980's by former US President Ronald Reagan (if I'm not mistaken):
"Trust, but verify."
There is no philosophical message in my quoting of that statement, but I did find the contradiction a bit amusing...
The frequency of the integrand becomes higher, so the slivers of area get everskinnier. 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. The stated integrand does not attenuate or approach any limiting values, and the upper limit is infinite. However, there is apparently a closedform solution, according to Rodger  that should suffice to prove the result.
On I_{6}:
I've already found your post from December 2005, and will examine it within the context of material presented in the HP15C AFH. No need for me to repeat the solution, until I fully understand the approach.
Thanks again for your unfailing interest in my 'challenges' and your always interesting (and didactic!) contributions.
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"
Your examples certainly prove the point that handheld calculators can't be expected to cope with every numericalintegration problem they might be posed, even if some performanceenhancing modifications are made. Sometimes, an enlightened transformation of the problem is required in order to achieve success.
My "arsenal" for these problems consists primarily of the standard firstyear calculus techniques, but the complexplane methods are certainly worth a closer look for my own edification.
Thanks and best regards,
 KS
Edited: 7 July 2006, 1:53 a.m.
▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
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 closedform solution, according to Rodger"
Yes, there's a generalized closedform 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 wellknown logarithmic series: 11/2+1/31/4+...).
This fact (i.e, alternating series and general term tending to zero) is enough to prove convergence.
Also, using wellknown 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 quadraticequationlike substitution to render them both as particular cases of improper Fresnel integrals, i.e, integral between 0 and infinite of sin(x^{2}) or cos(x^{2}), both of which have very wellknown, 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
 x^{2}
/ 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.
▼
Posts: 325
Threads: 18
Joined: Jul 2006
I_{7} ~ 0.4228 +/ 0.00005
This may not be the best way, but here goes...
/Inf
I_{7} =  FRAC(x)/x^{2} dx
/1
/2 /3
=  (x1)/x^{2} dx +  (x2)/x^{2} dx + ...
/1 /2
/2 /2 /3 /3
=  x/x^{2} dx   1/x^{2} dx +  x/x^{2} dx   2/x^{2} dx + ...
/1 /1 /2 /2
Summing the x/x ^{2} terms gives ln(Inf)ln(1)=ln(Inf). Deferring the infinite limit for now, we go on to the other terms.
/k+1
 k/x^{2} 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/x^{2}. Actually, since FRAC(x) averages 0.5 and the larger parts get divided by larger quotients, it's bounded by 0.5/x^{2}. (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/x^{2} 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
I_{7} ~ 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 1gamma, where gamma is Euler's constant.
Edited: 7 July 2006, 6:23 p.m.
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, Kiyoshi Akima:
Kiyoshi Akima posted:
"Does that qualify as an "enlightened" transformation?"
10 S=0 @ FOR I=1 TO 6 @ S=S+INTEGRAL(10^(I1),10^I,1E12,FP(IVAR)/IVAR/IVAR)
20 NEXT I @ FIX 5 @ DISP S
>RUN
0.42278
so no sweat with this one.
Best regards from V.
Posts: 306
Threads: 3
Joined: Sep 2009
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/(84x)
This gives the integral
/\ inf

 (21/(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
▼
Posts: 306
Threads: 3
Joined: Sep 2009
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^2a*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.
Posts: 325
Threads: 18
Joined: Jul 2006
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)...
▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
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.
▼
Posts: 325
Threads: 18
Joined: Jul 2006
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 :)
▼
Posts: 325
Threads: 18
Joined: Jul 2006
To save space I won't repeat the methods used by other people, but only describe what I did differently.
I_{1} ~ 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 builtin 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 x^{2} x^{4} x^{6}
 (1   +    + ... ) ln(x) dx
/0 2! 4! 6!
Using the elementary integral
/1 x^{(i+1)} 1 1 1
 x^{i} 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 +    +    + ...
3^{2} 2! 5^{2} 4! 7^{2} 6! 9^{2} 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 ~ 3E7.
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 tailormade for Filon's method for
/b
 g(x) sin(kx) dx
/a
which gives the same result with only 16 evaluations of g(x).
Posts: 305
Threads: 17
Joined: Jun 2007
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 1E9 in the INTEGRAL command. It took 915 seconds on my spedup HP71.
No. 4 .4916997
No. 5 .62471
No. 6 .323367+, the same result obtained in your posting of some months ago.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi all:
Thanks for all your interesting inputs, whether actual results or comments, to
my question about the actual efectiveness of HP calcs' builtin methods of integration
when applied to these notsoeasy 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 HP71B, which
uses much the same Romberg's quadrature method as former (HP34C) 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,1E12,SIN(IVAR)*SIN(IVAR*IVAR))
.356719258
Insisting on blindly going on by using a plain 4096point 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 nonelementary 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=1E12
>(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,1E12,COS(LN(IVAR)/IVAR)/IVAR)
.169114784125
Also, this time even a ponderous 8192point 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"
Selfquoting 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/x1)))
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*(1t)*i
which takes us far from the oscillations of the real line.
Using this path, your faithful HP71B can compute the integral to full
12digit precision in a few minutes, as follows:
10 DESTROY ALL @ COMPLEX Z @ SFLAG 1 @ DISP INTEGRAL(0,1,1E9,FNY(IVAR))
20 DEF FNY(T) @ Z=(T,TT*T) @ FNY=REPT(Z^((0,1)/Z1)*(1,12*T))
>RUN
.323367431678
which, as stated, is the correct result rounded to the 12 decimal digits shown.
This technique of computing difficult realline 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 HP15C 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,1E12,SIN(1/IVAR)/IVAR)
1.6116828454
Also, the 8192point 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 wellknown 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 12digit accuracy:
>PI/2INTEGRAL(0,1,1E12,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,1E12,COS(IVAR)*LN(IVAR))
.946083070191
Using just one subinterval gets us 10 moreorless correct digits.
We can get full 12digit 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,1E12,LN(1+IVAR)*SIN(10*IVAR))
.197626807716
and produces nearly 12digit accuracy.
/ 2
I6 =  x^(Pi/4)*sin(Pi/4/(2x) .dx = 1.01123909053
/ 0
The naive approach also works, but with reduced accuracy and
taking very long:
>INTEGRAL(0,2,1E12,IVAR^(PI/4)*SIN(PI/4/(2IVAR)))
1.0112556592
In search of increased accuracy, we can make it last twice as long by using two subintervals:
>P=PI/4 @ E=1E12
>INTEGRAL(0,1,E,IVAR^P*SIN(P/(2IVAR)))+INTEGRAL(1,2,E,IVAR^P*SIN(P/(2IVAR)))
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/(2IVAR))) @ 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,1E10,FNY(IVAR))
20 DEF FNY(T) @ Z=(T,T+TT*T) @ FNY=REPT((0,1)*EXP(PI/4*(LOG(Z)+(0,1)/(2Z)))*(1,22*T))
>RUN
1.01123909053
Best regards from V.
▼
Posts: 4,587
Threads: 105
Joined: Jul 2005
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!
▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
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 lifelong 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.
Posts: 1,792
Threads: 62
Joined: Jan 2005
Hi, Valentin 
Thank you for the numerical answers and informative discussion  and particularly for the HP71B code. I'll take some time to try your solutions.
 KS
▼
Posts: 1,755
Threads: 112
Joined: Jan 2005
