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
python类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`.
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(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
