Trying to find Integration Challenge - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: Trying to find Integration Challenge (/thread-103994.html) |
Trying to find Integration Challenge - Les Wright - 12-10-2006 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
Re: Trying to find Integration Challenge - Gerson W. Barbosa - 12-10-2006 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.
Re: Trying to find Integration Challenge - Les Wright - 12-10-2006 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
Re: Trying to find Integration Challenge - Gerson W. Barbosa - 12-10-2006 Hi again Les, I don't want to spoil your searching fun but I think this is one of threads: Please let us know how your TI83 handles those! :-) Regards, Gerson. ------------
Edited: 10 Dec 2006, 3:19 p.m.
Gaussian Quadrature - Chuck - 12-11-2006 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). Re: Gaussian Quadrature - Chuck - 12-11-2006 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. Re: Gaussian Quadrature - Valentin Albillo - 12-11-2006 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."
Best regards from V.
Re: Gaussian Quadrature - Chuck - 12-11-2006 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. :) Re: Gaussian Quadrature - Les Wright - 12-11-2006 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
Re: Gaussian Quadrature - Valentin Albillo - 12-11-2006 Hi again, Chuck: Chuck wrote:
"I've been eyeing the 71B for a few years but just haven't
"Anyway, are you able to actually calculate the
Thanks for your interest and
Best regards from V.
Re: Gaussian Quadrature - Les Wright - 12-11-2006 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.
Re: Gaussian Quadrature - Les Wright - 12-11-2006 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
Re: Gaussian Quadrature - Les Wright - 12-11-2006 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
Re: Gaussian Quadrature - Chuck - 12-12-2006 Hi Les. There are several different formulas for the weights. The form I used was 2 (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) 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.
Re: Gaussian Quadrature - Les Wright - 12-12-2006 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
Re: Gaussian Quadrature - Les Wright - 12-12-2006 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
Re: Gaussian Quadrature - Marcus von Cube, Germany - 12-12-2006 Quote: 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
Re: Gaussian Quadrature - Les Wright - 12-12-2006 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
Re: Gaussian Quadrature - David Smith - 12-12-2006 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.
Re: Gaussian Quadrature - Les Wright - 12-12-2006 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
Re: Gaussian Quadrature - Les Wright - 12-12-2006 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
Re: Gaussian Quadrature - Chuck - 12-12-2006 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.... Re: Gaussian Quadrature - Les Wright - 12-12-2006 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
Re: Gaussian Quadrature - Chuck - 12-12-2006 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. :(
Re: Gaussian Quadrature - Les Wright - 12-16-2006 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
|