So your HP can INTEGRATE ...



#20

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.

#21

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


#22

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.

#23

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.


#24

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.


#25

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


#26

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.


#27

Hello, Valentin --

On I3:

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 I5:

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 arms-limits 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 ever-skinnier. 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 closed-form solution, according to Rodger -- that should suffice to prove the result.

On I6:

I've already found your post from December 2005, and will examine it within the context of material presented in the HP-15C 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 numerical-integration problem they might be posed, even if some performance-enhancing 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 first-year calculus techniques, but the complex-plane methods are certainly worth a closer look for my own edification.

Thanks and best regards,

-- KS


Edited: 7 July 2006, 1:53 a.m.


#28

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:

    1. Using elementary trigonometric identities, express the product as a sum of two cosines.

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

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

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

#29

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.


#30

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.

#31

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

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


#33

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.

#34

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


#35

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.


#36

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


#37

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).
#38

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.

#39

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.

#40

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!


#41

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.

#42

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


#43

Best regards from V.


Possibly Related Threads...
Thread Author Replies Views Last Post
  HP Prime: Integrate Won't Integrate James Williams 17 3,398 11-15-2013, 01:41 PM
Last Post: CR Haeger
  PSE bug also in Integrate and Solve Namir 0 564 09-18-2011, 04:05 PM
Last Post: Namir
  15C LE expanded a bit on SOLVE & INTEGRATE? Namir 6 1,414 09-13-2011, 01:43 AM
Last Post: Marcus von Cube, Germany
  re-release of the voyager with matrix, solve, integrate & complex numbers exschr 17 3,288 06-08-2011, 02:00 PM
Last Post: Dave Shaffer (Arizona)
  Solve/Integrate Algorithms Patrick 3 973 02-25-2003, 09:35 AM
Last Post: Ex-PPC

Forum Jump: