find_root

def find_root (   self,   lower_bracket,   upper_bracket,   start,   value,   trace = False   ) solve f(x)==value very efficiently, with explicit knowledge of derivatives of the function find_root solves by iterated inverse quadratic extrapolation for a solution to f(x)=y. It includes checks against bad convergence, so it should never be able to fail. Unlike typical secant method or fancier Brent's method finders, this does not depend in any strong wasy on the brackets, unless the finder has to resort to successive approximations to close in on a root. Often, it is possible to make the brackets equal to the domain of the function, if there is any clue as to where the root lies, as given by the parameter start. Parameters: lower_bracket  the lower bound which is ever permitted in the search. upper_bracket  the upper bound. Note that the function have opposite signs at these two points start  an initial guess for where to start the search value  the target value of the function trace  print debugging info to sys.stderr if True Definition at line 157 of file C2Functions.py. : "solve f(x)==value very efficiently, with explicit knowledge of derivatives of the function" # can use local f(x)=a*x**2 + b*x + c # and solve quadratic to find root, then iterate ftol=(5e-16*abs(value)+2e-308); xtol=(5e-16*(abs(upper_bracket)+abs(lower_bracket))+2e-308); root=start #start looking in the middle cupper, yp, ypp=self.value_with_derivatives(upper_bracket) cupper -= value if abs(cupper) < ftol: return upper_bracket clower, yp, ypp=self.value_with_derivatives(lower_bracket) clower -= value if abs(clower) < ftol: return lower_bracket; if cupper*clower >0: # argh, no sign change in here! raise C2Exception("unbracketed root in find_root at xlower=%g, xupper=%g, bailing" % (lower_bracket, upper_bracket)) delta=upper_bracket-lower_bracket #first error step c, b, ypp=self.value_with_derivatives(root) # compute initial values c -= value while abs(delta) > xtol: # can allow quite small steps, since we have exact derivatives! a=ypp/2 #second derivative is 2*a disc=b*b-4*a*c if disc >= 0: if b>=0: q=-0.5*(b+math.sqrt(disc)) else: q=-0.5*(b-math.sqrt(disc)) if q*q > abs(a*c): delta=c/q #since x1=q/a, x2=c/q, x1/x2=q^2/ac, this picks smaller step else: delta=q/a root+=delta; if disc < 0 or root<=lower_bracket or root>=upper_bracket or abs(delta) >= 0.5*(upper_bracket-lower_bracket): #if we jump out of the bracket, or if the step is not smaller than half the bracket, bisect root=0.5*(lower_bracket+upper_bracket) delta=upper_bracket-lower_bracket c, b, ypp=self.value_with_derivatives(root) #get local curvature and value c -= value if trace: import sys print >> sys.stderr, "find_root x, dx, c, b, a", ( (5*"%10.4g ") % (root, delta, c, b, ypp/2) ) if abs(c) < ftol or abs(c) < xtol*abs(b): return root #got it close enough # now, close in bracket on whichever side this still brackets if c*clower < 0.0: cupper=c upper_bracket=root else: clower=c lower_bracket=root return root

---