Trying to find Integration Challenge



#26

Hi all,

In recent months Valentin Albillo posted one of his famous challenges, this one regarding tough integrals, and there was a lot of discussion regarding which calc was faster and how to structure the problems in order for the algorithm to find the integral more efficiently, etc.

I didn't pay a lot of attention at the time because my attentions where focussed on issues other than quadrature. Now, the subject has drawn my attention.

The problem is, I can't seem to find the thread. Maybe I am using the search engine wrongly, but right now I can just seem to turn up very recent posts, and nothing older.

Can anyone who remembers this thread help me find the link?

Also, about the same time there was a lot of discussion about Gaussian quadrature vs. Romberg integration, etc. I think either Valentin or Gerson Barbosa headed up that discussion too.

Many thanks in advance,

Les


#27

Hi Les,

I prefer to use Google when searching for past threads. One of the threads you are looking for may be found with the following search key:

site:hpmuseum.org  Romberg + Simpson + Valentin + Gerson

This leads to http://www.hpmuseum.org/cgi-sys/cgiwrap/hpmuseum/archv016.cgi?read=94446, the only result found by Google. As you can see, site:hpmuseum.org makes the engine to search only in the museum site, the others words were included because I remember they were in that thread. The other thread you have mentioned may be found this way.

I hope this helps,

Regards,

Gerson.


Edited: 10 Dec 2006, 6:21 a.m.


#28

Gerson, than you so much! I was particularly interested in revisiting the Gaussian quadrature routine of Valentin. Most such routines that I have seen require precomputing and storing the weights and sample points, but Valentin's routine seems to compute these as required and does it with such little code--I really should write him and ask him to send me the 1980 article, since I can't discern exactly what happens from looking at the code alone. Back-engineering is not one of my strengths!

As for locating the Integration Challenge, my motivation is a little sinister, I am afraid. I few months ago, I slammed the TI83 Plus I had just acquired. Recently, as I find myself wanting to give TI products their due, I have gone back to the 83 and have tried to learned a little more about it.

In the process I have discovered the built-in numerical integration feature (the fnInt() function for those who know the calculator), and I must admit that I am really impressed. I have been working with a particular integral that takes a goodly time to complete to reasonable precision on the 15C, 42S, and 49G+, but seems to complete quickly on the lowly TI83.

So, I want to revisit the integration challenge and see how the TI83 handles the problems Valentin posed for the group. IIRC, I don't think this particular calc was mentioned, but then again it was a massively long thread and I could be misremembering.

When I find the thread I will let you know how I do.

Les


#29

Hi again Les,

I don't want to spoil your searching fun but I think this is one of threads:

So your HP can INTEGRATE ...

Please let us know how your TI83 handles those! :-)

Regards,

Gerson.

------------


It seems John Smitherman was able to solve them with the TI83+SE, except the last three.


Edited: 10 Dec 2006, 3:19 p.m.

#30

Les, here's a little information on Gaussian Quadrature. Usually it's called Gauss-Legendre Quadrature. Gauss showed that the optimal points (nodes) for evaluating the integrand on an interval are the zeros to a Legendre polynomial (easily found on the HP50 with the Legendre and Zeros commands). However, the zeros are on the interval (-1,1), so the integrating function needs to be translated so that it is centered at the origin, and then horizontally "compressed" to fit between [-1,1] (and also vertically "stretched" to maintain the "area" under the function).



Example, suppose we want to integrate f(x) = x^2 on the interval [1,9]. First, shift f(x) left 5, and then compressH and stretchV by 4. This gives the function g(x) = 4(4x+5)^2 on [-1,1] which will have the same area as f(x) on [1,9]. This function is then evaluated at the nodes.



These function values are then multiplied by appropriate weights. The weights are calculated using each node value and another Legendre polynomial (I can explain later if you'd like). Finally, the values are summed up giving the value of the integral.



The nice thing about the Gauss-Legendre Quadrature is its accuracy: a GL-N (using an Nth degree Legendre polynomial) can integrate a 2N-1 degree polynomial exactly. So, with N=15, you have 15 nodes to evaluate your function (multiplied by 15 weights), and this will integrate a 29 degree polynomial exactly. So any function that can be approximated by a 29th degree polynomial can be integrated using only 15 function calls. Nice.



Many calculators today actually use a Gauss-Kronrod method for integrating, which is an adaptive quadrature method that gives the ability for error estimation (the Gauss-Legendre method doesn't provide an error estimate). Kronrod inserts additional nodes (also using the original nodes) and new weights to get a better approximattion, thus giving an estimate of the error.



Now, if we choose to use N=3, the 3rd degree Legendre polynomial has exact zeros of 0, and +/-sqrt(.6), with weights of 8/9, 5/9, 5/9. Using the function g(x) above with these values gives an integral of 242.666666667.



If you look at Valentin's HP11C program, you can see the nodes of 0 and sqrt(.6), and also the weights of 8/9 and 5/9, however the program uses 8 and 5 and then divides by 9 at the end (actually 18 at the end, but this has to do with the transformations). All of the other pieces of the program are the transformations on f(x) and keeping a running total of the intermediate values. Valentin's program will be able to integrate a 5th degree polynomial exactly.



Last night I figured out how to get the HP50G to do the GL quadrature for any N (within reason). Quite a litte machine. For more information do a Google search with Gauss Legendre Mathworld, or Gauss-Kronrod. NEAT STUFF!!!!!!!!!



By the way, the TI83 uses the Gauss Kronrod method with I think up to G31-K63 points, meaning the TI83 should be able to integrate a 125 (2*63-1) degree polynomial exactly. But it might be only a G15-K31 for up to a 61 degree polynomial. I'll have to do a little more sleuthing.



Cheers, CHUCK


#31

Also, I just ran across reference to a 1980 article in the HP Journal that mentioned the HP34 using the Romberg method for integration. Not sure if this is what current HP's use or not.
CHUCK

#32

Hi, Chuck:

Chuck posted:

"Last night I figured out how to get the HP50G to do the GL quadrature for any N (within reason). Quite a litte machine."

    Congratulations but it's actually pretty easy. I can do that in the bare-bones HP-71B (i.e., no Math ROM needed) in less than 15 lines without needing those fancy "Legendre" and "Zeros" commands at all.

Best regards from V.


#33

Hi Valentin. I'm intrigued by the HP71B and your routine. I've been eyeing the 71B for a few years but just haven't gotten around to getting one. Maybe a Christmas present to myself. Anyway, are you able to actually calculate the Gauss-Legendre nodes and weights for any N on the fly in the routine? If so, I'm highly impressed. As far as I know the actual zeros for the Legendre polynomials and the weights have to be precalculated and stored in ROM (at least for N>5) in a table. Do you have the 71B routine available anywhere to look at. I haven't looked at this stuff for quite some time, so I might be still thinking in the dark ages. :)

CHUCK


#34

Hi again, Chuck:

Chuck wrote:

"I've been eyeing the 71B for a few years but just haven't
gotten around to getting one. Maybe a Christmas present to myself."

    I concur. It would be a most excellent present, as long as you can get one in pristine condition (easy), with the HP-IL ROM (easy), at a good price (easy), and with a Math ROM (not so easy). In the meantime you can get Emu71, an excellent, free emulator for it, and try some of the many programs for it featured in my Datafile articles, many of which you can get online for free in PDF format from my web page (google for 'albillo calculator'). They're quite varied. Another very good source would be to search for my challenges and long posts about it in this MoHP Forum.

"Anyway, are you able to actually calculate the
Gauss-Legendre nodes and weights for any N on the fly in the routine?"

    Why, of course ! Else, it wouldn't be funny ... :-)


If so, I'm highly impressed. As far as I know
the actual zeros for the Legendre polynomials and the weights have to be precalculated and stored in ROM (at least for
N>5) in a table. Do you have the 71B routine available anywhere to look at.

    It's intended to be published as part of a next article for Datafile, so you can have a look at it there in the very next issue, assuming you're a subscriber like many people in this Forum. Failing that, I'll make the article available for free online in my web page, within four months or so, as always. But I insist it's actually quite easy, just 13 or 14 lines of pretty simple BASIC code (though, as a bonus, it actually accepts a maximum error estimate for the requested integral and makes sure the computed value meets the specified accuracy). I've seen much more complicated programming tasks in my days, believe me.

Thanks for your interest and

Best regards from V.

#35

Valentin, you tease, you have talked about these economical routines that compute tha Gaussian weights and sampled abcissas "on the fly" in less than two dozen lines, and I am really intrigued about how you do it. Your 46 line routine in the thread from June is very compact and elegant, but as I understand it it is basically a 3-point Gaussian quadrature applied across N subdivisions of the domain of integration, so the function value is computed 3N times. But it sounds like you have something on hand that computes N-point Gaussian integration "on the fly" without prestoring the weights and abscissas in registers or variables. That fascinates me and I would love to see the details.

Chuck, I have a 49G+ but have no idea how to use its Legendre function, which I presume is similar to that of the 50G. (My understanding is that the 49G+ and 50G are basically the same calculator, with the newer calc have a better keyboard and more pleasing colour scheme, plus, I think, differences in speed and RAM.) I will try to educate myself and reproduce your approach.

Les

#36

Thanks for educating me about the TI83's routine.

I have long been interested in the error function, and it so happens that erfc(x) is equal to the function 1/Sqrt[-Pi*Log[t]] integrated from 0 to Exp[-x^2]. (I think, actually, that the 15C Advanced manual discusses this very transformation of the more familiar erfc(x) integral, the function 2/Sqrt[Pi]*Exp[-t^2] integrated from x to Infinity.)

Anyway, as x gets smaller, the bigger the area of integration gets and the longer quadrature on the integral takes. The singularity at the left end of the domain of integration makes life colourful too.

Using the Advantage Pac INTEG or the HP15C to compute erfc(2) using this integral really inspires a bit of a workout--my calcs seem to chug along a goodly couple of minutes to give 5 or 6 digits. The HP33S is a bit faster, and the 48G and 49G+ are faster still, but asking for more than 5 or 6 digits really drags things out to several minutes at least. Likewise with the 42S, though if you only want 4 or 5 digits it seems fast indeed.

The TI83's fnInt function is a whole different ball of wax. By adjusting the tolerance parameter, I can get anywhere from 4 digits almost instantly to 13 digits within 1 ULP (confirmed by digit shifting to see the extra digits) in less than a minute.

A few months ago I denigrated the TI83 as an ostensibly cheaply made piece of junk that looked, felt, and acted like the adolescent market it is pitched at. Well, this little area of strength really heightened my respect for the gizmo. I am even getting used to the algebraic operating system and the on-screen editing, so that when I revert to an HP I get confused!

I believe the TI83 fared poorly in Valentin's integration challenge of the summer, but then again so did a lot of the HPs. It had to do with how the problem was framed in many cases (e.g. the oscillatory nature of sin(x)*sin(x^2) dooms typical quadrature to ignominious failure, since the sampled absicissas miss so much information due to the oscillations and the interpolative nature of of classical quadrature approaches is mismatched).

I am familiar with the basics of Gauss-Legendre quadrature, but right now I am dying to know how Valentin computes the weights and sample points for a given arbitrary N.

Les


Edited: 11 Dec 2006, 9:19 p.m.

#37

Hi Chuck,

I have confirmed indeedy that my 49G+ has the Legendre and Zeros functions so if you can walk we through how to compute the weights (the actual node determination is easy) I would be much obliged.

Les


#38

Actually, cancel that--the formula is in Abramowitz and Stegun.

Basically, the weights are roots of the first deriviative of the Legendre polynomial in question, if I understand the formula correctly.

Les


#39

Hi Les. There are several different formulas for the weights. The form I used was


              2
2( 1-xi )
wi = ----------------
2 2
(n+1) (P(xi))


(I hope that looks correct.) Where xi are the node values and the P is the legendre polynomial of degree n+1. So here are my steps..


1)  7 LEGENDRE     (for n=7 polynomial)
2) 'x' Zeros (creates a list of nodes)
3) obj-> ->array (changes the list to a vector)
4) duplicate the vector
5) Next, I have a tiny function (program) with the above
weight formula called wi, which takes 7 as the argument
and spits out a large formula for xi
6) Use the function MAP which will take the node vector and
"map" the function onto it. Then ->NUM to get a nice
looking weight vector.
7) Next, SWAP the stack to get the nodes back to 1: and the
weights on 2:.
8) Next, I have another little function that transforms f(x) to
the origin. It looks like this...

<< -> F A B << F 'X=(A-B)/2*X+(A+B)/2' EVAL SUBST '(B-A)/2' EVAL * >> >>

(I realize this could be more efficient, but hey.)
This takes as input f(x) and the limits A and B and returns g(x).
9) Apply MAP which maps the nodes to the function. ->NUM make nice.
10) Finally!!!, DOT product the two vectors, and voila, the integral value.



All of these steps could be put into a single fairly small program, but I kept them separate for now. Tomorrow I have a three-hour road trip, so maybe I'll try and get it all down to the infamous "12-line program." :) :)



Hope this makes sense. If not, I'll be back tomorrow.

By the way, wi looks like:

<< -> N '2*(1-x^2)/((n+1)^2*LEGENDRE(N+1)^2)' >>


CHUCK

Edited: 12 Dec 2006, 1:33 a.m.


#40

I prefer your weight function over the one in A&S since yours, which of course is equivalent due to the relationship between the n+1-th Legendre polynomial and the derivative of the n-th one, doesn't require the derivative to be computed.

The only problem I have with the 49G+ so far is that ZEROS is darn slow once n gets above 4 or 5. I understand that the 50G is much faster--one of the reported advantages, I understand.

Right now I am working on extending Valentin's repeated 3 pt quadrature routine to a higher order. Just embedding the constants in the program is no hardship, since on 41, 42S, and 33S, for example, numeric constants in a program take just one line, no one line for each digit like in the 11C, 15C, and 67/97.

Will keep you posted.

Les

#41

The HP50G has clearly been improved in this regard.

On the HP49G+, ZEROS only seems to locate roots for a Legendre polynomial of degree less than 6. Above that, I either get zero returned, or an error "unable to factor".

Unless I am doing something wrong....

Les


#42

Quote:
The HP50G has clearly been improved in this regard.

The 49g+ and the 50g share exactly the same firmware. Both can be upgraded to version 2.09. There shouldn't be any difference in behavior between these two.

Marcus


#43

I will do the upgrade.

I have a demo of Hrast's HP41X for 49G+ that doesn't work with the newer firmware so when I attempted the upgrade a few months ago and experienced the incompatibility I reverted back to the older version.

Hrast has never been happy with the performance of his emulator for the 49G+ and frankly I don't seem to use it. I will take it off and upgrade my firmware.

Les


#44

The upgrade is done, and the VERSION command displays the right information.

Alas, unlike for Chuck and his 50g, the ZEROS command on my 49G+ will not give all of the roots of a 7th degree Legendre polynomial. The display flashes briefly "unable to factor" and only the 0 root is returned.

Chuck, could you double check this on your 50G? It really has me scratching my head.

Les


#45

Hi Les. I'm back from Olympia. The version on my 50G is 2.09. I'll run a little timed-test here with N=31....



It took about 4 seconds to make the polynomial.
Now the fun part...the zeros.....about 29 seconds. I have the modes set to RPN, Approx, Rigorous, Sim Non-Rational.


Now the bad news, the zeros closest to zero are accurate to 8 decimal places, while the zeros closest to +/-1 are only good to 3 decimal places. Hmmm, not good. (I checked these with Mathematica.) Thre decimal place accuracy on the zeros and weights will significantly alter the value. I'll look into this later.


CHUCK


#46

Jeez, I figured out my problem--I was trying to find ZEROS in exact mode (i.e., approx mode off). It seems that in symbolic mode ZEROS will only find roots if it can rationally factor the polynomial into a product of quadratics or linear factors. Legendre(7) factors into X and a 6th degree polynomial that is basically a cubic in x^2. You would think the 49G+ would symbolically solve the cubic in radicals. Doesn't. But, I set the approx mode and it works like a charm. How Rigourous and Simp Non Rational are set seems irrelevant--in APPROX mode, ZEROS finds the roots in decimal form, and for smaller polynomials seems quicker about it.

It would be interesting to see how slight imprecisions in the values of the Legendre zeros propagate in computing the quadrature.

Thanks for helping redirect me. I must admit, I don't know the 49G+ very well and the CAS features (i.e. symbolic vs. approx mode, etc.) sort of still bewilder me.

Les


#47

Interesting. Graph the N=31 polynomial, and zoombox around the far left roots. All heck breaks out due to the step size, creating a lot of false zeros. Changing the step to about .01 helps, but changing to .0001 causes more problems. It looks like the size of the polynomial is just beyond the accuracy of a handheld device. My hunch is that the numerical solving capabilities would parallel the graphing ability. It's just too much for the calculator. :(


#48

Hi Chuck,

I have had internet problems and haven't been able to post for a few days.

I did a little experimentation with the 49G+ and would appear that the keystroke sequence N LEGENDRE 'X' ZEROS begins loosing digits significantly once one goes beyond N = 10, at least when I compare the results with the tabulated values in Abramowitz and Stegun. (The HoMF computes Gaussian abscissas and weights up to N = 96, I think!) Even though the 49G+ gives 12 full digits for the abscissae for N=10, with rounding error computing the weights leads to a lost digit here and there in those values. For N=12 and above, it falls all apart, especially for the roots near unity.

In an infinite precsion enviroment, Gauss-Legendre quadrature looks insanely accurate. This tells me that threats to actual accuracy are the nature of function (i.e., is it reasonable to assume that a degree 2N-1 polynomial is a good approximation on the interval?), and the precision of the computed parameters. The latter issue is the Achilles heel of your proposed high-order routine, I think.

I wonder if Valentin's original routine is the best compromise after all. If rounding error is a huge issue in computing the higher order parameters, it makes more sense to me to do a 3-point quadrature on 10 subintervals rather than a 30-point quadrature on the whole interval. Maybe on the 49G+/50G there is a middle road here--we divide up the range of integration into subintervls but do say a 5 or 6 point quadrature on each. Would that give an improvement over Valentin's approach?

Still you can't beat Valentin's approach for sheer brevity and elegance. I am eager to see his work on computing quadrature weights and abscissae "on the fly". There is a routine in Numerical Recipes (gauleg), but it is a little longer than 11 lines!

best,

Les

#49

Hi Chuck, hope your road trip went or is going well.

When you can, can you confirm the ROM version on your 50g for me? I took the advice of others in this thread and updated my 49G+ ROM to 2.09, and still I can't get ZEROS to yield roots to Legendre polynomials of degree greater than 5. If you can do up to degree 7, like your suggested algorithm implies, I need to know what I am doing wrong!

Best,

Les

#50

Dig out your 1970's era IBM Scientific Subroutie Package manual for the IBM1130. There are several Fortran examples of Gaussian quadrature routines. It is my favorite integration method.


Possibly Related Threads…
Thread Author Replies Views Last Post
  Integration question and "RPN" mode comment Craig Thomas 16 5,842 12-05-2013, 02:32 AM
Last Post: Nick_S
  WP34s integration trapped in infinite loop Bernd Grubert 25 6,972 10-17-2013, 08:50 AM
Last Post: Dieter
  Prime cant find Memory available Dougggg 5 2,031 10-07-2013, 07:24 PM
Last Post: Han
  HP Prime integration Richard Berler 1 1,205 10-06-2013, 10:52 PM
Last Post: Helge Gabert
  HP Prime can't find how to use partfrac function ! dg1969 3 1,636 10-04-2013, 09:25 PM
Last Post: Joe Horn
  OT - A lucky find - Busicom Handy-LE LE-120A Cristian Arezzini 2 1,838 09-26-2013, 04:43 AM
Last Post: Cristian Arezzini
  integration on 39gII emulator Wes Loewer 29 7,089 06-07-2013, 05:58 PM
Last Post: Chris Smith
  WP-34S Integration Richard Berler 15 3,755 03-08-2013, 02:29 AM
Last Post: Walter B
  HP 34S integration Richard Berler 16 4,183 02-18-2013, 04:42 PM
Last Post: Marcus von Cube, Germany
  Craigslist find and fix - HP 19C Matt Fegenbush 12 3,749 08-13-2012, 02:57 PM
Last Post: Matt Fegenbush

Forum Jump: