/usr/lib/python2.7/dist-packages/Scientific/Functions/FindRoot.py is in python-scientific 2.9.4-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 | # 'Safe' Newton-Raphson for numerical root-finding
#
# Written by Scott M. Ransom <ransom@cfa.harvard.edu>
# last revision: 14 Nov 98
#
# Cosmetic changes by Konrad Hinsen <hinsen@cnrs-orleans.fr>
# last revision: 2006-11-23
#
"""
Newton-Raphson for numerical root finding
Example::
>>>from Scientific.Functions.FindRoot import newtonRaphson
>>>from Scientific.N import pi, sin, cos
>>>def func(x):
>>> return (2*x*cos(x) - sin(x))*cos(x) - x + pi/4.0
>>>newtonRaphson(func, 0.0, 1.0, 1.0e-12)
yields 0.952847864655.
"""
from FirstDerivatives import DerivVar
def newtonRaphson(function, lox, hix, xacc):
"""
X{Newton-Raphson} algorithm for X{root} finding.
The algorithm used is a safe version of Newton-Raphson
(see page 366 of Numerical Recipes in C, 2ed).
@param function: function of one numerical variable that
uses only those operations that are defined for DerivVar
objects in the module L{Scientific.Functions.FirstDerivatives}
@type function: callable
@param lox: lower limit of the search interval
@type lox: C{float}
@param hix: upper limit of the search interval
@type hix: C[{loat}
@param xacc: requested absolute precision of the root
@type xacc: C{float}
@returns: a root of function between lox and hix
@rtype: C{float}
@raises ValueError: if C{function(lox)} and C{function(hix)} have
the same sign, or if there is no convergence after 500 iterations
"""
maxit = 500
tmp = function(DerivVar(lox))
fl = tmp[0]
tmp = function(DerivVar(hix))
fh = tmp[0]
if ((fl > 0.0 and fh > 0.0) or (fl < 0.0 and fh < 0.0)):
raise ValueError("Root must be bracketed")
if (fl == 0.0): return lox
if (fh == 0.0): return hix
if (fl < 0.0):
xl=lox
xh=hix
else:
xh=lox
xl=hix
rts=0.5*(lox+hix)
dxold=abs(hix-lox)
dx=dxold
tmp = function(DerivVar(rts))
f = tmp[0]
df = tmp[1][0]
for j in range(maxit):
if ((((rts-xh)*df-f)*((rts-xl)*df-f) > 0.0)
or (abs(2.0*f) > abs(dxold*df))):
dxold=dx
dx=0.5*(xh-xl)
rts=xl+dx
if (xl == rts): return rts
else:
dxold=dx
dx=f/df
temp=rts
rts=rts-dx
if (temp == rts): return rts
if (abs(dx) < xacc): return rts
tmp = function(DerivVar(rts))
f = tmp[0]
df = tmp[1][0]
if (f < 0.0):
xl=rts
else:
xh=rts
raise ValueError("Maximum number of iterations exceeded")
# Test code
if __name__ == '__main__':
from Scientific.Numeric import pi, sin, cos
def _func(x):
return ((2*x*cos(x) - sin(x))*cos(x) - x + pi/4.0)
_r = newtonRaphson(_func, 0.0, 1.0, 1.0e-12)
_theo = 0.9528478646549419474413332
print ''
print 'Finding the root (between 0.0 and 1.0) of:'
print ' (2*x*cos(x) - sin(x))*cos(x) - x + pi/4 = 0'
print ''
print 'Safe-style Newton-Raphson gives (xacc = 1.0e-12) =', _r
print 'Theoretical result (correct to all shown digits) = %15.14f' % _theo
print ''
|