[SciPy-Dev] Comments on optimize.newton function

[josef.pktd at gmail.com]

On Sun, May 22, 2011 at 1:55 AM, 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))

Wouldn't this be easier (for derivatives), and maybe be more stable, taking logs

np.log(rh) -  kelvin/x + np.log(..) ...  ?

(independently of any improvement to the solvers)

Josef