我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用sympy.integrate()。
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`. pyra = numpy.array([ [-1, -1, -1], [+1, -1, -1], [+1, +1, -1], [-1, +1, -1], [0, 0, 1], ]) degree = check_degree( lambda poly: quadpy.pyramid.integrate(poly, pyra, scheme), lambda k: _integrate_exact(k, pyra), 3, scheme.degree + 1 ) assert degree == scheme.degree return
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`. tetrahedron = 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 = check_degree( lambda poly: quadpy.tetrahedron.integrate( poly, tetrahedron, scheme ), integrate_monomial_over_unit_simplex, 3, scheme.degree + 1, ) assert degree == scheme.degree, \ 'Observed: {}, expected: {}'.format(degree, scheme.degree) return
def _integrate_exact(f, quadrilateral): xi = sympy.DeferredVector('xi') pxi = quadrilateral[0] * 0.25*(1.0 + xi[0])*(1.0 + xi[1]) \ + quadrilateral[1] * 0.25*(1.0 - xi[0])*(1.0 + xi[1]) \ + quadrilateral[2] * 0.25*(1.0 - xi[0])*(1.0 - xi[1]) \ + quadrilateral[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)
def __init__(self, index): self.points = numpy.linspace(-1.0, 1.0, index+1) self.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 self.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)) \ / (math.factorial(r) * math.factorial(n-r)) \ / index self.weights[r] = alpha return
def test_units(): assert convert_to((5*m/s * day) / km, 1) == 432 assert convert_to(foot / meter, meter) == Rational('0.3048') # amu is a pure mass so mass/mass gives a number, not an amount (mol) # TODO: need better simplification routine: assert str(convert_to(grams/amu, grams).n(2)) == '6.0e+23' # Light from the sun needs about 8.3 minutes to reach earth t = (1*au / speed_of_light) / minute # TODO: need a better way to simplify expressions containing units: t = convert_to(convert_to(t, meter / minute), meter) assert t == 49865956897/5995849160 # TODO: fix this, it should give `m` without `Abs` assert sqrt(m**2) == Abs(m) assert (sqrt(m))**2 == m t = Symbol('t') assert integrate(t*m/s, (t, 1*s, 5*s)) == 12*m*s assert (t * m/s).integrate((t, 1*s, 5*s)) == 12*m*s
def compute_moments(w, a, b, n, polynomial_class=lambda k, x: x**k): '''Symbolically calculate the first n moments int_a^b w(x) P_k(x) dx where `P_k` is the `k`th polynomials of a specified class. The default settings are monomials, i.e., `P_k(x)=x^k`, but you can provide any function with the signature `p(k, x)`, e.g., `sympy.polys.orthopolys.legendre_poly` scaled by the inverse of its leading coefficient `(2n)! / 2^n / (n!)^2`. ''' x = sympy.Symbol('x') return numpy.array([ sympy.integrate(w(x) * polynomial_class(k, x), (x, a, b)) for k in range(n) ])
def test_integral0(n=4): polar = sympy.Symbol('theta', real=True) azimuthal = sympy.Symbol('phi', real=True) tree = numpy.concatenate( orthopy.sphere.sph_tree( n, polar, azimuthal, normalization='quantum mechanic', symbolic=True )) assert sympy.integrate( tree[0] * sympy.sin(polar), (polar, 0, pi), (azimuthal, 0, 2*pi) ) == 2*sqrt(pi) for val in tree[1:]: assert sympy.integrate( val * sympy.sin(polar), (azimuthal, 0, 2*pi), (polar, 0, pi) ) == 0 return
def test_normality(n=3): '''Make sure that the polynomials are normal. ''' polar = sympy.Symbol('theta', real=True) azimuthal = sympy.Symbol('phi', real=True) tree = numpy.concatenate( orthopy.sphere.sph_tree( n, polar, azimuthal, normalization='quantum mechanic', symbolic=True )) for val in tree: integrand = sympy.simplify( val * sympy.conjugate(val) * sympy.sin(polar) ) assert sympy.integrate( integrand, (azimuthal, 0, 2*pi), (polar, 0, pi) ) == 1 return
def test_orthogonality(normalization, n=4): polar = sympy.Symbol('theta', real=True) azimuthal = sympy.Symbol('phi', real=True) tree = numpy.concatenate( orthopy.sphere.sph_tree( n, polar, azimuthal, normalization=normalization, symbolic=True )) vals = tree * sympy.conjugate(numpy.roll(tree, 1, axis=0)) for val in vals: integrand = sympy.simplify(val * sympy.sin(polar)) assert sympy.integrate( integrand, (azimuthal, 0, 2*pi), (polar, 0, pi) ) == 0 return
def test_integral0(n=4): '''Make sure that the polynomials are orthonormal ''' x = sympy.Symbol('x') y = sympy.Symbol('y') vals = numpy.concatenate( orthopy.e2r2.tree(n, numpy.array([x, y]), symbolic=True) ) assert sympy.integrate( vals[0] * sympy.exp(-x**2-y**2), (x, -oo, +oo), (y, -oo, +oo) ) == sympy.sqrt(sympy.pi) for val in vals[1:]: assert sympy.integrate( val * sympy.exp(-x**2-y**2), (x, -oo, +oo), (y, -oo, +oo) ) == 0 return
def green(f, curveXY): f = -sp.integrate(sp.sympify(f), y) f = f.subs({x:curveXY[0], y:curveXY[1]}) f = sp.integrate(f * sp.diff(curveXY[0], t), (t, 0, 1)) return f
def integrate(self, *, equation : str): ''' Integrate an equation with respect to x (dx) ''' x = sympy.symbols('x') try: await self.bot.embed_reply("`{}`".format(sympy.integrate(equation.strip('`'), x)), title = "Integral of {}".format(equation)) except Exception as e: await self.bot.embed_reply(py_code_block.format("{}: {}".format(type(e).__name__, e)), title = "Error")
def integrate_definite(self, lower_limit : str, upper_limit : str, *, equation : str): ''' Definite integral of an equation with respect to x (dx) ''' x = sympy.symbols('x') try: await self.bot.embed_reply("`{}`".format(sympy.integrate(equation.strip('`'), (x, lower_limit, upper_limit))), title = "Definite Integral of {} from {} to {}".format(equation, lower_limit, upper_limit)) except Exception as e: await self.bot.embed_reply(py_code_block.format("{}: {}".format(type(e).__name__, e)), title = "Error")
def sol(self): # Do the inner product! - we are in cyl coordinates! j, Sig = self.fcts() jTSj = j.T*Sig*j # we are integrating in cyl coordinates ans = sympy.integrate(sympy.integrate(sympy.integrate(r * jTSj, (r, 0, 1)), (t, 0, 2*sympy.pi)), (z, 0, 1))[0] # The `[0]` is to make it an int. return ans
def sol(self): h, Sig = self.fcts() # Do the inner product! - we are in cyl coordinates! hTSh = h.T*Sig*h ans = sympy.integrate(sympy.integrate(sympy.integrate(r * hTSh, (r, 0, 1)), (t, 0, 2*sympy.pi)), (z, 0, 1))[0] # The `[0]` is to make it an int. return ans
def _integrate_exact(f, tetrahedron): # # 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 tetrahedron [(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 tetrahedron. (See, e.g., # <http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF>). # xi = sympy.DeferredVector('xi') x_xi = \ + tetrahedron[0] * (1 - xi[0] - xi[1] - xi[2]) \ + tetrahedron[1] * xi[0] \ + tetrahedron[2] * xi[1] \ + tetrahedron[3] * xi[2] abs_det_J = 6 * quadpy.tetrahedron.volume(tetrahedron) 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)
def test_scheme(scheme, tol): # Test integration until we get to a polynomial degree `d` that can no # longer be integrated exactly. The scheme's degree is `d-1`. def eval_orthopolys(x): return numpy.concatenate( orthopy.quadrilateral.tree(scheme.degree+1, x, symbolic=False) ) quad = quadpy.quadrilateral.rectangle_points([-1.0, +1.0], [-1.0, +1.0]) vals = quadpy.quadrilateral.integrate(eval_orthopolys, quad, scheme) # 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 = check_degree_ortho(approximate, exact, abs_tol=tol) assert degree >= scheme.degree, \ 'Observed: {}, expected: {}'.format(degree, scheme.degree) return
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.triangle.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)
def test_scheme(scheme, tol): triangle = numpy.array([ [0.0, 0.0], [1.0, 0.0], [0.0, 1.0] ]) def eval_orthopolys(x): bary = numpy.array([x[0], x[1], 1.0-x[0]-x[1]]) out = numpy.concatenate(orthopy.triangle.tree( scheme.degree+1, bary, 'normal', symbolic=False )) return out vals = quadpy.triangle.integrate(eval_orthopolys, triangle, scheme) # 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] = numpy.sqrt(2.0) / 2 degree = check_degree_ortho(approximate, exact, abs_tol=tol) assert degree >= scheme.degree, \ 'Observed: {}, expected: {}'.format(degree, scheme.degree) return
def test_norm(n=1): # Maximum "n" which is tested: for i in range(n + 1): assert integrate( wavefunction(i, x) * wavefunction(-i, x), (x, 0, 2 * pi)) == 1
def test_orthogonality(n=1): # Maximum "n" which is tested: for i in range(n + 1): for j in range(i+1, n+1): assert integrate( wavefunction(i, x) * wavefunction(j, x), (x, 0, 2 * pi)) == 0
def test_norm(): # Maximum "n" which is tested: n_max = 2 # you can test any n and it works, but it's slow, so it's commented out: #n_max = 4 for n in range(n_max + 1): for l in range(n): assert integrate(R_nl(n, l, r)**2 * r**2, (r, 0, oo)) == 1
def test_norm(n=1): # Maximum "n" which is tested: for i in range(n + 1): assert integrate(psi_n(i, x, 1, 1)**2, (x, -oo, oo)) == 1
def test_orthogonality(n=1): # Maximum "n" which is tested: for i in range(n + 1): for j in range(i + 1, n + 1): assert integrate( psi_n(i, x, 1, 1)*psi_n(j, x, 1, 1), (x, -oo, oo)) == 0
def test_issue_5981(): u = symbols('u') assert integrate(1/(u**2 + 1)) == atan(u)
def test_issue_4023(): sage.var("a x") log = sage.log i = sympy.integrate(log(x)/a, (x, a, a + 1)) i2 = sympy.simplify(i) s = sage.SR(i2) assert s == (a*log(1 + a) - a*log(a) + log(1 + a) - 1)/a # This string contains Sage doctests, that execute all the functions above. # When you add a new function, please add it here as well.
def test_adaptive_limits(): """ensures that the n value given is sufficient for the curve function. This function has calculated n of 366 on the interval 0 to 2""" t = symbols('t') curve_func = t**2 + t + 1 integ = integrate(curve_func, t) value = lambdify([t], integ) value = value(2) print value apt = abs(adaptive_trapzint(curve, 0, 2)[0]-value)<1E-5 msg = "N is too small." assert apt, msg
def bench_integrate(): x, y = symbols('x y') f = (1 / tan(x)) ** 10 return integrate(f, x)
def test_norm(): # Maximum "n" which is tested: n_max = 2 # it works, but is slow, for n_max > 2 for n in range(n_max + 1): for l in range(n): assert integrate(R_nl(n, l, r)**2 * r**2, (r, 0, oo)) == 1
def test_issue_10798(): from sympy import integrate, pi, I, log, polylog, exp_polar, Piecewise, meijerg, Abs from sympy.abc import x, y assert integrate(1/(1-(x*y)**2), (x, 0, 1), y) == \ -Piecewise((I*pi*log(y) - polylog(2, y), Abs(y) < 1), (-I*pi*log(1/y) - polylog(2, y), Abs(1/y) < 1), \ (-I*pi*meijerg(((), (1, 1)), ((0, 0), ()), y) + I*pi*meijerg(((1, 1), ()), ((), (0, 0)), y) - polylog(2, y), True))/2 \ - log(y)*log(1 - 1/y)/2 + log(y)*log(1 + 1/y)/2 + log(y)*log(y - 1)/2 \ - log(y)*log(y + 1)/2 + I*pi*log(y)/2 - polylog(2, y*exp_polar(I*pi))/2
def test_issue_10488(): a,b,c,x = symbols('a b c x', real=True, positive=True) assert integrate(x/(a*x+b),x) == x/a - b*log(a*x + b)/a**2
def test_issue_4023(): sage.var("a x") log = sage.log i = sympy.integrate(log(x)/a, (x, a, a + 1)) i2 = sympy.simplify(i) s = sage.SR(i2) assert s == (a*log(1 + a) - a*log(a) + log(1 + a) - 1)/a
def stieltjes(w, a, b, n): 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)) alpha[k] = sympy.integrate(t * pi[k]**2 * w(t), (t, a, b)) / mu[k] beta[k] = mu[0] # not used, by convection 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)) alpha[k] = sympy.integrate(t * pi[k]**2 * w(t), (t, a, b)) / 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)) alpha[k] = sympy.integrate(t * pi[k]**2 * w(t), (t, a, b)) / mu[k] beta[k] = mu[k] / mu[k-1] return alpha, beta
def test_integral0(n=4): b0 = sympy.Symbol('b0') b1 = sympy.Symbol('b1') b = numpy.array([b0, b1, 1-b0-b1]) tree = numpy.concatenate( orthopy.triangle.tree(n, b, 'normal', symbolic=True) ) assert \ sympy.integrate(tree[0], (b0, 0, 1-b1), (b1, 0, 1)) == sympy.sqrt(2)/2 for val in tree[1:]: assert sympy.integrate(val, (b0, 0, 1-b1), (b1, 0, 1)) == 0 return
def test_normality(n=4): '''Make sure that the polynomials are orthonormal ''' b0 = sympy.Symbol('b0') b1 = sympy.Symbol('b1') b = numpy.array([b0, b1, 1-b0-b1]) tree = numpy.concatenate( orthopy.triangle.tree(n, b, 'normal', symbolic=True) ) for val in tree: assert sympy.integrate(val**2, (b0, 0, 1-b1), (b1, 0, 1)) == 1 return
def test_orthogonality(n=4): b0 = sympy.Symbol('b0') b1 = sympy.Symbol('b1') b = numpy.array([b0, b1, 1-b0-b1]) tree = numpy.concatenate( orthopy.triangle.tree(n, b, 'normal', symbolic=True) ) shifts = tree * numpy.roll(tree, 1, axis=0) for val in shifts: assert sympy.integrate(val, (b0, 0, 1-b1), (b1, 0, 1)) == 0 return
def test_orthogonality(n=4): x = sympy.Symbol('x') y = sympy.Symbol('y') z = sympy.Symbol('z') tree = numpy.concatenate( orthopy.enr2.tree(n, numpy.array([x, y, z]), symbolic=True) ) vals = tree * numpy.roll(tree, 1, axis=0) for val in vals: assert sympy.integrate( val * sympy.exp(-x**2 - y**2 - z**2), (x, -oo, +oo), (y, -oo, +oo), (z, -oo, +oo) ) == 0 return
def test_normality(n=4): x = sympy.Symbol('x') y = sympy.Symbol('y') z = sympy.Symbol('z') tree = numpy.concatenate( orthopy.enr2.tree(n, numpy.array([x, y, z]), symbolic=True) ) for val in tree: assert sympy.integrate( val**2 * sympy.exp(-x**2 - y**2 - z**2), (x, -oo, +oo), (y, -oo, +oo), (z, -oo, +oo) ) == 1 return
def test_integral0(n=4): '''Make sure that the polynomials are orthonormal ''' x = sympy.Symbol('x') y = sympy.Symbol('y') vals = numpy.concatenate( orthopy.quadrilateral.tree(n, numpy.array([x, y]), symbolic=True) ) assert sympy.integrate(vals[0], (x, -1, +1), (y, -1, +1)) == 2 for val in vals[1:]: assert sympy.integrate(val, (x, -1, +1), (y, -1, +1)) == 0 return
def test_orthogonality(n=4): x = sympy.Symbol('x') y = sympy.Symbol('y') tree = numpy.concatenate( orthopy.quadrilateral.tree(n, numpy.array([x, y]), symbolic=True) ) vals = tree * numpy.roll(tree, 1, axis=0) for val in vals: assert sympy.integrate(val, (x, -1, +1), (y, -1, +1)) == 0 return
def test_normality(n=4): x = sympy.Symbol('x') y = sympy.Symbol('y') tree = numpy.concatenate( orthopy.quadrilateral.tree(n, numpy.array([x, y]), symbolic=True) ) for val in tree: assert sympy.integrate(val**2, (x, -1, +1), (y, -1, +1)) == 1 return
def test_orthogonality(n=4): x = sympy.Symbol('x') tree = numpy.concatenate( orthopy.e1r2.tree(n, numpy.array([x]), symbolic=True) ) vals = tree * numpy.roll(tree, 1, axis=0) for val in vals: assert sympy.integrate(val * sympy.exp(-x**2), (x, -oo, +oo)) == 0 return
def test_normality(n=4): x = sympy.Symbol('x') tree = numpy.concatenate( orthopy.e1r2.tree(n, numpy.array([x]), symbolic=True) ) for val in tree: assert sympy.integrate( val**2 * sympy.exp(-x**2), (x, -oo, +oo) ) == 1 return
def test_integral0(n=4, dim=5): '''Make sure that the polynomials are orthonormal ''' variables = numpy.array([ sympy.Symbol('x{}'.format(k)) for k in range(dim) ]) vals = numpy.concatenate(orthopy.ncube.tree(n, variables, symbolic=True)) integration_limits = [(variable, -1, +1) for variable in variables] assert sympy.integrate(vals[0], *integration_limits) == sympy.sqrt(2)**dim for val in vals[1:]: assert sympy.integrate(val, *integration_limits) == 0 return
def test_orthogonality(n=4, dim=5): variables = numpy.array([ sympy.Symbol('x{}'.format(k)) for k in range(dim) ]) tree = numpy.concatenate(orthopy.ncube.tree(n, variables, symbolic=True)) vals = tree * numpy.roll(tree, 1, axis=0) integration_limits = [(variable, -1, +1) for variable in variables] for val in vals: assert sympy.integrate(val, *integration_limits) == 0 return
def test_normality(n=4, dim=5): variables = numpy.array([ sympy.Symbol('x{}'.format(k)) for k in range(dim) ]) tree = numpy.concatenate(orthopy.ncube.tree(n, variables, symbolic=True)) integration_limits = [(variable, -1, +1) for variable in variables] for val in tree: assert sympy.integrate(val**2, *integration_limits) == 1 return
