[Scipy-svn] r6501 - branches/0.8.x/scipy/optimize

scipy-svn@scip... scipy-svn@scip...
Tue Jun 15 21:31:21 CDT 2010


Author: charris
Date: 2010-06-15 21:31:20 -0500 (Tue, 15 Jun 2010)
New Revision: 6501

Modified:
   branches/0.8.x/scipy/optimize/minpack.py
   branches/0.8.x/scipy/optimize/zeros.py
Log:
DEP: Deprecate minpack.bisection and move minpack.newton to zeros.newton.
Also fix ticket #1198.


Modified: branches/0.8.x/scipy/optimize/minpack.py
===================================================================
--- branches/0.8.x/scipy/optimize/minpack.py	2010-06-15 01:33:56 UTC (rev 6500)
+++ branches/0.8.x/scipy/optimize/minpack.py	2010-06-16 02:31:20 UTC (rev 6501)
@@ -7,8 +7,7 @@
 
 error = _minpack.error
 
-__all__ = ['fsolve', 'leastsq', 'newton', 'fixed_point',
-           'bisection', 'curve_fit']
+__all__ = ['fsolve', 'leastsq', 'fixed_point', 'bisection', 'curve_fit']
 
 def check_func(thefunc, x0, args, numinputs, output_shape=None):
     res = atleast_1d(thefunc(*((x0[:numinputs],)+args)))
@@ -518,6 +517,9 @@
     interval has been found.
 
     """
+    msg  = "minpack.newton is moving to zeros.newton"
+    warnings.warn(msg, DeprecationWarning)
+
     if fprime is not None:
         # Newton-Rapheson method
         p0 = x0
@@ -620,10 +622,13 @@
     func(a) * func(b) < 0, find the root between a and b.
 
     """
+    msg = "minpack.bisection is deprecated, use zeros.bisect instead"
+    warnings.warn(msg, DeprecationWarning)
+
     i = 1
     eva = func(a,*args)
     evb = func(b,*args)
-    if eva*evb < 0:
+    if eva*evb >= 0:
         msg = "Must start with interval where func(a) * func(b) < 0"
         raise ValueError(msg)
     while i <= maxiter:

Modified: branches/0.8.x/scipy/optimize/zeros.py
===================================================================
--- branches/0.8.x/scipy/optimize/zeros.py	2010-06-15 01:33:56 UTC (rev 6500)
+++ branches/0.8.x/scipy/optimize/zeros.py	2010-06-16 02:31:20 UTC (rev 6501)
@@ -1,4 +1,3 @@
-
 import _zeros
 from numpy import finfo
 
@@ -7,13 +6,14 @@
 # not actually used at the moment
 _rtol = finfo(float).eps * 2
 
-__all__ = ['bisect','ridder','brentq','brenth']
+__all__ = ['newton', 'bisect', 'ridder', 'brentq', 'brenth']
 
 CONVERGED = 'converged'
 SIGNERR = 'sign error'
 CONVERR = 'convergence error'
 flag_map = {0 : CONVERGED, -1 : SIGNERR, -2 : CONVERR}
 
+
 class RootResults(object):
     def __init__(self, root, iterations, function_calls, flag):
         self.root = root
@@ -25,6 +25,7 @@
         except KeyError:
             self.flag = 'unknown error %d' % (flag,)
 
+
 def results_c(full_output, r):
     if full_output:
         x, funcalls, iterations, flag = r
@@ -36,6 +37,101 @@
     else:
         return r
 
+
+# Newton-Raphson method
+def newton(func, x0, fprime=None, args=(), tol=1.48e-8, maxiter=50):
+    """Find a zero using the Newton-Raphson or secant method.
+
+    Find a zero of the function `func` given a nearby starting point `x0`.
+    The Newton-Rapheson method is used if the derivative `fprime` of `func`
+    is provided, otherwise the secant method is used.
+
+    Parameters
+    ----------
+    func : function
+        The function whose zero is wanted. It must be a function of a
+        single variable of the form f(x,a,b,c...), where a,b,c... are extra
+        arguments that can be passed in the `args` parameter.
+    x0 : float
+        An initial estimate of the zero that should be somewhere near the
+        actual zero.
+    fprime : {None, function}, optional
+        The derivative of the function when available and convenient. If it
+        is None, then the secant method is used. The default is None.
+    args : tuple, optional
+        Extra arguments to be used in the function call.
+    tol : float, optional
+        The allowable error of the zero value.
+    maxiter : int, optional
+        Maximum number of iterations.
+
+    Returns
+    -------
+    zero : float
+        Estimated location where function is zero.
+
+    See Also
+    --------
+    brentq, brenth, ridder, bisect -- find zeroes in one dimension.
+    fsolve -- find zeroes in n dimensions.
+
+    Notes
+    -----
+    The convergence rate of the Newton-Rapheson method is quadratic while
+    that of the secant method is somewhat less. This means that if the
+    function is  well behaved the actual error in the estimated zero is
+    approximatly the square of the requested tolerance up to roundoff
+    error. However, the stopping criterion used here is the step size and
+    there is no quarantee that a zero has been found. Consequently the
+    result should be verified. Safer algorithms are brentq, brenth, ridder,
+    and bisect, but they all require that the root first be bracketed in an
+    interval where the function changes sign. The brentq algorithm is
+    recommended for general use in one dimemsional problems when such an
+    interval has been found.
+
+    """
+    if fprime is not None:
+        # Newton-Rapheson method
+        p0 = x0
+        for iter in range(maxiter):
+            myargs = (p0,) + args
+            fval = func(*myargs)
+            fder = fprime(*myargs)
+            if fder == 0:
+                msg = "derivative was zero."
+                warnings.warn(msg, RuntimeWarning)
+                return p0
+            p = p0 - func(*myargs)/fprime(*myargs)
+            if abs(p - p0) < tol:
+                return p
+            p0 = p
+    else:
+        # Secant method
+        p0 = x0
+        if x0 >= 0:
+            p1 = x0*(1 + 1e-4) + 1e-4
+        else:
+            p1 = x0*(1 + 1e-4) - 1e-4
+        q0 = func(*((p0,) + args))
+        q1 = func(*((p1,) + args))
+        for iter in range(maxiter):
+            if q1 == q0:
+                if p1 != p0:
+                    msg = "Tolerance of %s reached" % (p1 - p0)
+                    warnings.warn(msg, RuntimeWarning)
+                return (p1 + p0)/2.0
+            else:
+                p = p1 - q1*(p1 - p0)/(q1 - q0)
+            if abs(p - p1) < tol:
+                return p
+            p0 = p1
+            q0 = q1
+            p1 = p
+            q1 = func(*((p1,) + args))
+    msg = "Failed to converge after %d iterations, value is %s" % (maxiter, p)
+    raise RuntimeError(msg)
+
+
 def bisect(f, a, b, args=(),
            xtol=_xtol, rtol=_rtol, maxiter=_iter,
            full_output=False, disp=True):
@@ -91,6 +187,7 @@
     r = _zeros._bisect(f,a,b,xtol,maxiter,args,full_output,disp)
     return results_c(full_output, r)
 
+
 def ridder(f, a, b, args=(),
            xtol=_xtol, rtol=_rtol, maxiter=_iter,
            full_output=False, disp=True):
@@ -160,6 +257,7 @@
     r = _zeros._ridder(f,a,b,xtol,maxiter,args,full_output,disp)
     return results_c(full_output, r)
 
+
 def brentq(f, a, b, args=(),
            xtol=_xtol, rtol=_rtol, maxiter=_iter,
            full_output=False, disp=True):
@@ -262,6 +360,7 @@
     r = _zeros._brentq(f,a,b,xtol,maxiter,args,full_output,disp)
     return results_c(full_output, r)
 
+
 def brenth(f, a, b, args=(),
            xtol=_xtol, rtol=_rtol, maxiter=_iter,
            full_output=False, disp=True):



More information about the Scipy-svn mailing list