"J.Pietschmann" <[EMAIL PROTECTED]> wrote in message news:[EMAIL PROTECTED] > Phil Steitz wrote: > > Can you provide a math reference desribing this? > > H.M.Antia: Num. Methods for Scientists and Engineers. > > > I have been referring > > to Atkinson and I am having a hard time following the implementation. > > What exactly do you mean by "only inverse quadratic interpolation"? > > Brent's method requires a bracketed root as a start. The primary > method for shrinking the bracket is inverse quadratic interpolation. > This means you get three points > x0,y0 x1,y1 x2,y2 x0<x2<x1 > and interpolate a parabola > x=a*y^2+b*y+c > and because your goal is to find the x for y=0, your next estimate > for the root is c (set y=0). The convergence is of the order 1.8, > which is better than the secant method (~1.4 if non bracketing, 1..1.2 > on average if bracketing). > The full Brent algorithm measures how well the interval shrinks, and > resorts to bisection if progress is too slow. I didn't implement this > because I forgot to take the book with me and had only a hazy memory > of the control parameters. All in all the algorithm combines > near-guaranteed convergence (occasionally problems might falsely > detected as ill-conditioned) with near Newton-speed. > > > What do you have in mind re: plausibility? > If the interval to search is (x0,x1), then absAccuracy should be of > the order of max(abs(x0),abs(x1))*relAccuracy. Achievable relative > accuracy depends on underlying hardware, although nowadays basically > everyone uses IEEE 754 (means, uh, 52 bit for double? Damn, have to > look it up!). Both relative and absolute accuracy of function > evaluation are also important, which is the reason to let the user > override defaults. > This means if relAcc is given then reject absAcc if > max(abs(x0),abs(x1))*relAcc+c*absAcc == max(abs(x0),abs(x1))*relAcc > for some predermined constant c in 0.2..1. > > > I guess I like your rootfinding framework better than Brent's (Here I > > mean Brent Worden, the famous commons-math contributor, not the > > numerical analyst). I suggest that we agree to use your interfaces and > > UnivariateRealSolverImpl base,include Brent's bisection implementation > > strategy and modify his CDF inversion to use the new rootfinding framework. > > No argument against :-) I took special care for the JavaDocs, which > seems to pay off...
I agree. The this looks like a very solid framework. One suggestion I would like to make, is instead of a both a firstDirevative and secondDerivate method to evaluate the derivates. Create a single getDerivate() method that returns a UnivariateRealFunction. That way if a user needs the n-th derivate, they just call the getDerivate() method n times, once on the original function and once on each of the returned functions. That way, common-math supports creating whatever degree derivate a method might need without ever having to change the framework. We provide the maximum amout of derivate creation support to the user while providing us with the minimual amount of maintenance. > > > I do have a few small questions/comments on the framework: > > > > 1. Having solve() *require* brackets makes it impossible (or at least, > > un-natural) to add Newton's method or any other method that just starts > > with one initial value. > > Why? Start with -4000000,+4000000000 or so. The algorithm is not > *required* to start with a bracket, just with an interval. Individual > solver implementations should document whether they require a bracket. > Apart from this, there is always the possiblity to add another method. Why not allow the user to supply any number of initial values and design the solvers to compensate as best they can when not enough values are provided. Each algorithm has criteria about the initial values that must be met before root finding can be attempted. If the user provided initial values do not satisfy the criteria, the solver should first attempt to morph the initial values into a set of values that do satisfy the criteria. If this can not be accomplish, then the solver would raise an exception. This would take away some of the user's responsibility of knowing how these algorithms actually work and place it on us to create more robust components. > > > 2. I am curious as to why the "result" property is there. How exactly > > do we expect clients to make use of this? > > The client can ask for the last result as long as no further attempt > to solve for another root was made. I found this handy. I don't insist > on keeping it. > > > 3. What were you thinking about putting into the base solve() > > implementation as diagnostics? Do you mean things like checking to make > > sure the bracket is good?\ > > No, just adding a message indicating that an unimplemented method > was called. Some frameworks don't display the type of the exception > to the user and rely on the message. > > > 4. Thinking of the developers who will use this stuff, I feel like we > > need a way to insulate them from having to think about rootfinding > > strategy choice. As Al has pointed out, we are not (at least yet ;-)) > > in the AI business, so we can't automagically make the best choice for > > them; but it might be a good idea to somehow specify a default strategy. A simple and flexible way to make the "best" choice for the user, is to control the creation of solver instances. This could be done via a factory or factory method. The creation could be based on the class of function the user want to evaluate (real, polynomial, rational, etc.) Also, we could create a solver chain that provides the ability to link various solvers together so if a solver fails to find a root, the next solver in the chain is given a chance to solve it. Either, eventually one of the solvers finds a root or all solvers fail and an exception is raised. The complexity of creating the "best chain" would be contained in the factory or factory method, totally hidden from the user. This again would relieve users of having some numerical method knowledge about root finding methods. > > Unfortunately, the "safe" choice would likely be bisection. > > Brent's algorithm is quite safe in this respect. I'm by no means beholden to bisection. I just used it because I needed something to use for the distribution work. Brent's algorithm is guarranteed to converge to a root, provided one is bracketed by the initial values. So, if a bracket can be provided or generated, it is just as good as bisection as far as safety and it should almost always beat it in speed of convergence. > I'll see whether I can complete the implementation near term. Unfortunately > I'll go on vacation on saturday and will be offline until 2003-06-09. > > > J.Pietschmann Brent Worden http://www.brent.worden.org --------------------------------------------------------------------- To unsubscribe, e-mail: [EMAIL PROTECTED] For additional commands, e-mail: [EMAIL PROTECTED]
