Python mpmath.findroot函数代码示例

 def _findRootAttempt(self, startingEne, multipler, type):
     try:
         fun = lambda e: self._getDet(e, multipler)
         if type == 'muller':
             return complex(mpmath.findroot(fun, startingEne, solver='muller', maxsteps=10000, tol=2.16840434497100886801e-19))
         else:
             return complex(mpmath.findroot(fun, startingEne[0], solver='secant'))
     except ValueError:
         return None

def getZero(Func, init=(1.+1.j, 1.1+1.j, 1.+1.1j), solver='muller', ftol=1.e-3, tol=1.e-3, maxsteps=600):
  import mpmath as mp
  w0 = complex(float('nan'), float('nan'))
  try:
     w0 = mp.findroot(lambda w : Func(complex(w)), init, solver=solver, maxsteps=maxsteps, ftol=ftol, tol=tol)
  except:
   try:
     w0 = mp.findroot(lambda w : Func(complex(w)), init[0], solver='halley', maxsteps=maxsteps, ftol=ftol, tol=tol)
   except:   
     print "Could not find zero"

  return complex(w0)

 def draw_waiting_time(self, avail_devices, failed_devices):
     self.curr_avail_devices = avail_devices
     self.curr_failed_devices = failed_devices
     self.curr_uniform_variate = random.uniform(0,1)
     
     while self.curr_uniform_variate == 0:
             self.curr_uniform_variate = random.uniform(0,1)
     
     try:
         wait_time = findroot(self.func, [0,100], solver='secant')
     except:
         wait_time = findroot(self.func, [0,100], solver='secant', maxsteps=1000)
         #wait_time = findroot(self.func, [0, 87600] , maxsteps=1000, solver='secant', verify=False)
         
     return abs(wait_time)

def nodes(n):
    left  = mp.mpf(0)
    right = mp.mpf(n+(n-1)*sqrt(n)*1.01)

    i = 2
    factor = 2

    while True:
        l = [ (x,mp.laguerre(n,0,x)) for x in mp.linspace(left,right,n*factor**i)]

        intervals = []
        for j in range(len(l)-1):
            prod = l[j][1]*l[j+1][1]
            if prod < 0:
                intervals.append([l[j][0], l[j+1][0]])

        if len(intervals) == n:
            break
        i += 1

    roots = []
    f = lambda x: mp.laguerre(n,0,x)
    for ab in intervals:
        a,b = ab
        try:
            z = mp.findroot(f, (a, b), tol=1e-50, solver='bisect')
        except:
            z = bisect(f, a, b, tol=1e-50)
        roots.append( z )

    return roots

    def update_mix(self):
        if self.esolute.shape != self.esolvent.shape:
            return
        
        eeff=np.empty(self.esolute.shape, dtype=complex)
        for i in range(len(eeff)):
            e1 = self.esolute[i]#.real    
            em = self.esolvent[i]#.real   #THIS IS GENERALLY THE SOLVENT
            
            partialfunc = partial(self.fill_func_dielectric, em=em, e1=e1, v=self.Vfrac, w=self.w)
            
            
            eeff[i] = findroot(partialfunc, 
                               (0.5, 1, 2),  #This is a root that we are guessing for 3 params
                               solver='muller',
                               )            
            
            #eeff[i] = scipy.optimize.newton(self.fill_func_dielectric, 
                                            #self.Root,
                                            #fprime=None, 
                                            #args=(em, e1, self.Vfrac, self.w),
                                            #tol=.0001, 
                                            #maxiter=1000)
#			self.Root=eeff[i]  #Reset the root
        self.mixedarray=eeff

 def long_interval(w, n, t):
     """
     Return the interval of integration for longitudinal impedance, which
     is of the form [0, root, inf], where root satisfies
     d_l(root, w, n, t) = 0.
     
     """
     if w < mp.sqrt(1+n):
         return [0, mp.inf]
     if n == 0:
         guesses = [4, 6, 8, 10]
     else:
         guesses = [4, 6, 8, 10, 12, 14]
     int_range = [0, mp.inf]    
     for guess in guesses:
         try:
             root = mp.findroot(lambda z: BiMax.d_l(z, w, n, t), guess)
         except ValueError:
             continue
         unique = True
         for z in int_range:
             if mp.fabs(z - mp.fabs(root)) < 1e-4:
                 unique = False
         if unique:
             int_range += [mp.fabs(root)]
     int_range = np.sort(int_range)
     return int_range

    def generate_zeta(self):
        import mpmath
        import sys
        root_list		= {}
        for s in [complex(0.5,t) for t in range(10,100)]:
            try:
                root_complex	= mpmath.findroot(mpmath.zeta, s)
                if self.debug:
                    sys.stderr.write("[NOTICE] Found nontrivial zeta zero at %s.\n" % str(root_complex))
                root_imag	= float(root_complex.imag)
                try:
                    root_list[root_imag]
                    if self.debug:
                        sys.stderr.write("[NOTICE] Already found this zero, skipping.\n")
                    next
                except KeyError:
                    if root_imag > 1:
                        yield root_imag
                        root_list[root_imag] = True
                    else:
                        if self.debug:
                            sys.stderr.write("[NOTICE] Irrelevant zero, skipping.\n")
            except ValueError, e:
#DEBUG#                sys.stderr.write("Exception at s=%s: %s.\n" % (str(s), e.args[0])) #DEBUG#
                pass
            except OverflowError, e:
#DEBUG#                sys.stderr.write("Exception at s=%s: %s.\n" % (str(s), e.args[0])) #DEBUG#
                pass

 def int_interval(wrel, n, t, k):
     """
     find out the almost-singular point and
     divide the integration integrals.
     
     """
     if wrel < 1 or wrel > 1.3:
         return [0, mp.inf]
     else:
         guesses = [4, 6, 8, 10, 12, 14]
     wc = wrel * mp.sqrt(1+n)
     int_range = [0, mp.inf]
     for guess in guesses:
         try:
             root = mp.findroot(lambda zc: MaxKappa.e_l(zc, wc, n, t, k), guess)
         except ValueError:
             continue
         unique = True
         for z in int_range:
             if mp.fabs(z - mp.fabs(root)) < 1e-4:
                 unique = False
         if unique:
             int_range += [mp.fabs(root)]
     int_range = np.sort(int_range)
     return int_range

	def mini(r):
		nonlocal txi, tv
		xi = [r]
		y = f(r)
		v = r*f(r) + mpmath.quad(f, [r, mpmath.inf])
		if verbose:
			print('Trying r={0} (v={1})'.format(r, v), file=sys.stderr)
		for i in itertools.count():
			xm1 = xi[i]
			h = v / xm1
			y += h
			if y >= maximum or mpmath.almosteq(y, maximum, abs_eps=mpmath.mp.eps * 2**10):
				break
			# We solve for x via secant method instead of using f's inverse.
			x = mpmath.findroot(lambda x: f(x) - y, xm1)
			xi.append(x)
		xi.append(mpmath.mpf())
		if len(xi) == nseg:
			if mpmath.almosteq(y, maximum, abs_eps=mpmath.mp.eps * 2**10):
				txi, tv = xi[::-1], v
				return 0
			# If y > maximum, then v is too large, which means r is too far
			# left, so we want to return a negative value. The opposite holds
			# true when y < maximum.
			return maximum - y
		return len(xi) - nseg + h*mpmath.sign(len(xi) - nseg)

def solve(f, x0=None, nseg=128, verbose=False):
	"""Find r, v, and the x coordinates for f."""
	r = x0
	if r is None:
		if verbose:
			print('Calculating initial guess... ', end='', file=sys.stderr)
		# The area we seek is nseg equal-area rectangles surrounding f(x), not
		# f(x) itself, but we can get a good approximation from it.
		v = mpmath.quad(f, [0, mpmath.inf]) / nseg
		r = mpmath.findroot(lambda x: x*f(x) + mpmath.quad(f, [x, mpmath.inf]) - v, x0=1, x1=100, maxsteps=100)
		if verbose:
			print(r, file=sys.stderr)
	# We know that f(0) is the maximum because f must decrease monotonically.
	maximum = f(0)
	txi = []
	tv = mpmath.mpf()
	def mini(r):
		nonlocal txi, tv
		xi = [r]
		y = f(r)
		v = r*f(r) + mpmath.quad(f, [r, mpmath.inf])
		if verbose:
			print('Trying r={0} (v={1})'.format(r, v), file=sys.stderr)
		for i in itertools.count():
			xm1 = xi[i]
			h = v / xm1
			y += h
			if y >= maximum or mpmath.almosteq(y, maximum, abs_eps=mpmath.mp.eps * 2**10):
				break
			# We solve for x via secant method instead of using f's inverse.
			x = mpmath.findroot(lambda x: f(x) - y, xm1)
			xi.append(x)
		xi.append(mpmath.mpf())
		if len(xi) == nseg:
			if mpmath.almosteq(y, maximum, abs_eps=mpmath.mp.eps * 2**10):
				txi, tv = xi[::-1], v
				return 0
			# If y > maximum, then v is too large, which means r is too far
			# left, so we want to return a negative value. The opposite holds
			# true when y < maximum.
			return maximum - y
		return len(xi) - nseg + h*mpmath.sign(len(xi) - nseg)
	r = mpmath.findroot(mini, r)
	assert len(txi) == nseg
	if verbose:
		print('Done calculating r, v, x[i].', file=sys.stderr)
	return r, tv, txi

    def _eval_evalf(self, prec):
        """Evaluate this complex root to the given precision."""
        with workprec(prec):
            g = self.poly.gen
            if not g.is_Symbol:
                d = Dummy('x')
                func = lambdify(d, self.expr.subs({g: d}), "mpmath")
            else:
                func = lambdify(g, self.expr, "mpmath")

            try:
                interval = self.interval
            except DomainError:
                return super()._eval_evalf(prec)

            while True:
                if self.is_extended_real:
                    a = mpf(str(interval.a))
                    b = mpf(str(interval.b))
                    if a == b:
                        root = a
                        break
                    x0 = mpf(str(interval.center))
                else:
                    ax = mpf(str(interval.ax))
                    bx = mpf(str(interval.bx))
                    ay = mpf(str(interval.ay))
                    by = mpf(str(interval.by))
                    x0 = mpc(*map(str, interval.center))
                    if ax == bx and ay == by:
                        root = x0
                        break

                try:
                    root = findroot(func, x0)
                    # If the (real or complex) root is not in the 'interval',
                    # then keep refining the interval. This happens if findroot
                    # accidentally finds a different root outside of this
                    # interval because our initial estimate 'x0' was not close
                    # enough. It is also possible that the secant method will
                    # get trapped by a max/min in the interval; the root
                    # verification by findroot will raise a ValueError in this
                    # case and the interval will then be tightened -- and
                    # eventually the root will be found.
                    if self.is_extended_real:
                        if (a <= root <= b):
                            break
                    elif (ax <= root.real <= bx and ay <= root.imag <= by
                          and (interval.ay > 0 or interval.by < 0)):
                        break
                except (ValueError, UnboundLocalError):
                    pass
                self.refine()
                interval = self.interval

        return ((Float._new(root.real._mpf_, prec) if not self.is_imaginary else 0) +
                I*Float._new(root.imag._mpf_, prec))

    def _eval_evalf(self, prec):
        """Evaluate this complex root to the given precision. """
        with workprec(prec):
            g = self.poly.gen
            if not g.is_Symbol:
                d = Dummy('x')
                func = lambdify(d, self.expr.subs(g, d))
            else:
                func = lambdify(g, self.expr)

            interval = self._get_interval()
            if not self.is_real:
                # For complex intervals, we need to keep refining until the
                # imaginary interval is disjunct with other roots, that is,
                # until both ends get refined.
                ay = interval.ay
                by = interval.by
                while interval.ay == ay or interval.by == by:
                    interval = interval.refine()

            while True:
                if self.is_real:
                    x0 = mpf(str(interval.center))
                else:
                    x0 = mpc(*map(str, interval.center))
                try:
                    root = findroot(func, x0, verify=False)
                    # If the (real or complex) root is not in the 'interval',
                    # then keep refining the interval. This happens if findroot
                    # accidentally finds a different root outside of this
                    # interval because our initial estimate 'x0' was not close
                    # enough.
                    if self.is_real:
                        a = mpf(str(interval.a))
                        b = mpf(str(interval.b))
                        if a == b:
                            root = a
                            break
                        if not (a < root < b):
                            raise ValueError("Root not in the interval.")
                    else:
                        ax = mpf(str(interval.ax))
                        bx = mpf(str(interval.bx))
                        ay = mpf(str(interval.ay))
                        by = mpf(str(interval.by))
                        if ax == bx and ay == by:
                            root = ax + S.ImaginaryUnit*by
                            break
                        if not (ax < root.real < bx and ay < root.imag < by):
                            raise ValueError("Root not in the interval.")
                except ValueError:
                    interval = interval.refine()
                    continue
                else:
                    break

        return Float._new(root.real._mpf_, prec) + I*Float._new(root.imag._mpf_, prec)

 def solver(args):
     equation, lower_boundary, upper_boundary, initial_guess, f, use_mpmath = args
     if not use_mpmath:
         warnings.filterwarnings('ignore')
         solution = root(equation, initial_guess, args=f, method='hybr')
         solution = solution.x[0]
         warnings.resetwarnings()
     else:
         equation_wrap = partial(equation, f=f)
         try:
             solution = mp.findroot(equation_wrap, (lower_boundary, upper_boundary),
                                    maxsteps=1000, solver='anderson', tol=5e-16)
         except ValueError as err:
             print err
             print 'Lowering tolerance to 5e-6'
             solution = mp.findroot(equation_wrap, (lower_boundary, upper_boundary),
                                    maxsteps=1000, solver='anderson', tol=5e-6)
         solution = np.float(solution)
     return solution

def besselJZeros(m, a, b):
  require_mpmath()
  if not hasattr(mpmath, 'besseljzero'):
    besseljn = lambda x: mpmath.besselj(m, x)
    results = [mpmath.findroot(besseljn, mpmath.pi*(kp - 1./4 + 0.5*m)) for kp in range(a, b+1)]
  else:
    results = [mpmath.besseljzero(m, i) for i in xrange(a, b+1)]
  # Check that we haven't found double roots or missed a root. All roots should be separated by approximately pi
  assert all([0.6*mpmath.pi < (b - a) < 1.4*mpmath.pi for a, b in zip(results[:-1], results[1:])]), "Separation of Bessel zeros was incorrect."
  return results

def printAnaRoots(V, fndStates):
  print "\n\nFor V=" + str(V)
  def denumWrap(k):
    return denum(k,a,V)
  for state in fndStates:
    try:
      root = mpm.findroot(denumWrap, state)
      diff = abs(root-state)
      print str(root.real) + "\t" + str(root.imag) + "\t" + str(diff)
    except ValueError:
      print "Not Fnd"