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)
python类expand()的实例源码
def test_use():
assert use(0, expand) == 0
f = (x + y)**2*x + 1
assert use(f, expand, level=0) == x**3 + 2*x**2*y + x*y**2 + + 1
assert use(f, expand, level=1) == x**3 + 2*x**2*y + x*y**2 + + 1
assert use(f, expand, level=2) == 1 + x*(2*x*y + x**2 + y**2)
assert use(f, expand, level=3) == (x + y)**2*x + 1
f = (x**2 + 1)**2 - 1
kwargs = {'gaussian': True}
assert use(f, factor, level=0, kwargs=kwargs) == x**2*(x**2 + 2)
assert use(f, factor, level=1, kwargs=kwargs) == (x + I)**2*(x - I)**2 - 1
assert use(f, factor, level=2, kwargs=kwargs) == (x + I)**2*(x - I)**2 - 1
assert use(f, factor, level=3, kwargs=kwargs) == (x**2 + 1)**2 - 1
def _get_constant_subexpressions(expr, Cs):
Cs = set(Cs)
Ces = []
def _recursive_walk(expr):
expr_syms = expr.free_symbols
if len(expr_syms) > 0 and expr_syms.issubset(Cs):
Ces.append(expr)
else:
if expr.func == exp:
expr = expr.expand(mul=True)
if expr.func in (Add, Mul):
d = sift(expr.args, lambda i : i.free_symbols.issubset(Cs))
if len(d[True]) > 1:
x = expr.func(*d[True])
if not x.is_number:
Ces.append(x)
elif isinstance(expr, C.Integral):
if expr.free_symbols.issubset(Cs) and \
all(map(lambda x: len(x) == 3, expr.limits)):
Ces.append(expr)
for i in expr.args:
_recursive_walk(i)
return
_recursive_walk(expr)
return Ces
def test_one_dof():
# This is for a 1 dof spring-mass-damper case.
# It is described in more detail in the KanesMethod docstring.
q, u = dynamicsymbols('q u')
qd, ud = dynamicsymbols('q u', 1)
m, c, k = symbols('m c k')
N = ReferenceFrame('N')
P = Point('P')
P.set_vel(N, u * N.x)
kd = [qd - u]
FL = [(P, (-k * q - c * u) * N.x)]
pa = Particle('pa', P, m)
BL = [pa]
KM = KanesMethod(N, [q], [u], kd)
KM.kanes_equations(FL, BL)
MM = KM.mass_matrix
forcing = KM.forcing
rhs = MM.inv() * forcing
assert expand(rhs[0]) == expand(-(q * k + u * c) / m)
assert KM.linearize() == \
(Matrix([[0, 1], [-k, -c]]), Matrix([]), Matrix([]))
def test_pend():
q, u = dynamicsymbols('q u')
qd, ud = dynamicsymbols('q u', 1)
m, l, g = symbols('m l g')
N = ReferenceFrame('N')
P = Point('P')
P.set_vel(N, -l * u * sin(q) * N.x + l * u * cos(q) * N.y)
kd = [qd - u]
FL = [(P, m * g * N.x)]
pa = Particle('pa', P, m)
BL = [pa]
KM = KanesMethod(N, [q], [u], kd)
KM.kanes_equations(FL, BL)
MM = KM.mass_matrix
forcing = KM.forcing
rhs = MM.inv() * forcing
rhs.simplify()
assert expand(rhs[0]) == expand(-g / l * sin(q))
def _guess_expansion(f, x):
""" Try to guess sensible rewritings for integrand f(x). """
from sympy import expand_trig
from sympy.functions.elementary.trigonometric import TrigonometricFunction
from sympy.functions.elementary.hyperbolic import HyperbolicFunction
res = [(f, 'originial integrand')]
expanded = expand_mul(res[-1][0])
if expanded != res[-1][0]:
res += [(expanded, 'expand_mul')]
expanded = expand(res[-1][0])
if expanded != res[-1][0]:
res += [(expanded, 'expand')]
if res[-1][0].has(TrigonometricFunction, HyperbolicFunction):
expanded = expand_mul(expand_trig(res[-1][0]))
if expanded != res[-1][0]:
res += [(expanded, 'expand_trig, expand_mul')]
return res
def test_recursive():
from sympy import symbols, exp_polar, expand
a, b, c = symbols('a b c', positive=True)
r = exp(-(x - a)**2)*exp(-(x - b)**2)
e = integrate(r, (x, 0, oo), meijerg=True)
assert simplify(e.expand()) == (
sqrt(2)*sqrt(pi)*(
(erf(sqrt(2)*(a + b)/2) + 1)*exp(-a**2/2 + a*b - b**2/2))/4)
e = integrate(exp(-(x - a)**2)*exp(-(x - b)**2)*exp(c*x), (x, 0, oo), meijerg=True)
assert simplify(e) == (
sqrt(2)*sqrt(pi)*(erf(sqrt(2)*(2*a + 2*b + c)/4) + 1)*exp(-a**2 - b**2
+ (2*a + 2*b + c)**2/8)/4)
assert simplify(integrate(exp(-(x - a - b - c)**2), (x, 0, oo), meijerg=True)) == \
sqrt(pi)/2*(1 + erf(a + b + c))
assert simplify(integrate(exp(-(x + a + b + c)**2), (x, 0, oo), meijerg=True)) == \
sqrt(pi)/2*(1 - erf(a + b + c))
def test_use():
assert use(0, expand) == 0
f = (x + y)**2*x + 1
assert use(f, expand, level=0) == x**3 + 2*x**2*y + x*y**2 + + 1
assert use(f, expand, level=1) == x**3 + 2*x**2*y + x*y**2 + + 1
assert use(f, expand, level=2) == 1 + x*(2*x*y + x**2 + y**2)
assert use(f, expand, level=3) == (x + y)**2*x + 1
f = (x**2 + 1)**2 - 1
kwargs = {'gaussian': True}
assert use(f, factor, level=0, kwargs=kwargs) == x**2*(x**2 + 2)
assert use(f, factor, level=1, kwargs=kwargs) == (x + I)**2*(x - I)**2 - 1
assert use(f, factor, level=2, kwargs=kwargs) == (x + I)**2*(x - I)**2 - 1
assert use(f, factor, level=3, kwargs=kwargs) == (x**2 + 1)**2 - 1
def test_recursive():
from sympy import symbols
a, b, c = symbols('a b c', positive=True)
r = exp(-(x - a)**2)*exp(-(x - b)**2)
e = integrate(r, (x, 0, oo), meijerg=True)
assert simplify(e.expand()) == (
sqrt(2)*sqrt(pi)*(
(erf(sqrt(2)*(a + b)/2) + 1)*exp(-a**2/2 + a*b - b**2/2))/4)
e = integrate(exp(-(x - a)**2)*exp(-(x - b)**2)*exp(c*x), (x, 0, oo), meijerg=True)
assert simplify(e) == (
sqrt(2)*sqrt(pi)*(erf(sqrt(2)*(2*a + 2*b + c)/4) + 1)*exp(-a**2 - b**2
+ (2*a + 2*b + c)**2/8)/4)
assert simplify(integrate(exp(-(x - a - b - c)**2), (x, 0, oo), meijerg=True)) == \
sqrt(pi)/2*(1 + erf(a + b + c))
assert simplify(integrate(exp(-(x + a + b + c)**2), (x, 0, oo), meijerg=True)) == \
sqrt(pi)/2*(1 - erf(a + b + c))
def find_functional_units(gpr_str):
"""
Return an iterator of gene IDs grouped by boolean rules from the gpr_str.
The gpr_str is first transformed into an algebraic expression, replacing
the boolean operators 'or' with '+' and 'and' with '*'. Treating the
gene IDs as sympy.symbols this allows a mathematical expansion of the
algebraic expression. The expanded form is then split again producing sets
of gene IDs that in the gpr_str had an 'and' relationship.
Parameters
----------
gpr_str : string
A string consisting of gene ids and the boolean expressions 'and'
and 'or'
"""
algebraic_form = re.sub('[Oo]r', '+', re.sub('[Aa]nd', '*', gpr_str))
expanded = str(expand(algebraic_form))
for unit in expanded.replace('+', ',').split(' , '):
yield unit.split('*')
def isReversible(self, formula):
formula = expand(formula)
# If we had an addition
if formula.func == SympyAdd:
# And one of the terms
for arg in formula.args:
# is *(-1)
if (arg.func == SympyMul
and (arg.args[0] == SympyInteger(-1)
or arg.args[1] == SympyInteger(-1))):
return True
return False
def can_do_meijer(a1, a2, b1, b2, numeric=True):
"""
This helper function tries to hyperexpand() the meijer g-function
corresponding to the parameters a1, a2, b1, b2.
It returns False if this expansion still contains g-functions.
If numeric is True, it also tests the so-obtained formula numerically
(at random values) and returns False if the test fails.
Else it returns True.
"""
from sympy import unpolarify, expand
r = hyperexpand(meijerg(a1, a2, b1, b2, z))
if r.has(meijerg):
return False
# NOTE hyperexpand() returns a truly branched function, whereas numerical
# evaluation only works on the main branch. Since we are evaluating on
# the main branch, this should not be a problem, but expressions like
# exp_polar(I*pi/2*x)**a are evaluated incorrectly. We thus have to get
# rid of them. The expand heuristically does this...
r = unpolarify(expand(r, force=True, power_base=True, power_exp=False,
mul=False, log=False, multinomial=False, basic=False))
if not numeric:
return True
repl = {}
for n, a in enumerate(meijerg(a1, a2, b1, b2, z).free_symbols - set([z])):
repl[a] = randcplx(n)
return tn(meijerg(a1, a2, b1, b2, z).subs(repl), r.subs(repl), z)
def __mul__(self, v):
if not isinstance(v, Vector):
self_x_v = Vector()
self_x_v.obj = self.obj * v
return self_x_v
else:
result = expand(self.obj * v.obj)
result = bilinear_product(result, Vector.dot)
return result
def doit(self, **hints):
"""Expand the density operator into an outer product format.
Examples
=========
>>> from sympy.physics.quantum.state import Ket
>>> from sympy.physics.quantum.density import Density
>>> from sympy.physics.quantum.operator import Operator
>>> A = Operator('A')
>>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
>>> d.doit()
0.5*|0><0| + 0.5*|1><1|
"""
terms = []
for (state, prob) in self.args:
state = state.expand() # needed to break up (a+b)*c
if (isinstance(state, Add)):
for arg in product(state.args, repeat=2):
terms.append(prob *
self._generate_outer_prod(arg[0], arg[1]))
else:
terms.append(prob *
self._generate_outer_prod(state, state))
return Add(*terms)
def test_two_dof():
# This is for a 2 d.o.f., 2 particle spring-mass-damper.
# The first coordinate is the displacement of the first particle, and the
# second is the relative displacement between the first and second
# particles. Speeds are defined as the time derivatives of the particles.
q1, q2, u1, u2 = dynamicsymbols('q1 q2 u1 u2')
q1d, q2d, u1d, u2d = dynamicsymbols('q1 q2 u1 u2', 1)
m, c1, c2, k1, k2 = symbols('m c1 c2 k1 k2')
N = ReferenceFrame('N')
P1 = Point('P1')
P2 = Point('P2')
P1.set_vel(N, u1 * N.x)
P2.set_vel(N, (u1 + u2) * N.x)
kd = [q1d - u1, q2d - u2]
# Now we create the list of forces, then assign properties to each
# particle, then create a list of all particles.
FL = [(P1, (-k1 * q1 - c1 * u1 + k2 * q2 + c2 * u2) * N.x), (P2, (-k2 *
q2 - c2 * u2) * N.x)]
pa1 = Particle('pa1', P1, m)
pa2 = Particle('pa2', P2, m)
BL = [pa1, pa2]
# Finally we create the KanesMethod object, specify the inertial frame,
# pass relevant information, and form Fr & Fr*. Then we calculate the mass
# matrix and forcing terms, and finally solve for the udots.
KM = KanesMethod(N, q_ind=[q1, q2], u_ind=[u1, u2], kd_eqs=kd)
KM.kanes_equations(FL, BL)
MM = KM.mass_matrix
forcing = KM.forcing
rhs = MM.inv() * forcing
assert expand(rhs[0]) == expand((-k1 * q1 - c1 * u1 + k2 * q2 + c2 * u2)/m)
assert expand(rhs[1]) == expand((k1 * q1 + c1 * u1 - 2 * k2 * q2 - 2 *
c2 * u2) / m)
def bench_expand():
x, y, z = symbols('x y z')
expand((1 + x + y + z) ** 20)
def bench_str():
x, y, z = symbols('x y z')
str(expand((x + 2 * y + 3 * z) ** 30))
def convert_to_dict(node: Node) -> dict:
children = OrderedDict()
for node_property in node.properties:
children[node_property] = convert_to_dict(node[node_property][0])
simplified = str(sympy.expand(parse_expr(''.join(to_token_sequence(node, [])))))
if len(children) > 0:
return dict(Name=node.name, Children=children, Symbol=simplified)
else:
return dict(Name=node.name, Symbol=simplified)
def non_commutative_symexpand(expr_string):
from sympy.parsing.sympy_parser import parse_expr
parsed_expr = parse_expr(expr_string, evaluate=False)
new_locals = {sym.name:sympy.Symbol(sym.name, commutative=False)
for sym in parsed_expr.atoms(sympy.Symbol)}
return sympy.expand(eval(expr_string, {}, new_locals))
def can_do_meijer(a1, a2, b1, b2, numeric=True):
"""
This helper function tries to hyperexpand() the meijer g-function
corresponding to the parameters a1, a2, b1, b2.
It returns False if this expansion still contains g-functions.
If numeric is True, it also tests the so-obtained formula numerically
(at random values) and returns False if the test fails.
Else it returns True.
"""
from sympy import unpolarify, expand
r = hyperexpand(meijerg(a1, a2, b1, b2, z))
if r.has(meijerg):
return False
# NOTE hyperexpand() returns a truly branched function, whereas numerical
# evaluation only works on the main branch. Since we are evaluating on
# the main branch, this should not be a problem, but expressions like
# exp_polar(I*pi/2*x)**a are evaluated incorrectly. We thus have to get
# rid of them. The expand heuristically does this...
r = unpolarify(expand(r, force=True, power_base=True, power_exp=False,
mul=False, log=False, multinomial=False, basic=False))
if not numeric:
return True
repl = {}
for n, a in enumerate(meijerg(a1, a2, b1, b2, z).free_symbols - {z}):
repl[a] = randcplx(n)
return tn(meijerg(a1, a2, b1, b2, z).subs(repl), r.subs(repl), z)
def doit(self, **hints):
"""Expand the density operator into an outer product format.
Examples
========
>>> from sympy.physics.quantum.state import Ket
>>> from sympy.physics.quantum.density import Density
>>> from sympy.physics.quantum.operator import Operator
>>> A = Operator('A')
>>> d = Density([Ket(0), 0.5], [Ket(1),0.5])
>>> d.doit()
0.5*|0><0| + 0.5*|1><1|
"""
terms = []
for (state, prob) in self.args:
state = state.expand() # needed to break up (a+b)*c
if (isinstance(state, Add)):
for arg in product(state.args, repeat=2):
terms.append(prob *
self._generate_outer_prod(arg[0], arg[1]))
else:
terms.append(prob *
self._generate_outer_prod(state, state))
return Add(*terms)
def _integrate_exact(k, pyra):
def f(x):
return x[0]**int(k[0]) * x[1]**int(k[1]) * x[2]**int(k[2])
# map the reference hexahedron [-1,1]^3 to the pyramid
xi = sympy.DeferredVector('xi')
pxi = (
+ pyra[0] * (1 - xi[0])*(1 - xi[1])*(1 - xi[2]) / 8
+ pyra[1] * (1 + xi[0])*(1 - xi[1])*(1 - xi[2]) / 8
+ pyra[2] * (1 + xi[0])*(1 + xi[1])*(1 - xi[2]) / 8
+ pyra[3] * (1 - xi[0])*(1 + xi[1])*(1 - xi[2]) / 8
+ pyra[4] * (1 + xi[2]) / 2
)
pxi = [
sympy.expand(pxi[0]),
sympy.expand(pxi[1]),
sympy.expand(pxi[2]),
]
# determinant of the transformation matrix
J = sympy.Matrix([
[sympy.diff(pxi[0], xi[0]),
sympy.diff(pxi[0], xi[1]),
sympy.diff(pxi[0], xi[2])],
[sympy.diff(pxi[1], xi[0]),
sympy.diff(pxi[1], xi[1]),
sympy.diff(pxi[1], xi[2])],
[sympy.diff(pxi[2], xi[0]),
sympy.diff(pxi[2], xi[1]),
sympy.diff(pxi[2], xi[2])],
])
det_J = sympy.det(J)
# 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))
# This is quite the leap of faith, but sympy will cowardly bail out
# otherwise.
abs_det_J = det_J
exact = sympy.integrate(
sympy.integrate(
sympy.integrate(abs_det_J * f(pxi), (xi[2], -1, 1)),
(xi[1], -1, +1)
),
(xi[0], -1, +1)
)
return float(exact)
def _integrate_exact(f, hexa):
xi = sympy.DeferredVector('xi')
pxi = \
+ hexa[0] * 0.125*(1.0 - xi[0])*(1.0 - xi[1])*(1.0 - xi[2]) \
+ hexa[1] * 0.125*(1.0 + xi[0])*(1.0 - xi[1])*(1.0 - xi[2]) \
+ hexa[2] * 0.125*(1.0 + xi[0])*(1.0 + xi[1])*(1.0 - xi[2]) \
+ hexa[3] * 0.125*(1.0 - xi[0])*(1.0 + xi[1])*(1.0 - xi[2]) \
+ hexa[4] * 0.125*(1.0 - xi[0])*(1.0 - xi[1])*(1.0 + xi[2]) \
+ hexa[5] * 0.125*(1.0 + xi[0])*(1.0 - xi[1])*(1.0 + xi[2]) \
+ hexa[6] * 0.125*(1.0 + xi[0])*(1.0 + xi[1])*(1.0 + xi[2]) \
+ hexa[7] * 0.125*(1.0 - xi[0])*(1.0 + xi[1])*(1.0 + xi[2])
pxi = [
sympy.expand(pxi[0]),
sympy.expand(pxi[1]),
sympy.expand(pxi[2]),
]
# determinant of the transformation matrix
J = sympy.Matrix([
[sympy.diff(pxi[0], xi[0]),
sympy.diff(pxi[0], xi[1]),
sympy.diff(pxi[0], xi[2])],
[sympy.diff(pxi[1], xi[0]),
sympy.diff(pxi[1], xi[1]),
sympy.diff(pxi[1], xi[2])],
[sympy.diff(pxi[2], xi[0]),
sympy.diff(pxi[2], xi[1]),
sympy.diff(pxi[2], xi[2])],
])
det_J = sympy.det(J)
# 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(
sympy.integrate(abs_det_J * g_xi, (xi[2], -1, 1)),
(xi[1], -1, 1)
),
(xi[0], -1, 1)
)
return float(exact)
# pylint: disable=too-many-arguments
def char_coefficients(poles):
"""
Calculate the coefficients of a characteristic polynomial.
Args:
poles (list or :obj:`numpy.ndarray`): pol configuration
Return:
:obj:`numpy.ndarray`: coefficients
"""
poles = np.array(poles) # convert to numpy array
poles = np.ravel(poles) # transform to 1d array
s = sp.symbols("s")
poly = 1
for s_i in poles:
poly = (s - s_i) * poly
poly = sp.expand(poly)
# calculate the coefficient of characteristic polynomial
n = len(poles)
p = []
for i in range(n):
p.append(poly.subs([(s, 0)]))
poly = poly - p[i]
poly = poly / s
poly = sp.expand(poly)
# convert numbers and complex objects from multiplication to a complex
# number
p = [complex(x) for x in p]
# if imaginary part is greater than the boundary, set imaginary part to zero
boundary = 1e-12
for idx, val in enumerate(p):
if abs(val.imag) > boundary:
msg = "Imaginary Part of the coefficient p[{}] "
"is unequal zero ({})) for a given tolerance of {}".format(
str(idx), str(boundary), str(val.imag))
warnings.warn(msg)
p[idx] = val.real
return np.array(p, dtype=float) # [a_0, a_1, ... , a_n-1]
def test_aux():
# Same as above, except we have 2 auxiliary speeds for the ground contact
# point, which is known to be zero. In one case, we go through then
# substitute the aux. speeds in at the end (they are zero, as well as their
# derivative), in the other case, we use the built-in auxiliary speed part
# of KanesMethod. The equations from each should be the same.
q1, q2, q3, u1, u2, u3 = dynamicsymbols('q1 q2 q3 u1 u2 u3')
q1d, q2d, q3d, u1d, u2d, u3d = dynamicsymbols('q1 q2 q3 u1 u2 u3', 1)
u4, u5, f1, f2 = dynamicsymbols('u4, u5, f1, f2')
u4d, u5d = dynamicsymbols('u4, u5', 1)
r, m, g = symbols('r m g')
N = ReferenceFrame('N')
Y = N.orientnew('Y', 'Axis', [q1, N.z])
L = Y.orientnew('L', 'Axis', [q2, Y.x])
R = L.orientnew('R', 'Axis', [q3, L.y])
w_R_N_qd = R.ang_vel_in(N)
R.set_ang_vel(N, u1 * L.x + u2 * L.y + u3 * L.z)
C = Point('C')
C.set_vel(N, u4 * L.x + u5 * (Y.z ^ L.x))
Dmc = C.locatenew('Dmc', r * L.z)
Dmc.v2pt_theory(C, N, R)
Dmc.a2pt_theory(C, N, R)
I = inertia(L, m / 4 * r**2, m / 2 * r**2, m / 4 * r**2)
kd = [dot(R.ang_vel_in(N) - w_R_N_qd, uv) for uv in L]
ForceList = [(Dmc, - m * g * Y.z), (C, f1 * L.x + f2 * (Y.z ^ L.x))]
BodyD = RigidBody('BodyD', Dmc, R, m, (I, Dmc))
BodyList = [BodyD]
KM = KanesMethod(N, q_ind=[q1, q2, q3], u_ind=[u1, u2, u3, u4, u5],
kd_eqs=kd)
(fr, frstar) = KM.kanes_equations(ForceList, BodyList)
fr = fr.subs({u4d: 0, u5d: 0}).subs({u4: 0, u5: 0})
frstar = frstar.subs({u4d: 0, u5d: 0}).subs({u4: 0, u5: 0})
KM2 = KanesMethod(N, q_ind=[q1, q2, q3], u_ind=[u1, u2, u3], kd_eqs=kd,
u_auxiliary=[u4, u5])
(fr2, frstar2) = KM2.kanes_equations(ForceList, BodyList)
fr2 = fr2.subs({u4d: 0, u5d: 0}).subs({u4: 0, u5: 0})
frstar2 = frstar2.subs({u4d: 0, u5d: 0}).subs({u4: 0, u5: 0})
frstar.simplify()
frstar2.simplify()
assert (fr - fr2).expand() == Matrix([0, 0, 0, 0, 0])
assert (frstar - frstar2).expand() == Matrix([0, 0, 0, 0, 0])
def test_inverse_laplace_transform():
from sympy import (expand, sinh, cosh, besselj, besseli, exp_polar,
unpolarify, simplify, factor_terms)
ILT = inverse_laplace_transform
a, b, c, = symbols('a b c', positive=True)
t = symbols('t')
def simp_hyp(expr):
return factor_terms(expand_mul(expr)).rewrite(sin)
# just test inverses of all of the above
assert ILT(1/s, s, t) == Heaviside(t)
assert ILT(1/s**2, s, t) == t*Heaviside(t)
assert ILT(1/s**5, s, t) == t**4*Heaviside(t)/24
assert ILT(exp(-a*s)/s, s, t) == Heaviside(t - a)
assert ILT(exp(-a*s)/(s + b), s, t) == exp(b*(a - t))*Heaviside(-a + t)
assert ILT(a/(s**2 + a**2), s, t) == sin(a*t)*Heaviside(t)
assert ILT(s/(s**2 + a**2), s, t) == cos(a*t)*Heaviside(t)
# TODO is there a way around simp_hyp?
assert simp_hyp(ILT(a/(s**2 - a**2), s, t)) == sinh(a*t)*Heaviside(t)
assert simp_hyp(ILT(s/(s**2 - a**2), s, t)) == cosh(a*t)*Heaviside(t)
assert ILT(a/((s + b)**2 + a**2), s, t) == exp(-b*t)*sin(a*t)*Heaviside(t)
assert ILT(
(s + b)/((s + b)**2 + a**2), s, t) == exp(-b*t)*cos(a*t)*Heaviside(t)
# TODO sinh/cosh shifted come out a mess. also delayed trig is a mess
# TODO should this simplify further?
assert ILT(exp(-a*s)/s**b, s, t) == \
(t - a)**(b - 1)*Heaviside(t - a)/gamma(b)
assert ILT(exp(-a*s)/sqrt(1 + s**2), s, t) == \
Heaviside(t - a)*besselj(0, a - t) # note: besselj(0, x) is even
# XXX ILT turns these branch factor into trig functions ...
assert simplify(ILT(a**b*(s + sqrt(s**2 - a**2))**(-b)/sqrt(s**2 - a**2),
s, t).rewrite(exp)) == \
Heaviside(t)*besseli(b, a*t)
assert ILT(a**b*(s + sqrt(s**2 + a**2))**(-b)/sqrt(s**2 + a**2),
s, t).rewrite(exp) == \
Heaviside(t)*besselj(b, a*t)
assert ILT(1/(s*sqrt(s + 1)), s, t) == Heaviside(t)*erf(sqrt(t))
# TODO can we make erf(t) work?
assert ILT(1/(s**2*(s**2 + 1)),s,t) == (t - sin(t))*Heaviside(t)
def test_fourier_transform():
from sympy import simplify, expand, expand_complex, factor, expand_trig
FT = fourier_transform
IFT = inverse_fourier_transform
def simp(x):
return simplify(expand_trig(expand_complex(expand(x))))
def sinc(x):
return sin(pi*x)/(pi*x)
k = symbols('k', real=True)
f = Function("f")
# TODO for this to work with real a, need to expand abs(a*x) to abs(a)*abs(x)
a = symbols('a', positive=True)
b = symbols('b', positive=True)
posk = symbols('posk', positive=True)
# Test unevaluated form
assert fourier_transform(f(x), x, k) == FourierTransform(f(x), x, k)
assert inverse_fourier_transform(
f(k), k, x) == InverseFourierTransform(f(k), k, x)
# basic examples from wikipedia
assert simp(FT(Heaviside(1 - abs(2*a*x)), x, k)) == sinc(k/a)/a
# TODO IFT is a *mess*
assert simp(FT(Heaviside(1 - abs(a*x))*(1 - abs(a*x)), x, k)) == sinc(k/a)**2/a
# TODO IFT
assert factor(FT(exp(-a*x)*Heaviside(x), x, k), extension=I) == \
1/(a + 2*pi*I*k)
# NOTE: the ift comes out in pieces
assert IFT(1/(a + 2*pi*I*x), x, posk,
noconds=False) == (exp(-a*posk), True)
assert IFT(1/(a + 2*pi*I*x), x, -posk,
noconds=False) == (0, True)
assert IFT(1/(a + 2*pi*I*x), x, symbols('k', negative=True),
noconds=False) == (0, True)
# TODO IFT without factoring comes out as meijer g
assert factor(FT(x*exp(-a*x)*Heaviside(x), x, k), extension=I) == \
1/(a + 2*pi*I*k)**2
assert FT(exp(-a*x)*sin(b*x)*Heaviside(x), x, k) == \
b/(b**2 + (a + 2*I*pi*k)**2)
assert FT(exp(-a*x**2), x, k) == sqrt(pi)*exp(-pi**2*k**2/a)/sqrt(a)
assert IFT(sqrt(pi/a)*exp(-(pi*k)**2/a), k, x) == exp(-a*x**2)
assert FT(exp(-a*abs(x)), x, k) == 2*a/(a**2 + 4*pi**2*k**2)
# TODO IFT (comes out as meijer G)
# TODO besselj(n, x), n an integer > 0 actually can be done...
# TODO are there other common transforms (no distributions!)?
def test_expint():
""" Test various exponential integrals. """
from sympy import (expint, unpolarify, Symbol, Ci, Si, Shi, Chi,
sin, cos, sinh, cosh, Ei)
assert simplify(unpolarify(integrate(exp(-z*x)/x**y, (x, 1, oo),
meijerg=True, conds='none'
).rewrite(expint).expand(func=True))) == expint(y, z)
assert integrate(exp(-z*x)/x, (x, 1, oo), meijerg=True,
conds='none').rewrite(expint).expand() == \
expint(1, z)
assert integrate(exp(-z*x)/x**2, (x, 1, oo), meijerg=True,
conds='none').rewrite(expint).expand() == \
expint(2, z).rewrite(Ei).rewrite(expint)
assert integrate(exp(-z*x)/x**3, (x, 1, oo), meijerg=True,
conds='none').rewrite(expint).expand() == \
expint(3, z).rewrite(Ei).rewrite(expint).expand()
t = Symbol('t', positive=True)
assert integrate(-cos(x)/x, (x, t, oo), meijerg=True).expand() == Ci(t)
assert integrate(-sin(x)/x, (x, t, oo), meijerg=True).expand() == \
Si(t) - pi/2
assert integrate(sin(x)/x, (x, 0, z), meijerg=True) == Si(z)
assert integrate(sinh(x)/x, (x, 0, z), meijerg=True) == Shi(z)
assert integrate(exp(-x)/x, x, meijerg=True).expand().rewrite(expint) == \
I*pi - expint(1, x)
assert integrate(exp(-x)/x**2, x, meijerg=True).rewrite(expint).expand() \
== expint(1, x) - exp(-x)/x - I*pi
u = Symbol('u', polar=True)
assert integrate(cos(u)/u, u, meijerg=True).expand().as_independent(u)[1] \
== Ci(u)
assert integrate(cosh(u)/u, u, meijerg=True).expand().as_independent(u)[1] \
== Chi(u)
assert integrate(expint(1, x), x, meijerg=True
).rewrite(expint).expand() == x*expint(1, x) - exp(-x)
assert integrate(expint(2, x), x, meijerg=True
).rewrite(expint).expand() == \
-x**2*expint(1, x)/2 + x*exp(-x)/2 - exp(-x)/2
assert simplify(unpolarify(integrate(expint(y, x), x,
meijerg=True).rewrite(expint).expand(func=True))) == \
-expint(y + 1, x)
assert integrate(Si(x), x, meijerg=True) == x*Si(x) + cos(x)
assert integrate(Ci(u), u, meijerg=True).expand() == u*Ci(u) - sin(u)
assert integrate(Shi(x), x, meijerg=True) == x*Shi(x) - cosh(x)
assert integrate(Chi(u), u, meijerg=True).expand() == u*Chi(u) - sinh(u)
assert integrate(Si(x)*exp(-x), (x, 0, oo), meijerg=True) == pi/4
assert integrate(expint(1, x)*sin(x), (x, 0, oo), meijerg=True) == log(2)/2
def reciprocal_frame_test():
Print_Function()
metric = '1 # #,' + \
'# 1 #,' + \
'# # 1,'
(e1, e2, e3) = MV.setup('e1 e2 e3', metric)
print('g_{ij} =', MV.metric)
E = e1 ^ e2 ^ e3
Esq = (E*E).scalar()
print('E =', E)
print('%E^{2} =', Esq)
Esq_inv = 1/Esq
E1 = (e2 ^ e3)*E
E2 = (-1)*(e1 ^ e3)*E
E3 = (e1 ^ e2)*E
print('E1 = (e2^e3)*E =', E1)
print('E2 =-(e1^e3)*E =', E2)
print('E3 = (e1^e2)*E =', E3)
w = (E1 | e2)
w = w.expand()
print('E1|e2 =', w)
w = (E1 | e3)
w = w.expand()
print('E1|e3 =', w)
w = (E2 | e1)
w = w.expand()
print('E2|e1 =', w)
w = (E2 | e3)
w = w.expand()
print('E2|e3 =', w)
w = (E3 | e1)
w = w.expand()
print('E3|e1 =', w)
w = (E3 | e2)
w = w.expand()
print('E3|e2 =', w)
w = (E1 | e1)
w = (w.expand()).scalar()
Esq = expand(Esq)
print('%(E1\\cdot e1)/E^{2} =', simplify(w/Esq))
w = (E2 | e2)
w = (w.expand()).scalar()
print('%(E2\\cdot e2)/E^{2} =', simplify(w/Esq))
w = (E3 | e3)
w = (w.expand()).scalar()
print('%(E3\\cdot e3)/E^{2} =', simplify(w/Esq))
return
def entropy(density):
"""Compute the entropy of a matrix/density object.
This computes -Tr(density*ln(density)) using the eigenvalue decomposition
of density, which is given as either a Density instance or a matrix
(numpy.ndarray, sympy.Matrix or scipy.sparse).
Parameters
==========
density : density matrix of type Density, sympy matrix,
scipy.sparse or numpy.ndarray
Examples
========
>>> from sympy.physics.quantum.density import Density, entropy
>>> from sympy.physics.quantum.represent import represent
>>> from sympy.physics.quantum.matrixutils import scipy_sparse_matrix
>>> from sympy.physics.quantum.spin import JzKet, Jz
>>> from sympy import S, log
>>> up = JzKet(S(1)/2,S(1)/2)
>>> down = JzKet(S(1)/2,-S(1)/2)
>>> d = Density((up,0.5),(down,0.5))
>>> entropy(d)
log(2)/2
"""
if isinstance(density, Density):
density = represent(density) # represent in Matrix
if isinstance(density, scipy_sparse_matrix):
density = to_numpy(density)
if isinstance(density, Matrix):
eigvals = density.eigenvals().keys()
return expand(-sum(e*log(e) for e in eigvals))
elif isinstance(density, numpy_ndarray):
import numpy as np
eigvals = np.linalg.eigvals(density)
return -np.sum(eigvals*np.log(eigvals))
else:
raise ValueError(
"numpy.ndarray, scipy.sparse or sympy matrix expected")
