Python optimize.fsolve方法代码示例

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def _fitstart(self, data):
        g1 = _skew(data)
        g2 = _kurtosis(data)

        def func(x):
            a, b = x
            sk = 2*(b-a)*np.sqrt(a + b + 1) / (a + b + 2) / np.sqrt(a*b)
            ku = a**3 - a**2*(2*b-1) + b**2*(b+1) - 2*a*b*(b+2)
            ku /= a*b*(a+b+2)*(a+b+3)
            ku *= 6
            return [sk-g1, ku-g2]
        a, b = optimize.fsolve(func, (1.0, 1.0))
        return super(beta_gen, self)._fitstart(data, args=(a, b)) 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def _xinf_ND(xdot,x0,args=(),xddot=None,xtol=1.49012e-8):
    """Private function for wrapping the fsolving for x_infinity
    for a variable x in N dimensions"""
    try:
        result = fsolve(xdot,x0,args,fprime=xddot,xtol=xtol,full_output=1)
    except (ValueError, TypeError, OverflowError):
        xinf_val = NaN
    except:
        print "Error in fsolve:", sys.exc_info()[0], sys.exc_info()[1]
        xinf_val = NaN
    else:
        if result[2] in (1,2,3): #,4,5):
            # 4,5 means "not making good progress" (see fsolve docstring)
            xinf_val = float(result[0])
        else:
            xinf_val = NaN
    return xinf_val 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def computeShiftedGears(m, alpha, t1, t2, x1, x2):
    """Summary

    Args:
        m (float): common module of both gears [length]
        alpha (float): pressure-angle [rad]
        t1 (int): number of teeth of gear1
        t2 (int): number of teeth of gear2
        x1 (float): relative profile-shift of gear1
        x2 (float): relative profile-shift of gear2

    Returns:
        (float, float): distance between gears [length], pressure angle of the assembly [rad]
    """
    def inv(x): return np.tan(x) - x
    inv_alpha_w = inv(alpha) + 2 * np.tan(alpha) * (x1 + x2) / (t1 + t2)

    def root_inv(x): return inv(x) - inv_alpha_w
    alpha_w = opt.fsolve(root_inv, 0.)
    dist = m * (t1 + t2) / 2 * np.cos(alpha) / np.cos(alpha_w)
    return dist, alpha_w 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def calculate_psi_goal(psi_baseline, Au_baseline, Au_goal, deltaz,
                       c_baseline, c_goal):
    """Find the value for psi that has the same location w on the upper
    surface of the goal as psi_baseline on the upper surface of the
    baseline"""

    def integrand(psi_baseline, Au, deltaz, c):
        return c*np.sqrt(1 + dxi_u(psi_baseline, Au, deltaz/c)**2)

    def equation(psi_goal, L_baseline, Au_goal, deltaz, c):
        y, err = quad(integrand, 0, psi_goal, args=(Au_goal, deltaz, c))
        return y - L_baseline

    L_baseline, err = quad(integrand, 0, psi_baseline,
                           args=(Au_baseline, deltaz, c_baseline))
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        y = fsolve(equation, psi_baseline, args=(L_baseline, Au_goal, deltaz,
                                                 c_goal))
    return y[0] 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def interp_isentrope(interp, pressures, entropy, T_guess):
    def _deltaS(args, S, P):
        T = args[0]
        return interp(P, T)[0] - S
    
    sol = [T_guess]
    temperatures = np.empty_like(pressures)
    for i, P in enumerate(pressures):
        sol = fsolve(_deltaS, sol, args=(entropy, P))
        temperatures[i] = sol[0]

    return temperatures

# Define function to self consistently calculate depth and gravity profiles
# from pressure and density profiles. 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def equilibrium_pressure(minerals, stoichiometry, temperature, pressure_initial_guess=1.e5):
    """
    Given a list of minerals, their reaction stoichiometries
    and a temperature of interest, compute the
    equilibrium pressure of the reaction.

    Parameters
    ----------
    minerals : list of minerals
        List of minerals involved in the reaction.

    stoichiometry : list of floats
        Reaction stoichiometry for the minerals provided.
        Reactants and products should have the opposite signs [mol]

    temperature : float
        Temperature of interest [K]

    pressure_initial_guess : optional float
        Initial pressure guess [Pa]

    Returns
    -------
    pressure : float
        The equilibrium pressure of the reaction [Pa]
    """
    def eqm(P, T):
        gibbs = 0.
        for i, mineral in enumerate(minerals):
            mineral.set_state(P[0], T)
            gibbs = gibbs + mineral.gibbs * stoichiometry[i]
        return gibbs

    pressure = fsolve(eqm, [pressure_initial_guess], args=(temperature))[0]

    return pressure 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def equilibrium_temperature(minerals, stoichiometry, pressure, temperature_initial_guess=1000.):
    """
    Given a list of minerals, their reaction stoichiometries
    and a pressure of interest, compute the
    equilibrium temperature of the reaction.

    Parameters
    ----------
    minerals : list of minerals
        List of minerals involved in the reaction.

    stoichiometry : list of floats
        Reaction stoichiometry for the minerals provided.
        Reactants and products should have the opposite signs [mol]

    pressure : float
        Pressure of interest [Pa]

    temperature_initial_guess : optional float
        Initial temperature guess [K]

    Returns
    -------
    temperature : float
        The equilibrium temperature of the reaction [K]
    """
    def eqm(T, P):
        gibbs = 0.
        for i, mineral in enumerate(minerals):
            mineral.set_state(P, T[0])
            gibbs = gibbs + mineral.gibbs * stoichiometry[i]
        return gibbs

    temperature = fsolve(eqm, [temperature_initial_guess], args=(pressure))[0]

    return temperature 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def _digammainv(y):
    # Inverse of the digamma function (real positive arguments only).
    # This function is used in the `fit` method of `gamma_gen`.
    # The function uses either optimize.fsolve or optimize.newton
    # to solve `sc.digamma(x) - y = 0`.  There is probably room for
    # improvement, but currently it works over a wide range of y:
    #    >>> y = 64*np.random.randn(1000000)
    #    >>> y.min(), y.max()
    #    (-311.43592651416662, 351.77388222276869)
    #    x = [_digammainv(t) for t in y]
    #    np.abs(sc.digamma(x) - y).max()
    #    1.1368683772161603e-13
    #
    _em = 0.5772156649015328606065120
    func = lambda x: sc.digamma(x) - y
    if y > -0.125:
        x0 = np.exp(y) + 0.5
        if y < 10:
            # Some experimentation shows that newton reliably converges
            # must faster than fsolve in this y range.  For larger y,
            # newton sometimes fails to converge.
            value = optimize.newton(func, x0, tol=1e-10)
            return value
    elif y > -3:
        x0 = np.exp(y/2.332) + 0.08661
    else:
        x0 = 1.0 / (-y - _em)

    value, info, ier, mesg = optimize.fsolve(func, x0, xtol=1e-11,
                                             full_output=True)
    if ier != 1:
        raise RuntimeError("_digammainv: fsolve failed, y = %r" % y)

    return value[0]


## Gamma (Use MATLAB and MATHEMATICA (b=theta=scale, a=alpha=shape) definition)

## gamma(a, loc, scale)  with a an integer is the Erlang distribution
## gamma(1, loc, scale)  is the Exponential distribution
## gamma(df/2, 0, 2) is the chi2 distribution with df degrees of freedom. 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def create_test_fn(gen, tmax, dist_pts):
    """Create integration test function using the supplied generator.

    The test function will return a dictionary of the two closest
    points (keys 1 and 2) mapping to their respective distances and
    pointset index positions.
    """

    gen.set(tdata=[0,tmax])

    def Tn_integ(ic):
        gen.set(ics=ic)
        try:
            test_traj = gen.compute('test')
        except:
            print "Problem integrating test trajectory at i.c. ", ic
            raise
        test_pts = test_traj.sample(coords=dist_pts.all_pts.coordnames)
        # distance of endpoint to pointset
        try:
            d_info = dist_pts(test_pts[-1], use_norm=True, minmax=['min'])
        except ValueError:
            # This error happens when fsolve tries initial conditions that
            # break the integrator
            return (_num_inf,NaN)
        pos = d_info['min'][1]['pos']
        return (test_pts[-1]-dist_pts.all_pts[pos], pos)
    return Tn_integ 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def create_test_fn_with_events(gen, tmax, dist_pts, iso_ev, other_evnames, pars_to_vars):
    """Create integration test function using the supplied generator,
    assuming it contains isochron-related events.

    The test function will return a dictionary of the two closest
    points (keys 1 and 2) mapping to their respective distances and
    pointset index positions.
    """

    def Tn_integ(ic):
        gen.set(ics=ic, tdata=[0,tmax])
        try:
            test_traj = gen.compute('test')
        except:
            print "Problem integrating test trajectory at i.c. ", ic
            raise
        test_pts = test_traj.sample(coords=dist_pts.all_pts.coordnames)
        # distance of endpoint to pointset
        try:
            d_info = dist_pts(test_pts[-1], use_norm=True, minmax=['min'])
        except ValueError:
            # This error happens when fsolve tries initial conditions that
            # break the integrator
            return (_num_inf,NaN)
        # refine min position using isochron-related events
        q=test_pts[-1]
        perp_ev, t_ev = _find_min_pt(gen, q,
                                     d_info['min'][1]['pos'], dist_pts.all_pts,
                                     pars_to_vars, iso_ev, other_evnames)
        ev_pt = perp_ev[0][q.coordnames]
        # different return format to version w/o events
        return (test_pts[-1]-ev_pt, t_ev, ev_pt)
    return Tn_integ 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def _xinf_1D(xdot,x0,args=(),xddot=None,xtol=1.49012e-8):
    """Private function for wrapping the solving for x_infinity
    for a variable x in 1 dimension"""
    try:
        if xddot is None:
            xinf_val = float(fsolve(xdot,x0,args,xtol=xtol))
        else:
            xinf_val = float(newton_meth(xdot,x0,fprime=xddot,args=args))
    except RuntimeError:
        xinf_val = NaN
    return xinf_val 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def test_pressure_network_no_gradient(self):
        """fsolve without gradient, equal pipes -> equal flows"""
        k = np.ones(4) * 0.5
        Qtot = 4
        initial_guess = array([2., 0., 2., 0.])
        final_flows, info, ier, mesg = optimize.fsolve(
            pressure_network, initial_guess, args=(Qtot, k),
            full_output=True)
        assert_array_almost_equal(final_flows, np.ones(4))
        assert_(ier == 1, mesg) 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def test_pressure_network_with_gradient(self):
        """fsolve with gradient, equal pipes -> equal flows"""
        k = np.ones(4) * 0.5
        Qtot = 4
        initial_guess = array([2., 0., 2., 0.])
        final_flows = optimize.fsolve(
            pressure_network, initial_guess, args=(Qtot, k),
            fprime=pressure_network_jacobian)
        assert_array_almost_equal(final_flows, np.ones(4)) 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def test_wrong_shape_func_callable(self):
        """The callable 'func' has no '__name__' attribute."""
        func = ReturnShape(1)
        # x0 is a list of two elements, but func will return an array with
        # length 1, so this should result in a TypeError.
        x0 = [1.5, 2.0]
        assert_raises(TypeError, optimize.fsolve, func, x0) 

# 需要导入模块: from scipy import optimize [as 别名]
# 或者: from scipy.optimize import fsolve [as 别名]
def test_wrong_shape_func_function(self):
        # x0 is a list of two elements, but func will return an array with
        # length 1, so this should result in a TypeError.
        x0 = [1.5, 2.0]
        assert_raises(TypeError, optimize.fsolve, dummy_func, x0, args=((1,),))