Python sympy.integrate方法代码示例

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def integrate(f, a, b, **kwargs):
    """Symbolically calculate the integrals

      int_a^b f_k(x) dx.

    Useful for computing the moments `w(x) * P_k(x)`, e.g.,

    moments = quadpy.tools.integrate(
            lambda x: [x**k for k in range(5)],
            -1, +1
            )

    Any keyword arguments are passed directly to `sympy.integrate` function.
    """
    x = sympy.Symbol("x")
    return numpy.array([sympy.integrate(fun, (x, a, b), **kwargs) for fun in f(x)]) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def newton_cotes_open(index, **kwargs):
    """
    Open Newton-Cotes formulae.
    <https://math.stackexchange.com/a/1959071/36678>
    """
    points = numpy.linspace(-1.0, 1.0, index + 2)[1:-1]
    degree = index if (index + 1) % 2 == 0 else index - 1
    #
    n = index + 1
    weights = numpy.empty(n - 1)
    t = sympy.Symbol("t")
    for r in range(1, n):
        # Compare with get_weights().
        f = sympy.prod([(t - i) for i in range(1, n) if i != r])
        alpha = (
            2
            * (-1) ** (n - r + 1)
            * sympy.integrate(f, (t, 0, n), **kwargs)
            / (math.factorial(r - 1) * math.factorial(n - 1 - r))
            / n
        )
        weights[r - 1] = alpha
    return C1Scheme("Newton-Cotes (open)", degree, weights, points) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def _integrate_exact(f, triangle):
    #
    # Note that
    #
    #     \int_T f(x) dx = \int_T0 |J(xi)| f(P(xi)) dxi
    #
    # with
    #
    #     P(xi) = x0 * (1-xi[0]-xi[1]) + x1 * xi[0] + x2 * xi[1].
    #
    # and T0 being the reference triangle [(0.0, 0.0), (1.0, 0.0), (0.0,
    # 1.0)].
    # The determinant of the transformation matrix J equals twice the volume of
    # the triangle. (See, e.g.,
    # <http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF>).
    #
    xi = sympy.DeferredVector("xi")
    x_xi = (
        +triangle[0] * (1 - xi[0] - xi[1]) + triangle[1] * xi[0] + triangle[2] * xi[1]
    )
    abs_det_J = 2 * quadpy.t2.volume(triangle)
    exact = sympy.integrate(
        sympy.integrate(abs_det_J * f(x_xi), (xi[1], 0, 1 - xi[0])), (xi[0], 0, 1)
    )
    return float(exact) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def test_scheme(scheme):
    assert scheme.points.dtype in [numpy.float64, numpy.int64], scheme.name
    assert scheme.weights.dtype in [numpy.float64, numpy.int64], scheme.name

    print(scheme)

    # Test integration until we get to a polynomial degree `d` that can no longer be
    # integrated exactly. The scheme's degree is `d-1`.
    pyra = numpy.array(
        [[-1, -1, -1], [+1, -1, -1], [+1, +1, -1], [-1, +1, -1], [0, 0, 1]]
    )
    degree, err = check_degree(
        lambda poly: scheme.integrate(poly, pyra),
        lambda k: _integrate_exact(k, pyra),
        3,
        scheme.degree + 1,
        tol=scheme.test_tolerance,
    )
    assert (
        degree >= scheme.degree
    ), "{} -- Observed: {}, expected: {} (max err: {:.3e})".format(
        scheme.name, degree, scheme.degree, err
    ) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def _integrate_exact(f, c2):
    xi = sympy.DeferredVector("xi")
    pxi = (
        c2[0] * 0.25 * (1.0 + xi[0]) * (1.0 + xi[1])
        + c2[1] * 0.25 * (1.0 - xi[0]) * (1.0 + xi[1])
        + c2[2] * 0.25 * (1.0 - xi[0]) * (1.0 - xi[1])
        + c2[3] * 0.25 * (1.0 + xi[0]) * (1.0 - xi[1])
    )
    pxi = [sympy.expand(pxi[0]), sympy.expand(pxi[1])]
    # determinant of the transformation matrix
    det_J = +sympy.diff(pxi[0], xi[0]) * sympy.diff(pxi[1], xi[1]) - sympy.diff(
        pxi[1], xi[0]
    ) * sympy.diff(pxi[0], xi[1])
    # we cannot use abs(), see <https://github.com/sympy/sympy/issues/4212>.
    abs_det_J = sympy.Piecewise((det_J, det_J >= 0), (-det_J, det_J < 0))

    g_xi = f(pxi)

    exact = sympy.integrate(
        sympy.integrate(abs_det_J * g_xi, (xi[1], -1, 1)), (xi[0], -1, 1)
    )
    return float(exact) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def test_scheme(scheme):
    assert scheme.points.dtype in [numpy.float64, numpy.int64], scheme.name
    assert scheme.weights.dtype in [numpy.float64, numpy.int64], scheme.name

    print(scheme)

    # Test integration until we get to a polynomial degree `d` that can no
    # longer be integrated exactly. The scheme's degree is `d-1`.
    t3 = numpy.array(
        [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
    )

    degree, err = check_degree(
        lambda poly: scheme.integrate(poly, t3),
        integrate_monomial_over_unit_simplex,
        3,
        scheme.degree + 1,
        scheme.test_tolerance,
    )

    assert (
        degree >= scheme.degree
    ), "{} -- observed: {}, expected: {} (max err: {:.3e})".format(
        scheme.name, degree, scheme.degree, err
    ) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def test_gaussian():
    """
    Make sure that symfit.distributions.Gaussians produces the expected
    sympy expression.
    """
    x0 = Parameter()
    sig = Parameter(positive=True)
    x = Variable()

    new = sympy.exp(-(x - x0)**2/(2*sig**2))/sympy.sqrt((2*sympy.pi*sig**2))
    assert isinstance(new, sympy.Expr)
    g = Gaussian(x, x0, sig)
    assert issubclass(g.__class__, sympy.Expr)
    assert new == g

    # A pdf should always integrate to 1 on its domain
    assert sympy.integrate(g, (x, -sympy.oo, sympy.oo)) == 1 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def test_exp():
    """
    Make sure that symfit.distributions.Exp produces the expected
    sympy expression.
    """
    l = Parameter(positive=True)
    x = Variable()

    new = l * sympy.exp(- l * x)
    assert isinstance(new, sympy.Expr)
    e = Exp(x, l)
    assert issubclass(e.__class__, sympy.Expr)
    assert new == e

    # A pdf should always integrate to 1 on its domain
    assert sympy.integrate(e, (x, 0, sympy.oo)) == 1 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def _pressure_approx(N):
    ζ, ζ_sl = sympy.symbols('ζ ζ_sl', real=True, positive=True)

    def coefficient(n):
        Sn = _legendre(n, ζ)
        norm_square = sympy.integrate(Sn**2, (ζ, 0, 1))
        return sympy.integrate((ζ_sl - ζ) * Sn, (ζ, 0, ζ_sl)) / norm_square

    polynomial = sum([coefficient(n) * _legendre(n, ζ) for n in range(N)])
    return sympy.lambdify((ζ, ζ_sl), sympy.simplify(polynomial)) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def quadrature_degree(self, u, h, **kwargs):
        r"""Return the quadrature degree necessary to integrate the action
        functional accurately

        Firedrake uses a very conservative algorithm for estimating the
        number of quadrature points necessary to integrate a given
        expression. By exploiting known structure of the problem, we can
        reduce the number of quadrature points while preserving accuracy.
        """
        xdegree_u, zdegree_u = u.ufl_element().degree()
        degree_h = h.ufl_element().degree()[0]
        return (3 * (xdegree_u - 1) + 2 * degree_h,
                3 * max(zdegree_u - 1, 0) + zdegree_u + 1) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def stieltjes(w, a, b, n, **kwargs):
    t = sympy.Symbol("t")

    alpha = n * [None]
    beta = n * [None]
    mu = n * [None]
    pi = n * [None]

    k = 0
    pi[k] = 1
    mu[k] = sympy.integrate(pi[k] ** 2 * w(t), (t, a, b), **kwargs)
    alpha[k] = sympy.integrate(t * pi[k] ** 2 * w(t), (t, a, b), **kwargs) / mu[k]
    beta[k] = mu[0]  # not used, by convention mu[0]

    k = 1
    pi[k] = (t - alpha[k - 1]) * pi[k - 1]
    mu[k] = sympy.integrate(pi[k] ** 2 * w(t), (t, a, b), **kwargs)
    alpha[k] = sympy.integrate(t * pi[k] ** 2 * w(t), (t, a, b), **kwargs) / mu[k]
    beta[k] = mu[k] / mu[k - 1]

    for k in range(2, n):
        pi[k] = (t - alpha[k - 1]) * pi[k - 1] - beta[k - 1] * pi[k - 2]
        mu[k] = sympy.integrate(pi[k] ** 2 * w(t), (t, a, b), **kwargs)
        alpha[k] = sympy.integrate(t * pi[k] ** 2 * w(t), (t, a, b), **kwargs) / mu[k]
        beta[k] = mu[k] / mu[k - 1]

    return alpha, beta 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def newton_cotes_closed(index, **kwargs):
    """
    Closed Newton-Cotes formulae.
    <https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas#Closed_Newton.E2.80.93Cotes_formulae>,
    <http://mathworld.wolfram.com/Newton-CotesFormulas.html>.
    """
    points = numpy.linspace(-1.0, 1.0, index + 1)
    degree = index + 1 if index % 2 == 0 else index

    # Formula (26) from
    # <http://mathworld.wolfram.com/Newton-CotesFormulas.html>.
    # Note that Sympy carries out all operations in rationals, i.e.,
    # _exactly_. Only at the end, the rational is converted into a float.
    n = index
    weights = numpy.empty(n + 1)
    t = sympy.Symbol("t")
    for r in range(n + 1):
        # Compare with get_weights().
        f = sympy.prod([(t - i) for i in range(n + 1) if i != r])
        alpha = (
            2
            * (-1) ** (n - r)
            * sympy.integrate(f, (t, 0, n), **kwargs)
            / (math.factorial(r) * math.factorial(n - r))
            / index
        )
        weights[r] = alpha
    return C1Scheme("Newton-Cotes (closed)", degree, weights, points) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def test_scheme(scheme):
    # Test integration until we get to a polynomial degree `d` that can no longer be
    # integrated exactly. The scheme's degree is `d-1`.
    assert scheme.points.dtype in [numpy.float64, numpy.int64], scheme.name
    assert scheme.weights.dtype in [numpy.float64, numpy.int64], scheme.name

    print(scheme)

    def eval_orthopolys(x):
        return numpy.concatenate(
            orthopy.quadrilateral.tree(x, scheme.degree + 1, symbolic=False)
        )

    quad = quadpy.c2.rectangle_points([-1.0, +1.0], [-1.0, +1.0])
    vals = scheme.integrate(eval_orthopolys, quad)
    # Put vals back into the tree structure:
    # len(approximate[k]) == k+1
    approximate = [
        vals[k * (k + 1) // 2 : (k + 1) * (k + 2) // 2]
        for k in range(scheme.degree + 2)
    ]

    exact = [numpy.zeros(k + 1) for k in range(scheme.degree + 2)]
    exact[0][0] = 2.0

    degree, err = check_degree_ortho(approximate, exact, abs_tol=scheme.test_tolerance)

    assert (
        degree >= scheme.degree
    ), "{} -- observed: {}, expected: {} (max err: {:.3e})".format(
        scheme.name, degree, scheme.degree, err
    ) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def _integrate_exact(f, t3):
    #
    # Note that
    #
    #     \int_T f(x) dx = \int_T0 |J(xi)| f(P(xi)) dxi
    #
    # with
    #
    #     P(xi) = x0 * (1-xi[0]-xi[1]) + x1 * xi[0] + x2 * xi[1].
    #
    # and T0 being the reference t3 [(0.0, 0.0), (1.0, 0.0), (0.0,
    # 1.0)].
    # The determinant of the transformation matrix J equals twice the volume of
    # the t3. (See, e.g.,
    # <http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF>).
    #
    xi = sympy.DeferredVector("xi")
    x_xi = (
        +t3[0] * (1 - xi[0] - xi[1] - xi[2])
        + t3[1] * xi[0]
        + t3[2] * xi[1]
        + t3[3] * xi[2]
    )
    abs_det_J = 6 * quadpy.t3.volume(t3)
    exact = sympy.integrate(
        sympy.integrate(
            sympy.integrate(abs_det_J * f(x_xi), (xi[2], 0, 1 - xi[0] - xi[1])),
            (xi[1], 0, 1 - xi[0]),
        ),
        (xi[0], 0, 1),
    )
    return float(exact) 

# 需要导入模块: import sympy [as 别名]
# 或者: from sympy import integrate [as 别名]
def phi(i, t):
    if i == 0:
        return 1  # Initial history x(t)=1 on [-1,0]
    else:
        return phi(i-1, i-1) - integrate(phi(i-1, xi-1), (xi, i-1, t))