Python integrate.quad方法代码示例

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def get_heat_capacity(temp, td):
        # credits to Dr. Joseph Montoya, LBNL
        t_ratio = temp / td

        def integrand(x):
            return (x ** 4 * pd.np.exp(x)) / (pd.np.exp(x) - 1) ** 2

        if isinstance(t_ratio, int) or isinstance(t_ratio, float):
            cv_p = 9 * R * (t_ratio ** 3) * quad(integrand, 0, t_ratio ** -1)[0]
        else:
            cv_p = []
            for i in range(len(t_ratio)):
                cv_i = (
                    9 * R * (t_ratio[i] ** 3) * quad(integrand, 0, t_ratio[i] ** -1)[0]
                )
                cv_p = pd.np.append(cv_p, cv_i)
        return cv_p * 5 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def _entropy(self, *args):
        def integ(x):
            val = self._pdf(x, *args)
            return entr(val)

        # upper limit is often inf, so suppress warnings when integrating
        olderr = np.seterr(over='ignore')
        h = integrate.quad(integ, self.a, self.b)[0]
        np.seterr(**olderr)

        if not np.isnan(h):
            return h
        else:
            # try with different limits if integration problems
            low, upp = self.ppf([1e-10, 1. - 1e-10], *args)
            if np.isinf(self.b):
                upper = upp
            else:
                upper = self.b
            if np.isinf(self.a):
                lower = low
            else:
                lower = self.a
            return integrate.quad(integ, lower, upper)[0] 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def mvstdtprob(a, b, R, df, ieps=1e-5, quadkwds=None, mvstkwds=None):
    '''probability of rectangular area of standard t distribution

    assumes mean is zero and R is correlation matrix

    Notes
    -----
    This function does not calculate the estimate of the combined error
    between the underlying multivariate normal probability calculations
    and the integration.

    '''
    kwds = dict(args=(a,b,R,df), epsabs=1e-4, epsrel=1e-2, limit=150)
    if not quadkwds is None:
        kwds.update(quadkwds)
    #print kwds
    res, err = integrate.quad(funbgh2, *chi.ppf([ieps,1-ieps], df),
                          **kwds)
    prob = res * bghfactor(df)
    return prob

#written by Enzo Michelangeli, style changes by josef-pktd
# Student's T random variable 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def cdf(self):
        """
        Returns the cumulative distribution function evaluated at the support.

        Notes
        -----
        Will not work if fit has not been called.
        """
        _checkisfit(self)
        kern = self.kernel
        if kern.domain is None: # TODO: test for grid point at domain bound
            a,b = -np.inf,np.inf
        else:
            a,b = kern.domain
        func = lambda x,s: kern.density(s,x)

        support = self.support
        support = np.r_[a,support]
        gridsize = len(support)
        endog = self.endog
        probs = [integrate.quad(func, support[i - 1], support[i],
                    args=endog)[0] for i in range(1, gridsize)]
        return np.cumsum(probs) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def entropy(self):
        """
        Returns the differential entropy evaluated at the support

        Notes
        -----
        Will not work if fit has not been called. 1e-12 is added to each
        probability to ensure that log(0) is not called.
        """
        _checkisfit(self)

        def entr(x,s):
            pdf = kern.density(s,x)
            return pdf*np.log(pdf+1e-12)

        kern = self.kernel

        if kern.domain is not None:
            a, b = self.domain
        else:
            a, b = -np.inf, np.inf
        endog = self.endog
        #TODO: below could run into integr problems, cf. stats.dist._entropy
        return -integrate.quad(entr, a, b, args=(endog,))[0] 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def length(self, t0=0, t1=1, error=LENGTH_ERROR, min_depth=LENGTH_MIN_DEPTH):
        """Calculate the length of the path up to a certain position"""
        if t0 == 0 and t1 == 1:
            if self._length_info['bpoints'] == self.bpoints() \
                    and self._length_info['error'] >= error \
                    and self._length_info['min_depth'] >= min_depth:
                return self._length_info['length']

        # using scipy.integrate.quad is quick
        if _quad_available:
            s = quad(lambda tau: abs(self.derivative(tau)), t0, t1,
                            epsabs=error, limit=1000)[0]
        else:
            s = segment_length(self, t0, t1, self.point(t0), self.point(t1),
                               error, min_depth, 0)

        if t0 == 0 and t1 == 1:
            self._length_info['length'] = s
            self._length_info['bpoints'] = self.bpoints()
            self._length_info['error'] = error
            self._length_info['min_depth'] = min_depth
            return self._length_info['length']
        else:
            return s 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def calculate_c_baseline(c_L, Au_C, Au_L, deltaz):
    """Equations in the New_CST.pdf. Calculates the upper chord in order for
       the cruise and landing airfoils ot have the same length."""

    def integrand(psi, Au, delta_xi):
        return np.sqrt(1 + dxi_u(psi, Au, delta_xi)**2)

    def f(c_C):
        """Function dependent of c_C and that outputs c_C."""
        y_C, err = quad(integrand, 0, 1, args=(Au_C, deltaz/c_C))
        y_L, err = quad(integrand, 0, 1, args=(Au_L, deltaz/c_L))
        return c_L*y_L/y_C
    c_C = optimize.fixed_point(f, [c_L])
    # In case the calculated chord is really close to the original, but the
    # algorithm was not able to make them equal
    if abs(c_L - c_C) < 1e-7:
        return c_L
    # The output is an array so it needs the extra [0]
    return c_C[0] 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [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 integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def verify_gauss_quad(root_func, eval_func, weight_func, a, b, N,
                      rtol=1e-15, atol=1e-14):
        # this test is copied from numpy's TestGauss in test_hermite.py
        x, w, mu = root_func(N, True)

        n = np.arange(N)
        v = eval_func(n[:,np.newaxis], x)
        vv = np.dot(v*w, v.T)
        vd = 1 / np.sqrt(vv.diagonal())
        vv = vd[:, np.newaxis] * vv * vd
        assert_allclose(vv, np.eye(N), rtol, atol)

        # check that the integral of 1 is correct
        assert_allclose(w.sum(), mu, rtol, atol)

        # compare the results of integrating a function with quad.
        f = lambda x: x**3 - 3*x**2 + x - 2
        resI = integrate.quad(lambda x: f(x)*weight_func(x), a, b)
        resG = np.vdot(f(x), w)
        rtol = 1e-6 if 1e-6 < resI[1] else resI[1] * 10
        assert_allclose(resI[0], resG, rtol=rtol) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def test_roots_hermitenorm():
    rootf = sc.roots_hermitenorm
    evalf = orth.eval_hermitenorm
    weightf = orth.hermitenorm(5).weight_func

    verify_gauss_quad(rootf, evalf, weightf, -np.inf, np.inf, 5)
    verify_gauss_quad(rootf, evalf, weightf, -np.inf, np.inf, 25, atol=1e-13)
    verify_gauss_quad(rootf, evalf, weightf, -np.inf, np.inf, 100, atol=1e-12)

    x, w = sc.roots_hermitenorm(5, False)
    y, v, m = sc.roots_hermitenorm(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -np.inf, np.inf)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_hermitenorm, 0)
    assert_raises(ValueError, sc.roots_hermitenorm, 3.3) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def test_roots_chebyt():
    weightf = orth.chebyt(5).weight_func
    verify_gauss_quad(sc.roots_chebyt, orth.eval_chebyt, weightf, -1., 1., 5)
    verify_gauss_quad(sc.roots_chebyt, orth.eval_chebyt, weightf, -1., 1., 25)
    verify_gauss_quad(sc.roots_chebyt, orth.eval_chebyt, weightf, -1., 1., 100, atol=1e-12)

    x, w = sc.roots_chebyt(5, False)
    y, v, m = sc.roots_chebyt(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -1, 1)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_chebyt, 0)
    assert_raises(ValueError, sc.roots_chebyt, 3.3) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def test_roots_chebyu():
    weightf = orth.chebyu(5).weight_func
    verify_gauss_quad(sc.roots_chebyu, orth.eval_chebyu, weightf, -1., 1., 5)
    verify_gauss_quad(sc.roots_chebyu, orth.eval_chebyu, weightf, -1., 1., 25)
    verify_gauss_quad(sc.roots_chebyu, orth.eval_chebyu, weightf, -1., 1., 100)

    x, w = sc.roots_chebyu(5, False)
    y, v, m = sc.roots_chebyu(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -1, 1)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_chebyu, 0)
    assert_raises(ValueError, sc.roots_chebyu, 3.3) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def test_roots_chebyc():
    weightf = orth.chebyc(5).weight_func
    verify_gauss_quad(sc.roots_chebyc, orth.eval_chebyc, weightf, -2., 2., 5)
    verify_gauss_quad(sc.roots_chebyc, orth.eval_chebyc, weightf, -2., 2., 25)
    verify_gauss_quad(sc.roots_chebyc, orth.eval_chebyc, weightf, -2., 2., 100, atol=1e-12)

    x, w = sc.roots_chebyc(5, False)
    y, v, m = sc.roots_chebyc(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -2, 2)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_chebyc, 0)
    assert_raises(ValueError, sc.roots_chebyc, 3.3) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def test_roots_chebys():
    weightf = orth.chebys(5).weight_func
    verify_gauss_quad(sc.roots_chebys, orth.eval_chebys, weightf, -2., 2., 5)
    verify_gauss_quad(sc.roots_chebys, orth.eval_chebys, weightf, -2., 2., 25)
    verify_gauss_quad(sc.roots_chebys, orth.eval_chebys, weightf, -2., 2., 100)

    x, w = sc.roots_chebys(5, False)
    y, v, m = sc.roots_chebys(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -2, 2)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_chebys, 0)
    assert_raises(ValueError, sc.roots_chebys, 3.3) 

# 需要导入模块: from scipy import integrate [as 别名]
# 或者: from scipy.integrate import quad [as 别名]
def vib_ent(temp, t_d_perov, t_d_brownm):
    """
    Vibrational entropy based on the Debye model
    :param temp:        temperature
    :param delta:       non-stoichiometry delta
    :return:            vibrational entropy
    """

    # integral for vibrational entropy using the Debye model
    def s_int(temp, t_d):
        def d_y(temp, t_d):
            y = t_d / temp

            def integrand(x):
                return x ** 3 / (pd.np.exp(x) - 1)

            if temp != 0:
                integral_y = quad(integrand, 0, y)[0]
                d = integral_y * (3 / (y ** 3))
            else:
                d = 0
            return d

        y = t_d / temp
        s = R * (-3 * pd.np.log(1 - pd.np.exp(-y)) + 4 * d_y(temp, t_d))

        return s

    s_perov = s_int(temp, t_d_perov)
    s_brownm = s_int(temp, t_d_brownm)

    return 2 * s_perov - (2 * s_brownm)