Re: It's All About Romberg Integration Message #23 Posted by Kiyoshi Akima on 23 Oct 2008, 2:27 p.m., in response to message #22 by Les Wright
The Romberg algorithm breaks down in this case in the spurious agreement of early results. Taking the "straight-up textbook way":
F0 = FRAC(0.0) = 0.0
F1 = FRAC(1.6) = 0.6
F2 = FRAC(3.2) = 0.2
F3 = FRAC(4.8) = 0.8
F4 = FRAC(6.4) = 0.4
T1 = (F0+F4)*(6.4/2) = 1.28
T2 = (F0+2F2+F4)*(6.4/4) = 1.28
Since these two values agree, the algorithm "decides" there's no need to continue. Had it continued, it would have found that:
T3 = (F0+2(F1+F2+F3)+F4)*(6.4/8) = 2.88
which doesn't agree with the previous value and thus gone on, using more and more function evaluations, until the estimates (aided by the Richardson Extrapolation) converge. This is the behavior seen when the upper limit is changed to 6.3, 6.39, 6.399, etc.
In the case of the built-in schemes of all HPs (with built-in schemes, of course), RPL and RPN alike, a nonlinear transformation is performed in an usually successful effort to avoid the effects of periodic functions. In this particular case it fails, as the third estimate also agrees with the first two, "convincing" the scheme that it has found the correct answer.
I didn't start out to find this. As I stated earlier, I was simply comparing various schemes on an integral that should have evaluated to approximately pi. It was only when the built-in scheme came out that wildly off that I dug in deeper.
Of course, the "correct" way to numerically evaluate this integral, assuming one had to do it numerically, is as:
6*INTG(0,1,FRAC(x),x) + INTG(0,0.4,FRAC(x),x)
exploiting the periodicity and at the same time removing the singularities.
|