[SciPy-Dev] Comments on optimize.newton function

[Charles R Harris]

On Sat, May 21, 2011 at 11:55 PM, Gökhan Sever <gokhansever at gmail.com>wrote:

> On Sat, May 21, 2011 at 7:50 PM, Charles R Harris
> <charlesr.harris at gmail.com> wrote:
> > Yeah, that left me thinking that we could really use a bracketing
> function.
> > The brent optimizer must have something like that it uses to bracket
> > minimums and perhaps that could be adapted. I know there are bracketing
> > methods out there, IIRC there is one in NR.
>
> I have managed to get those bracketed solvers working after plotting
> my function. brenth seems converging the fastest but still providing
> the search interval is a problem for me. I am probably skipping these
> solvers.
> >
> > IIRC, the old version blew up when the root was at zero, the problem was
> > posted on the list.
> >
> > ? I believe there is only one newton method in scipy, but we moved it
> into
> > the zeros module and deprecated the version at the old location. It has
> > since been removed.
>
> Yes, you are right. At the current source repository, there is only
> one newton in scipy.optimize which resides in zeros.py
>
> >
> > I would just stick to the secant method in the general case unless the
> > derivative is easy to come by. Note that newton is one of the 'original'
> > functions in scipy, the other 1d zero finders came later. So if you can
> make
> > an improved cythonized version I don't see any reason not to use it. If
> you
> > do make cythonized versions it might be worth implementing the Newton and
> > secant parts separately and make the current newton a driver function.
>
> I have figured out the derivative option and tested fsolve and newton
> with fprime arg provided. Still slower comparing to the secant method.
> Most likely, that the derivative function requires a bit calculation
> to be evaluated. As you can see, the function is:
>
> cpdef double petters_solve_for_rw(double x, double rd, double rh):
>    return rh - exp(kelvin/x) * (x**3 - rd**3) / (x**3 - rd**3 * (1.0 -
> kappa))
>
>
You could also try rewriting this in various ways. For instance

cpdef double petters_solve_for_rw(double x, double rd, double rh):
   return rh*( x**3 - rd**3 * (1.0 - kappa)) - exp(kelvin/x) * (x**3 -
rd**3)

or

cpdef double petters_solve_for_rw(double x, double rd, double rh):
   return rh*(1  - (1 - kappa)*y**3) - exp(y*kelvin/rd) * (1 - y**3)

or

cpdef double petters_solve_for_rw(double x, double rd, double rh):
   return rh*kappa*y**3 - (exp(y*kelvin/rd) - rh) * (1 - y**3)



where x = rd/y. The last might allow you to bracket things fairly easily,
i.e., (exp(y*kelvin/rd) - rh) * (1 - y**3) has to be >0 if you expect y>0

<snip>

Chuck
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20110522/abc7c280/attachment.html>