[SciPy-dev] genetic algorithm, number theory, filterdesign,zerofinding
oliphant at ee.byu.edu
Thu Apr 11 17:50:36 CDT 2002
> Another standard solution is to pass in two tolerances, one absolute and the other relative, say x_atol, x_rtol, and then do something like
> tol = x_atol + x_rtol*abs(current_estimate)
> x_rtol is usually bit more than the smallest number such that 1 + x_rtol != 1.
> This number is often computed during installation of large packages, such as scipy, and made available somehow to all the routines. With python it would be easy to assign it as a default. Or we could just check
> x == x + dx # oh oh (I don't really like this)
> > > > > 4. newton should probably check for blowup, as
> > this is not uncommon
> > > >
> > > > A blowup of what? The second derivative becoming very large?
> > >
> > > newton tends to head off to infinity if the function is
> > bumpy and the initial estimate is insufficiently close.
> This was a whole subject way back when. Weingarten,Dekker developed a routine that did all sorts of checks on convergence and fell back on bisection when things weren't going well. Brent added second order extrapolation, giving the Weingarten,Dekker,Brent zero finder. It's a total b*tch to understand, and even worse in the original Algol60. There is actually a fairly nice Fortran implementation --- with lots of gotos --- due to Moler floating around on the net, its called zero or some such, and I don't recall if it is the original, or Brent's improvement. This might be the way to go. The gotos actually make it easier to understand, as it is best seen as a finite state machine, and those can look pretty ugly in structured code.
Great, it would be easy to throw this in with an f2py-generated wrapper.
> The problem of doing this in python is that all the checks and whatnot add *lots* of overhead. The second order extrapolation does improve things though. I also prefer my own method of doing the extrapolation, but if its packaged up and works, who cares.
People are working on Python-compilers. At some-point we would like to be
able to compile the solvers using that. I wouldn't expect this in the
next year, but I believe it will happen to some degree.
I don't mind wrapping code, you can see that I've done a lot of it
throughout SciPy. Most of SciPy is just interfaces to other-people's
More information about the Scipy-dev