Romberg implementation
in [HP48]
|
Prev: Programming
Next: Uploading programs to 50g
From: Wes on 2 Dec 2007 12:00 John Meyers, I've been playing around with integration techniques the past few days and came across something you wrote about the Romberg method a few months back. You wrote: vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv The Romberg method of numeric integration is a "sampling" method, which makes successive passes through the interval of integration, each time evaluating the integrand at points which are between the points previously used (or closer to the endpoint at the two extremes). The points are selected by applying a third-degree polynomial to what would otherwise be interval halving, causing the "sampling" points to be more "bunched up" near the endpoints and most sparse in the middle of the interval, ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ What caught my eye was the part about the "third-degree polynomial". Can you be more specific? Is there a specific cubic that gives you the intervals? And how do you get the intervals from the cubic? the roots? It appears that HP's implementations may indeed be halving the interval. John Edry's documentation for his QAG48 program mentions that the 48's built-in integrator fails (as does the 50g's) for functions like integral(0,pi,cos(512*x),x) or any multiple of 512. This would certainly lead me to think that it is bisecting the interval. (Maxima's Romberg fails for multiples of 8.) This weekend, I tried my hand at writing a Gauss-Kronrod program that recursively bisects the interval if the error for that interval is too large. After reading what you wrote, I'm wondering if it would better to subdivide into uneven intervals. Then again, the Gauss-Kronrod method already samples points distributed more towards the outer edges of the interval, so perhaps this would be redundant. -wes
From: mjc on 3 Dec 2007 13:06 On Dec 2, 9:00 am, Wes <wjltemp...(a)yahoo.com> wrote: > John Meyers, > > I've been playing around with integration techniques the past few days > and came across something you wrote about the Romberg method a few > months back. You wrote: > > vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv > The Romberg method of numeric integration > is a "sampling" method, which makes > successive passes through the interval of integration, > each time evaluating the integrand at points > which are between the points previously used > (or closer to the endpoint at the two extremes). > > The points are selected by applying a third-degree polynomial > to what would otherwise be interval halving, causing the > "sampling" points to be more "bunched up" near the endpoints > and most sparse in the middle of the interval, > ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ > > What caught my eye was the part about the "third-degree polynomial". > Can you be more specific? Is there a specific cubic that gives you > the intervals? And how do you get the intervals from the cubic? the > roots? > > It appears that HP's implementations may indeed be halving the > interval. John Edry's documentation for his QAG48 program mentions > that the 48's built-in integrator fails (as does the 50g's) for > functions like integral(0,pi,cos(512*x),x) or any multiple of 512. > This would certainly lead me to think that it is bisecting the > interval. (Maxima's Romberg fails for multiples of 8.) > > This weekend, I tried my hand at writing a Gauss-Kronrod program that > recursively bisects the interval if the error for that interval is too > large. After reading what you wrote, I'm wondering if it would better > to subdivide into uneven intervals. Then again, the Gauss-Kronrod > method already samples points distributed more towards the outer edges > of the interval, so perhaps this would be redundant. > > -wes AFAIK, the standard Romberg method halves successive intervals, using extrapolation to estimate the integral. It assumes that the integral can be represented in the form a0 + a1*h^2 + a2*h^4 + ... + ak*h^(2k) + ... By evaluating trapazoidal sums for h, h/2, h/4, ..., we can successively eliminate a1, a2, ... to get better and better estimates for a1, the integral. Look up Romberg integration and Richardson extrapolation (the basis for it). In a sense, the opposite of Romberg integration is adaptive integration, which looks at individual intervals and tries to decide if the intervals need to be further subdivided.
From: John H Meyers on 3 Dec 2007 22:54 On Sun, 02 Dec 2007 11:00:33 -0600, Wes wrote: > I've been playing around with integration techniques the past few days > and came across something written about the Romberg method: >> The Romberg method of numeric integration >> is a "sampling" method, which makes >> successive passes through the interval of integration, >> each time evaluating the integrand at points >> which are between the points previously used >> (or closer to the endpoint at the two extremes). >> >> The points are selected by applying a third-degree polynomial >> to what would otherwise be interval halving, causing the >> "sampling" points to be more "bunched up" near the endpoints >> and most sparse in the middle of the interval, > What caught my eye was the part about the "third-degree polynomial". > Can you be more specific? Is there a specific cubic that gives you > the intervals? And how do you get the intervals from the cubic? There is only one 3rd degree polynomial mapping an interval's endpoints and midpoint to themselves, while having a zero derivative at the endpoints, which this information is sufficient to derive :) > [and] the roots? To map "halved intervals" to the actual sample points, you simply apply the polynomial to the "halfway" points; there are no "roots" to take. The function value at each sample point is then multiplied by a weight which is more or less the derivative of the polynomial at that point. The endpoints (where the derivative is zero) are never used -- and that's basically the whole idea, allowing the integration of functions which may "approach infinity" at the interval endpoint(s), although the integral itself may be finite. For example, 5 'N' STO 2 FIX '\.S(1,0,LN(X)^N,X)' \->NUM (Ans=119.35, IERR=1.19) Infinite as X->0, but exact integral is FACT(N) http://groups.google.com/group/comp.sys.hp48/msg/5c4fea79054a2bdf > It appears that HP's implementations may indeed > be halving the interval. Why not collect some data and see; on any HP48/49/50: { } 'L' STO \<< \-> X \<< X 'L' OVER STO+ \>> \>> 'F' STO '\.S(0,1,F(X),X)' \->NUM L \->Q\pi Surprised? An "HP Journal" article by William Kahan (reference below) somewhat explains the original HP34C implementation of both the numeric solver and integrator, along with a textbook reference for the Romberg formulas. All the same, there was a small, subtle bug, I think identical in HP34C and HP15C but later corrected, producing a bit more error than need have existed (often exceeding IERR for the above example, and not consistent for different number display settings that evaluated exactly the same number of function samples, at the same points, which proved the existence of some original bug). > John Edry's documentation See http://groups.google.com/group/comp.sys.hp48/msg/ad5bc70dddf634be http://groups.google.com/group/comp.sys.hp48/msg/8ee3e862d869b5ab "Working with INTEG." HP-15C Advanced Functions Handbook, Hewlett-Packard Co., 1982. Kahan, W.M., "Handheld Calculator Evaluates Integrals," Hewlett-Packard Journal, 31:8, August 1980. http://en.wikipedia.org/wiki/William_Kahan [r->] [OFF]
From: Wes on 6 Dec 2007 16:36 John, Thanks for the reply. It took me a little while to decipher what you were saying, but I think finally got it. For estimating integral(a,b,f(x),x)... > There is only one 3rd degree polynomial mapping > an interval's endpoints and midpoint to themselves, > while having a zero derivative at the endpoints, > which this information is sufficient to derive :) For anyone else who is interested, what he is saying is the cubic, p(x), will have the properties p(a)=a p(b)=b p'(a)=0 p'(b)=0 The stipulations that p(midpoint)=midpoint is not required to find the cubic, but will come about from the symmetry. (At first, I incorrectly thought you were saying: p(a)=f(a) p(b)=f(b) p(midpoint)=f(midpoint) p'(a)=0 p'(b)=0 which would require a 4th degree polynomial.) >> [and] the roots? > To map "halved intervals" to the actual sample points, > you simply apply the polynomial to the "halfway" points; > there are no "roots" to take. In other words, if you're working on a subinterval from q to r and you want to subdivide it into two pieces, instead of splitting it at (q+r)/ 2, split it at p((q+r)/2) where p(x) is the cubic above. > Why not collect some data and see; Best advice of all. > on any HP48/49/50: > > { } 'L' STO > \<< \-> X \<< X 'L' OVER STO+ \>> \>> 'F' STO > '\.S(0,1,F(X),X)' \->NUM L \->Q\pi > > Surprised? Indeed. Surprised at the numbers (which made everything make sense) and surprised that I hadn't thought of doing that myself as I've used something similar in the past to peek in on the solver. I used the same idea to see what Maxima's Romberg is doing and found that it's doing straight bisection. (I'd like to do the same for the TI-83+/89, but neither allows functions to have any side-effects what- so-ever. If anybody knows a way around this, please let me know.) > Kahan, W.M., "Handheld Calculator Evaluates Integrals," > Hewlett-Packard Journal, 31:8, August 1980. A good read. Thanks. And thanks to Don for posting the link to HP- Journal. -wes
From: John H Meyers on 7 Dec 2007 22:00
On Thu, 06 Dec 2007 15:36:25 -0600, Wes wrote: > I used the same idea to see what Maxima's Romberg is doing > and found that it's doing straight bisection. A little Google searching turned up, for me, only the same; whoever wrote the Wikipedia article even claims "Romberg's method evaluates the integrand at equally-spaced points." How well does it integrate 'LN(X)^N' (for fixed integer N) over [0,1] ? Many methods have a bit of difficulty calculating at one of those endpoints, but it seems to be no problem with Kahan's approach, which could be viewed as the integration of a transformed function at evenly spaced points, omitting the endpoints (although the limit of the transformed function could apparently also be infinite at those endpoints). Was that Kahan's own idea, or was it already in the reference he cited? (I once looked it up, but have long forgotten). [r->] [OFF] |