--- Al Chou <[EMAIL PROTECTED]> wrote:
> --- Phil Steitz <[EMAIL PROTECTED]> wrote:
> [deletia]
> > OK, long-winded disclaimer aside, here is how I see the task list ordered:
[deletia]
> > * Framework and implementation strategie(s) for finding roots or
> real-valued
> > functions of one (real) variable. Here again -- largely done. I would
> > prefer
> > to wait until J gets back and let him submit his framework and R. Brent's
> > algorithm. Then "our" Brent's implementation and usage can be integrated
> > (actually not much to do, from the looks of the current code) and I will
> add
> > my "bean equations" stuff (in progress).
>
> I may have time to submit my Ridders' method implementation using J.'s
> framework before he returns 2 days hence. Should I bother to try, or should
> I
> wait until he submits his code as a patch via Bugzilla?
Well, I've just spent some time over the past 3 days reminding myself pf some
of the things that are so hard about numerics.
I was testing my Ridders' method implementation and couldn't understand why it
took so many iterations to converge and still not be within the requested
accuracy of the known root I asked it to find. I used a simple quintic
(x+1)(x+0.5)(x)(x-0.5)(x-1) as the function whose roots I want to find, and I
made sure to give upper and lower bounds that I know bracket one and only one
root. When trying to find the roots at x = +- 0.5 my solver had no trouble
(though I didn't ask it how many of the 100 iterations it was allowed that it
actually used, until later), but the root at x = 0 was never within even a
factor of 15 of the requested 1e-6 accuracy even when allowed to take up to 200
iterations (actually, I used this test case first, which was what prompted me
to try the larger-valued roots in case I was seeing some loss of precision or
roundoff error effect).
BTW, in the process of using Herr Pietschmann's root finder framework, I
discovered a bug in setMaximalIterationCount (it sets
defaultMaximalIterationCount instead of maximalIterationCount).
I then decided to try Brent W.'s bisection solver, which converged to the
desired root to within its requested accuracy (1e-9) in 26 or 27 iterations
even for the root at x = 0. At this point I asked my Ridders' method how many
iterations it took to find x = 0.5, and it said 1, and I realized that was
probably because my bracket values were symmetric (or close enough) about the
root, so its midpoint evaluation of the function found the root by coincidence.
When I made sure the bracket values weren't symmetric about that root, I was
back to 146 iterations or more and not getting to within the requested accuracy
of the root location.
So I pulled out Herr Pietschmann's Brent method class and tested it, and it
threw an exception telling me, "Possibly multiple zeros in interval or ill
conditioned function."
The morals of the story are:
- More-sophisticated algorithms that are supposed to converge faster don't
always do so
- It's easy to outsmart yourself and create code that's too finicky for
non-numericist users.
As someone said recently on the list, a typical user probably is more
interested in an algorithm that's guaranteed to converge to a root (if there is
one) than in the rate of convergence, as long as it's not too ridiculously
slow. Given that we've repeatedly determined that commons-math is not to be a
general numerical mathematics library, I think now that we should provide only
a bisection method in the initial release (assuming we achieve one) and spend
time later making our implementations of the more sophisticated algorithms more
user-friendly, if we find they're even needed. I believe we've let ourselves
go down the path of as-yet-unjustified optimization in our designs, because we
know of algorithms that are supposed to be "better". I also have a greater,
first-hand, appreciation of the subtleties in NR's code to make it more robust
for the user, and I believe we can only achieve that level of robustness if we
take enough time -- which we should not prior to the initial release, because
that will be too much time.
Finally, having used the Pietschmann root finder framework, I think it needs
some modification to make it more user-friendly. As a lay user, I would have
been much happier dealing with Brent W.'s interface than Herr Pietschmann's,
which was kind of cumbersome. I think, though, with a little slimming down, it
would be quite workable.
Al
=====
Albert Davidson Chou
Get answers to Mac questions at http://www.Mac-Mgrs.org/ .
__________________________________
Do you Yahoo!?
Yahoo! Calendar - Free online calendar with sync to Outlook(TM).
http://calendar.yahoo.com
---------------------------------------------------------------------
To unsubscribe, e-mail: [EMAIL PROTECTED]
For additional commands, e-mail: [EMAIL PROTECTED]
