[SciPy-User] Is this weird answer a bug in optimize.newton?

[Thomas Baruchel]

Hi,

with my students, this afternoon, I encountered the following annoying
behaviour of Scipy:

>>> from scipy import optimize
>>> optimize.newton(lambda x: x**3-x+1, 0)
0.999999989980082

This is very weird, because the equation is rather simple with simple coefficients,
and the initial point is simple also. Any teacher or student could fall into this
case when trying to use optimize.newton.

Of course, the answer is absolutely wrong. First, some facts:

   * using the function with the derivative gives the correct answer (we are not
     in some tricky case);
   * my students wrote pieces of code with Newton-Raphson algorithm returning
     the expected answer even without the derivative; they guessed the value of
     the derivative either by some (f(x+eps)-f(x-eps))/(2*eps) or
     (f(x+eps)-f(x))/eps, etc. All these methods worked well.

I had a glance in the source code and found that:

   * the function is actually thinking the answer is correct (the answer has
     nothing to do with iterating too many times) because two consecutive values
     of x are close "enough".
   * the answer is returned after only 3 iterations (correct answer requires
     about 20 iterations).

I think it has to do with the way the derivative is computed by Scipy (in order to
call the f function only once in each iteration, the secant is computed by
using f(x) for consecutive values of x (rather than by using some epsilon);
furthermore, in order to start this method the function has to provide an
arbitrary "second" value of x).

Is it a bug? Regards.

-- 
Thomas Baruchel