Improper (difficult) integrals at SwissMicros Calculator Forum

[Valentin Albillo]

      
Hi, all,

Some six of my difficult improper integrals that I recently posted to Namir's thread here about 16-point Gaussian integration, have also been posted a few days ago (not by me) at the SwissMicros Calculator Forum, where people there have struggled to compute them accurately (if at all ! .) Some comments (I have no account there):

• Not all of them have been accurately computed, though some members (e.g. henrik) made heroic efforts; for #5 he provided an accurate plot, which looks terrific, and then he went on to provide 1,000 digits or so using Mathematica's NIntegrate, but regrettably it barely has * 3 * correct digits so all 997 subsequent ones are but garbage !. So much for Mathematica ...
        
• The link above points to the OP's message with his results (and some of mine,) but there are 4 pages in all with posts by some members and even ex-members of this (MoHPC) forum, several dozen messages in all.
  

You might want to have a look at them and see whether your favorite HP calcs can do better than SwissMicros models ... or even Mathematica for that matter !  

Regards.
V.

[peacecalc]

Hello Valentin,

I tried to integrate your fifth integral from your test-set numericaly:

\[ I = \int_{-1}^{1}\exp\left(x + \sin\left(\exp\left(\exp\left(\exp\left(x + \frac{1}{3}\right)\right)\right)\right)\right)dx \] .

The argument from the sin has a range from 5 to 1.9E19 , so its a wild oscillating story. I tried to find an integrand which has a similar value like this definite integral and I find something like this:

\[ I = \int_{-1}^{1}\exp\left(x\cdot\left(\frac{1+\sqrt{5}}{2}\right)\right)dx \] .

Why I search for this? The reason is simple: I programmed a simple midrule integrator. That program gives the first 4 valid digits of your integral above. The prime integration program gives a value over 3 and the simple program of mine does it better: 2.994... (but with 10^6 calls of function, but prime needs only 4 minutes for this). So the midpoint rule is "blind" for the behavier of the function values. If the oscillating is very wild, the rule will calculate a statistical mean value which is near by the true value.
That is only my assumption. I get the idea, the integrand function behave in a way like the function:
\[ \exp(a\cdot x) \] (the graph shows it clearly, it is an exponential form superimposed with a sinus) and a is the golden ratio (maybe?, it is a guess).

My questions are: from where you know so many valid digits for your test-integral?
and: in some posts you wrote about your own in basic programmed integrator for your HP-71. Do you like to share the code or if you have had already shared it, can you tell us the link?

Thank you very much for test suite in case of numeric integration, I learned a lot, even I failed for getting much more valid digits for your integral number 5. I like this experimental way of doing mathematics, because my mathematical skills, aren't that good, they should be.

[Valentin Albillo]

.
Hi, peacecalc,

(11-25-2023 08:29 PM)peacecalc Wrote:  I tried to integrate your fifth integral from your test-set numericaly:
\[ I = \int_{-1}^{1}\exp\left(x + \sin\left(\exp\left(\exp\left(\exp\left(x + \frac{1}{3}\right)\right)\right)\right)\right)dx \]

Good. The more people interested, the merrier.

Quote:The argument from the sin has a range from 5 to 1.9E19 , so its a wild oscillating story.

Yes, it is. Quite wildly.

        

Quote:[...] the simple program of mine does it better: 2.994... (but with 10^6 calls of function, but prime needs only 4 minutes for this).

Getting four digits (save 2 ulp) is no mean feat, congratulations. And speed is one of the redeeming hardware features the to.. er, the Prime has. Did you use its native programming language or its limited Python implementation, which I'm told is much faster ? Nevertheless, some 4,200 f(x) evals per sec. is no mean feat.

Quote:So the midpoint rule is "blind" for the behavier of the function values. If the oscillating is very wild, the rule will calculate a statistical mean value which is near by the true value. That is only my assumption.

Well, as seen in the plot above, the oscillations of f(x) are so fast that the graph becomes an impenetrably dense yellow blur, but if you subtract the mean value of the oscillations over one period of sin(x) (plotted in light blue) from f(x) itself (plotted in yellow), the positive and negative contributions to the area tend to cancel each other out and so contribute very little in all.

Thus, we can obtain a decently accurate value (about 3 correct digits) of I5 by simply computing the relevant mean value (a constant which we'll call C) over a full [0..2\(\pi\)] period, then we just compute the integral of the exponential function times C like this (because exp(x+sin(...)) = exp(x) * exp(sin(...)) ):
        
where both integrals are very easily computed (and very quickly at that) using e.g. the HP-71B:

>FIX 2 @ K=1E-3 @ INTEGRAL(0,2*PI,K,EXP(SIN(IVAR)))/(2*PI)*INTEGRAL(-1,1,K,EXP(IVAR))

      2.98

or the second integral can be simply replaced by its obvious value:

>FIX 2 @ INTEGRAL(0,2*PI,1E-3,EXP(SIN(IVAR)))/(2*PI)*(EXP(1)-EXP(-1))

      2.98

Quote:My questions are: from where you know so many valid digits for your test-integral? [...]

I got them from two different reputed sources using two different programs, whose results perfectly agree so it's quite safe to assume they're correct.

Quote:[...] in some posts you wrote about your own in basic programmed integrator for your HP-71. Do you like to share the code or if you have had already shared it, can you tell us the link?

I'd like to but my ancient integration program, which I initially wrote for the HP-71B several decades ago and then also ported it to UBASIC (thus providing me with up to 1,200-digit results and several orders of magnitude speedup,) was never completely finalized to my satisfaction, some improvements were still pending.

A few years ago I set up to complete its implementation and once done I was about to create an article in the format and style of my many articles already published, intending to publish it as well, but when I participated in some thread about such matters and announced its next availabilty, I was subjected out of the blue to unpolite attempts to ruin my announcement by outguessing my implementation, belittling previous contributions by me, and finally subjecting me to thinly-veiled attacks and uncalled-for disrespect, which was the last thing I expected as a reply to my highly technical posts in that thread. Thus, I ditched it all for good, dismissing any attempts at publication and left the forum for a time out of sheer disgust.

The bottom line is that I've misplaced its source code since then and I only have the carefully-crafted ancient tests nowadays, but in general it was as fast as the assembler INTEGRAL keyword for normal integrals and many times faster for difficult integrals. A pity. It was the second time something like this happened to me, with identical results.

Quote:Thank you very much for test suite in case of numeric integration, I learned a lot, even I failed for getting much more valid digits for your integral number 5. I like this experimental way of doing mathematics, because my mathematical skills, aren't that good, they should be.

You're welcome, that's my main goal here: to provide entertainment and knowledge whenever possible via my posts, challenges, articles, the works !

Thanks for your interest and regards.
V.

[Ajaja]

(11-26-2023 03:55 AM)Valentin Albillo Wrote:  Yes, it is. Quite wildly.

        

What about arc length of the curve?

[peacecalc]

Hello Valentin,

thank you very much for your quick and comprehensive answer. I didn't expect
this. I feel esteemed.

off topic:
I regret and feel sorry, that you had made bad experiences with some
members of this forum. Using fora is a tricky thing. As a member
in some other fora I saw that even one forum died because of a bad ping-pong
post game between only few members.

ontopic:
About your question:
I used the native programming language of the hp prime (PPL).
Maybe I give the builtin Python a chance, but for getting a
higher precision the number of calls have to raise with factor 10,
that means, instead of 4, I have to wait 40 minutes.
So Python has to be much more faster then PPL.

I carried out the estimated Integral:

\[ I = \int_{-1}^{1}\exp\left(x\cdot\left(\frac{1+\sqrt{5}}{2}\right)\right)dx = (\sqrt{5}-1)*\sinh\left(\frac{\sqrt{5}+1}{2}\right) \approx 2.9943\]

thats not bad for such an easy estimation.

Because there only a few valid digits,
it is clear, that I haven't found a secrect relation, it is only
a not bad guess.

@Ajaja: Just carry out the integral

\[ s = \int_{-1}^{1}\sqrt{\left(1+f'(x)^2\right)}dx =
\int_{-1}^{1} \sqrt{\left(1+\left((1+e^{x+\frac{1}{3}} e^{e^{x+\frac{1}{3}}} e^{e^{e^{x+\frac{1}{3}}}}\cos\left(e^{e^{e^{x+\frac{1}{3}}}}\right)) e^{x+\sin\left(e^{e^{e^{x+\frac{1}{3}}}}\right)} \right)^2\right)}dx \]

Maybe a substitution makes it easier..., but it looks ugly, very ugly...
By the way, for writing down the integral in my poste, I used the command "latex" of hp prime. I try to make it smoother...

\[ \int_{-1}^{1}\sqrt{1+f'(x)^2}dx =
\int_{-1}^{1} \sqrt{1+\left(\left(1+ \exp\left(x+\frac{1}{3}\right) \exp\left(\exp\left(x+\frac{1}{3}\right)\right)\exp\left(\exp\left(\exp\left(x+ \frac{1}{3}\right)\right)\right)
\cos\left(\exp\left(\exp\left(\exp\left(x+ \frac{1}{3}\right)\right)\right)\right) \right)
\exp\left(x+\sin\left(\exp\left(\exp\left(\exp\left(x+\frac{1}{3}
\right)\right)\right)\right)\right)\right)^2}dx \]

However I do, the integral is still ugly. Editing was necessary because of false brackets setting and typo at estimated value.

[Ajaja]

@peacecalc
I just wonder if it's even possible to calculate it using a programmable calculator. Looking at the graph I see how the numerical integration can produce more or less accurate results for the integral, but for the arc lenght it's seems impossible.

[lrdheat]

Even breaking it up into 4 integrals, from 10^-18 to 10^-17 + from 10^-17 to 10^-12 + from 10^-12 to 10^-6 + from 10^-6 to 1 yields 33.93+ on my $20 CASIO 991CW. Fascinating, illuminating integration problem for numerical integration capabilities!

[peacecalc]

Hello,

@Ajaja, yes the area integral is tricky, because the adaptive quadratur methods fails in the Intervall
from lets say x in 0.5 to 1. And the integral for the arc length will be for such methods good possible fom x in -1 to 0, but for x > 0 the quadratur algorithmus runs into trouble.

@Irdhat in principle it is a good idea to integrate such monsters in smaller areas. But I didn't understand , why you choose the greatest intervall from 10-6 to 1. In my view, je nearer you come to 1 you have to make smaller intervalls, for example:
I -1 to 0
II 0 to 0.5
III 0.5 to 0.7
IV 0.7 to 0.8
VI 0.8 to 0.85
VII 0.85 to 0.88
...
and so on, that can be managed with a program. For the arc length im not shure wether 34 is a correct estimation, because the function is oszillating with a very high frequency (which raises!) between (e^x)*(e^{-1}) and (e^x)*(e^{1}) for example x =0.8 you get a difference value 5.231. For x = 0.8 you get a very great gradient (near vertical plus or minus), so the arc length has to be very long, I assume. But I give hp prime a chance to manage it (with the method above).

[lrdheat]

Interestingly, integrating from 10^-18 to .99 + from .99 to 1-10^-6 + from 1-10^-6 to 1-10^-12 + from 1-10^-12 to 1-10^-15 + from 1-10^-15 to 1-10^-18 produced the exact same answer! Integration from 10^-18 to .99 produces 33.92+, or about .1 lower than the earlier result. The integrations above .99 only contribute another ~.1 The big contribution to the real answer of 100 comes from the very first part of the interval. My CAS unit shows 81+ of the 100 comes from the the first 1 millionth part (from 0 to 10^-6)! My CASIO 991CW correctly showed the result from .99 to 1.

[peacecalc]

Hello Irdhat,

it is no offence.

I'm confused, because I'm not shure that we talk about the same integral or did you use a transformation, which changes the startpoint of integration? If second case is right, can you show how you made the transformation?

Thank you in advance.

[lrdheat]

Your suspicion is correct! I am sorry for the confusion. I was computing x^(-.99) from 0 to 1. My apologies! This equals 100, but numerical methodologies are quite challenged by this integral!

[peacecalc]

Hello Irdheat,

thank you for clarification. The arc-length of the function handled in this thread is about 4E19 for x in [-1; 1] and the main part comes from [0,95 ; 1] the number of oszillations is in the same magnitude (3E18). If one has a look to the derivative, it is seen, that the magnitude of the arclength of the function is not only oszillating with the sin, it is superimposed by a cos, too.

Every physicist or engineer could be glad that his problem needn't such an integral for solution.