def test_matrixelement_diff():
dexpr = diff((D*w)[k,0], w[p,0])
assert w[k, p].diff(w[k, p]) == 1
assert w[k, p].diff(w[0, 0]) == KroneckerDelta(0, k)*KroneckerDelta(0, p)
assert str(dexpr) == "Sum(KroneckerDelta(_k, p)*D[k, _k], (_k, 0, n - 1))"
assert str(dexpr.doit()) == 'Piecewise((D[k, p], (0 <= p) & (p <= n - 1)), (0, True))'
python类diff()的实例源码
def test_MatrixElement_with_values():
x, y, z, w = symbols("x y z w")
M = Matrix([[x, y], [z, w]])
i, j = symbols("i, j")
Mij = M[i, j]
assert isinstance(Mij, MatrixElement)
Ms = SparseMatrix([[2, 3], [4, 5]])
msij = Ms[i, j]
assert isinstance(msij, MatrixElement)
for oi, oj in [(0, 0), (0, 1), (1, 0), (1, 1)]:
assert Mij.subs({i: oi, j: oj}) == M[oi, oj]
assert msij.subs({i: oi, j: oj}) == Ms[oi, oj]
A = MatrixSymbol("A", 2, 2)
assert A[0, 0].subs(A, M) == x
assert A[i, j].subs(A, M) == M[i, j]
assert M[i, j].subs(M, A) == A[i, j]
assert isinstance(M[3*i - 2, j], MatrixElement)
assert M[3*i - 2, j].subs({i: 1, j: 0}) == M[1, 0]
assert isinstance(M[i, 0], MatrixElement)
assert M[i, 0].subs(i, 0) == M[0, 0]
assert M[0, i].subs(i, 1) == M[0, 1]
assert M[i, j].diff(x) == Matrix([[1, 0], [0, 0]])[i, j]
raises(ValueError, lambda: M[i, 2])
raises(ValueError, lambda: M[i, -1])
raises(ValueError, lambda: M[2, i])
raises(ValueError, lambda: M[-1, i])
def test_diff():
from sympy.abc import x, y, z
md = MutableDenseNDimArray([[x, y], [x*z, x*y*z]])
assert md.diff(x) == MutableDenseNDimArray([[1, 0], [z, y*z]])
assert diff(md, x) == MutableDenseNDimArray([[1, 0], [z, y*z]])
sd = MutableSparseNDimArray(md)
assert sd == MutableSparseNDimArray([x, y, x*z, x*y*z], (2, 2))
assert sd.diff(x) == MutableSparseNDimArray([[1, 0], [z, y*z]])
assert diff(sd, x) == MutableSparseNDimArray([[1, 0], [z, y*z]])
def test_Ynm():
# http://en.wikipedia.org/wiki/Spherical_harmonics
th, ph = Symbol("theta", real=True), Symbol("phi", real=True)
from sympy.abc import n,m
assert Ynm(0, 0, th, ph).expand(func=True) == 1/(2*sqrt(pi))
assert Ynm(1, -1, th, ph) == -exp(-2*I*ph)*Ynm(1, 1, th, ph)
assert Ynm(1, -1, th, ph).expand(func=True) == sqrt(6)*sin(th)*exp(-I*ph)/(4*sqrt(pi))
assert Ynm(1, -1, th, ph).expand(func=True) == sqrt(6)*sin(th)*exp(-I*ph)/(4*sqrt(pi))
assert Ynm(1, 0, th, ph).expand(func=True) == sqrt(3)*cos(th)/(2*sqrt(pi))
assert Ynm(1, 1, th, ph).expand(func=True) == -sqrt(6)*sin(th)*exp(I*ph)/(4*sqrt(pi))
assert Ynm(2, 0, th, ph).expand(func=True) == 3*sqrt(5)*cos(th)**2/(4*sqrt(pi)) - sqrt(5)/(4*sqrt(pi))
assert Ynm(2, 1, th, ph).expand(func=True) == -sqrt(30)*sin(th)*exp(I*ph)*cos(th)/(4*sqrt(pi))
assert Ynm(2, -2, th, ph).expand(func=True) == (-sqrt(30)*exp(-2*I*ph)*cos(th)**2/(8*sqrt(pi))
+ sqrt(30)*exp(-2*I*ph)/(8*sqrt(pi)))
assert Ynm(2, 2, th, ph).expand(func=True) == (-sqrt(30)*exp(2*I*ph)*cos(th)**2/(8*sqrt(pi))
+ sqrt(30)*exp(2*I*ph)/(8*sqrt(pi)))
assert diff(Ynm(n, m, th, ph), th) == (m*cot(th)*Ynm(n, m, th, ph)
+ sqrt((-m + n)*(m + n + 1))*exp(-I*ph)*Ynm(n, m + 1, th, ph))
assert diff(Ynm(n, m, th, ph), ph) == I*m*Ynm(n, m, th, ph)
assert conjugate(Ynm(n, m, th, ph)) == (-1)**(2*m)*exp(-2*I*m*ph)*Ynm(n, m, th, ph)
assert Ynm(n, m, -th, ph) == Ynm(n, m, th, ph)
assert Ynm(n, m, th, -ph) == exp(-2*I*m*ph)*Ynm(n, m, th, ph)
assert Ynm(n, -m, th, ph) == (-1)**m*exp(-2*I*m*ph)*Ynm(n, m, th, ph)
def test_fd_space(derivative, space_order):
"""
This test compare the discrete finite-difference scheme against polynomials
For a given order p, the fiunite difference scheme should
be exact for polynomials of order p
:param derivative: name of the derivative to be tested
:param space_order: space order of the finite difference stencil
"""
clear_cache()
# dummy axis dimension
nx = 100
xx = np.linspace(-1, 1, nx)
dx = xx[1] - xx[0]
# Symbolic data
grid = Grid(shape=(nx,), dtype=np.float32)
x = grid.dimensions[0]
u = Function(name="u", grid=grid, space_order=space_order)
du = Function(name="du", grid=grid, space_order=space_order)
# Define polynomial with exact fd
coeffs = np.ones((space_order,), dtype=np.float32)
polynome = sum([coeffs[i]*x**i for i in range(0, space_order)])
polyvalues = np.array([polynome.subs(x, xi) for xi in xx], np.float32)
# Fill original data with the polynomial values
u.data[:] = polyvalues
# True derivative of the polynome
Dpolynome = diff(diff(polynome)) if derivative == 'dx2' else diff(polynome)
Dpolyvalues = np.array([Dpolynome.subs(x, xi) for xi in xx], np.float32)
# FD derivative, symbolic
u_deriv = getattr(u, derivative)
# Compute numerical FD
stencil = Eq(du, u_deriv)
op = Operator(stencil, subs={x.spacing: dx})
op.apply()
# Check exactness of the numerical derivative except inside space_brd
space_border = space_order
error = abs(du.data[space_border:-space_border] -
Dpolyvalues[space_border:-space_border])
assert np.isclose(np.mean(error), 0., atol=1e-3)
def buildJacobianMatrix(self, including_fast_reactions=True):
self.interactionMatrix = []
self.jacobianMatrix = []
for ode in self.ODEs:
partial_ders = []
signs = []
for var in self.ODE_vars:
t_part_der = diff(ode.getDeveloppedInternalMathFormula(), var.getDeveloppedInternalMathFormula())
t_part_der_formula = MathFormula(self)
t_part_der_formula.setInternalMathFormula(t_part_der)
partial_ders.append(t_part_der_formula)
if t_part_der == MathFormula.ZERO:
signs.append(0)
else:
for atom in t_part_der.atoms(SympySymbol):
t_part_der = t_part_der.subs({atom: SympyInteger(1)})
if t_part_der > MathFormula.ZERO:
signs.append(1)
elif t_part_der < MathFormula.ZERO:
signs.append(-1)
else:
signs.append(2)
self.jacobianMatrix.append(partial_ders)
self.interactionMatrix.append(signs)
def chapman_fit(alt,
ne,
x0=[1e6, 300, 50],
bounds=[(1, None),
(150, 500),
(30, 80)],
verbose=False,
**kwds):
"""
"""
# optimization setup
z, Nm, Hm, H_O = SYM.symbols('z Nm Hm H_O')
chapman = chapman_sym(z, Nm, Hm, H_O)
dNm = SYM.diff(chapman, Nm)
dHm = SYM.diff(chapman, Hm)
dH_O = SYM.diff(chapman, H_O)
chapman_f = SYM.lambdify((z, Nm, Hm, H_O),
chapman,
modules='numexpr')
dNm_f = SYM.lambdify((z, Nm, Hm, H_O),
dNm,
modules='numexpr')
dHm_f = SYM.lambdify((z, Nm, Hm, H_O),
dHm,
modules='numexpr')
dH_O_f = SYM.lambdify((z, Nm, Hm, H_O),
dH_O,
modules='numexpr')
# define cost function
y = NP.asarray(ne)
def J(x):
Nm, Hm, H_O = x
if verbose:
print('-' * 80)
print(x)
y_hat = NP.array([chapman_f(z, Nm, Hm, H_O) for z in alt])
diff = y - y_hat
J1 = NP.array([dNm_f(z, Nm, Hm, H_O) for z in alt])
J2 = NP.array([dHm_f(z, Nm, Hm, H_O) for z in alt])
J3 = NP.array([dH_O_f(z, Nm, Hm, H_O) for z in alt])
return (NP.dot(diff, diff),
NP.array([-2 * NP.sum(diff * J1),
-2 * NP.sum(diff * J2),
-2 * NP.sum(diff * J3)]))
# minimize cost function
x_star, f, d = scipy.optimize.fmin_l_bfgs_b(J,
x0,
bounds=bounds,
**kwds)
assert d['warnflag'] == 0
return x_star
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 __init__(self, Lagrangian, q_list, coneqs=None, forcelist=None, frame=None):
"""Supply the following for the initialization of LagrangesMethod
Lagrangian : Sympifyable
q_list : list
A list of the generalized coordinates
coneqs : list
A list of the holonomic and non-holonomic constraint equations.
VERY IMPORTANT NOTE- The holonomic constraints must be
differentiated with respect to time and then included in coneqs.
forcelist : list
Takes a list of (Point, Vector) or (ReferenceFrame, Vector) tuples
which represent the force at a point or torque on a frame. This
feature is primarily to account for the nonconservative forces
amd/or moments.
frame : ReferenceFrame
Supply the inertial frame. This is used to determine the
generalized forces due to non-sonservative forces.
"""
self._L = sympify(Lagrangian)
self.eom = None # initializing the eom Matrix
self._m_cd = Matrix([]) # Mass Matrix of differentiated coneqs
self._m_d = Matrix([]) # Mass Matrix of dynamic equations
self._f_cd = Matrix([]) # Forcing part of the diff coneqs
self._f_d = Matrix([]) # Forcing part of the dynamic equations
self.lam_coeffs = Matrix([]) # Initializing the coeffecients of lams
self.forcelist = forcelist
self.inertial = frame
self.lam_vec = Matrix([])
self._term1 = Matrix([])
self._term2 = Matrix([])
self._term3 = Matrix([])
self._term4 = Matrix([])
# Creating the qs, qdots and qdoubledots
q_list = list(q_list)
if not isinstance(q_list, list):
raise TypeError('Generalized coords. must be supplied in a list')
self._q = q_list
self._qdots = [diff(i, dynamicsymbols._t) for i in self._q]
self._qdoubledots = [diff(i, dynamicsymbols._t) for i in self._qdots]
self.coneqs = coneqs
def curl(vect, frame):
"""
Returns the curl of a vector field computed wrt the coordinate
symbols of the given frame.
Parameters
==========
vect : Vector
The vector operand
frame : ReferenceFrame
The reference frame to calculate the curl in
Examples
========
>>> from sympy.physics.vector import ReferenceFrame
>>> from sympy.physics.vector import curl
>>> R = ReferenceFrame('R')
>>> v1 = R[1]*R[2]*R.x + R[0]*R[2]*R.y + R[0]*R[1]*R.z
>>> curl(v1, R)
0
>>> v2 = R[0]*R[1]*R[2]*R.x
>>> curl(v2, R)
R_x*R_y*R.y - R_x*R_z*R.z
"""
_check_vector(vect)
if vect == 0:
return Vector(0)
vect = express(vect, frame, variables=True)
#A mechanical approach to avoid looping overheads
vectx = vect.dot(frame.x)
vecty = vect.dot(frame.y)
vectz = vect.dot(frame.z)
outvec = Vector(0)
outvec += (diff(vectz, frame[1]) - diff(vecty, frame[2])) * frame.x
outvec += (diff(vectx, frame[2]) - diff(vectz, frame[0])) * frame.y
outvec += (diff(vecty, frame[0]) - diff(vectx, frame[1])) * frame.z
return outvec
def scalar_potential(field, frame):
"""
Returns the scalar potential function of a field in a given frame
(without the added integration constant).
Parameters
==========
field : Vector
The vector field whose scalar potential function is to be
calculated
frame : ReferenceFrame
The frame to do the calculation in
Examples
========
>>> from sympy.physics.vector import ReferenceFrame
>>> from sympy.physics.vector import scalar_potential, gradient
>>> R = ReferenceFrame('R')
>>> scalar_potential(R.z, R) == R[2]
True
>>> scalar_field = 2*R[0]**2*R[1]*R[2]
>>> grad_field = gradient(scalar_field, R)
>>> scalar_potential(grad_field, R)
2*R_x**2*R_y*R_z
"""
#Check whether field is conservative
if not is_conservative(field):
raise ValueError("Field is not conservative")
if field == Vector(0):
return S(0)
#Express the field exntirely in frame
#Subsitute coordinate variables also
_check_frame(frame)
field = express(field, frame, variables=True)
#Make a list of dimensions of the frame
dimensions = [x for x in frame]
#Calculate scalar potential function
temp_function = integrate(field.dot(dimensions[0]), frame[0])
for i, dim in enumerate(dimensions[1:]):
partial_diff = diff(temp_function, frame[i + 1])
partial_diff = field.dot(dim) - partial_diff
temp_function += integrate(partial_diff, frame[i + 1])
return temp_function
def curl(vect, frame):
"""
Returns the curl of a vector field computed wrt the coordinate
symbols of the given frame.
Parameters
==========
vect : Vector
The vector operand
frame : ReferenceFrame
The reference frame to calculate the curl in
Examples
========
>>> from sympy.physics.vector import ReferenceFrame
>>> from sympy.physics.vector import curl
>>> R = ReferenceFrame('R')
>>> v1 = R[1]*R[2]*R.x + R[0]*R[2]*R.y + R[0]*R[1]*R.z
>>> curl(v1, R)
0
>>> v2 = R[0]*R[1]*R[2]*R.x
>>> curl(v2, R)
R_x*R_y*R.y - R_x*R_z*R.z
"""
_check_vector(vect)
if vect == 0:
return Vector(0)
vect = express(vect, frame, variables=True)
#A mechanical approach to avoid looping overheads
vectx = vect.dot(frame.x)
vecty = vect.dot(frame.y)
vectz = vect.dot(frame.z)
outvec = Vector(0)
outvec += (diff(vectz, frame[1]) - diff(vecty, frame[2])) * frame.x
outvec += (diff(vectx, frame[2]) - diff(vectz, frame[0])) * frame.y
outvec += (diff(vecty, frame[0]) - diff(vectx, frame[1])) * frame.z
return outvec
def scalar_potential(field, frame):
"""
Returns the scalar potential function of a field in a given frame
(without the added integration constant).
Parameters
==========
field : Vector
The vector field whose scalar potential function is to be
calculated
frame : ReferenceFrame
The frame to do the calculation in
Examples
========
>>> from sympy.physics.vector import ReferenceFrame
>>> from sympy.physics.vector import scalar_potential, gradient
>>> R = ReferenceFrame('R')
>>> scalar_potential(R.z, R) == R[2]
True
>>> scalar_field = 2*R[0]**2*R[1]*R[2]
>>> grad_field = gradient(scalar_field, R)
>>> scalar_potential(grad_field, R)
2*R_x**2*R_y*R_z
"""
#Check whether field is conservative
if not is_conservative(field):
raise ValueError("Field is not conservative")
if field == Vector(0):
return S(0)
#Express the field exntirely in frame
#Subsitute coordinate variables also
_check_frame(frame)
field = express(field, frame, variables=True)
#Make a list of dimensions of the frame
dimensions = [x for x in frame]
#Calculate scalar potential function
temp_function = integrate(field.dot(dimensions[0]), frame[0])
for i, dim in enumerate(dimensions[1:]):
partial_diff = diff(temp_function, frame[i + 1])
partial_diff = field.dot(dim) - partial_diff
temp_function += integrate(partial_diff, frame[i + 1])
return temp_function
def derive_by_array(expr, dx):
r"""
Derivative by arrays. Supports both arrays and scalars.
Given the array `A_{i_1, \ldots, i_N}` and the array `X_{j_1, \ldots, j_M}`
this function will return a new array `B` defined by
`B_{j_1,\ldots,j_M,i_1,\ldots,i_N} := \frac{\partial A_{i_1,\ldots,i_N}}{\partial X_{j_1,\ldots,j_M}}`
Examples
========
>>> from sympy import derive_by_array
>>> from sympy.abc import x, y, z, t
>>> from sympy import cos
>>> derive_by_array(cos(x*t), x)
-t*sin(t*x)
>>> derive_by_array(cos(x*t), [x, y, z, t])
[-t*sin(t*x), 0, 0, -x*sin(t*x)]
>>> derive_by_array([x, y**2*z], [[x, y], [z, t]])
[[[1, 0], [0, 2*y*z]], [[0, y**2], [0, 0]]]
"""
from sympy.matrices import MatrixBase
array_types = (collections.Iterable, MatrixBase, NDimArray)
if isinstance(dx, array_types):
dx = ImmutableDenseNDimArray(dx)
for i in dx:
if not i._diff_wrt:
raise ValueError("cannot derive by this array")
if isinstance(expr, array_types):
expr = ImmutableDenseNDimArray(expr)
if isinstance(dx, array_types):
new_array = [[y.diff(x) for y in expr] for x in dx]
return type(expr)(new_array, dx.shape + expr.shape)
else:
return expr.diff(dx)
else:
if isinstance(dx, array_types):
return ImmutableDenseNDimArray([expr.diff(i) for i in dx], dx.shape)
else:
return diff(expr, dx)
def scalar_potential(field, coord_sys):
"""
Returns the scalar potential function of a field in a given
coordinate system (without the added integration constant).
Parameters
==========
field : Vector
The vector field whose scalar potential function is to be
calculated
coord_sys : CoordSysCartesian
The coordinate system to do the calculation in
Examples
========
>>> from sympy.vector import CoordSysCartesian
>>> from sympy.vector import scalar_potential, gradient
>>> R = CoordSysCartesian('R')
>>> scalar_potential(R.k, R) == R.z
True
>>> scalar_field = 2*R.x**2*R.y*R.z
>>> grad_field = gradient(scalar_field, R)
>>> scalar_potential(grad_field, R)
2*R.x**2*R.y*R.z
"""
# Check whether field is conservative
if not is_conservative(field):
raise ValueError("Field is not conservative")
if field == Vector.zero:
return S(0)
# Express the field exntirely in coord_sys
# Subsitute coordinate variables also
if not isinstance(coord_sys, CoordSysCartesian):
raise TypeError("coord_sys must be a CoordSysCartesian")
field = express(field, coord_sys, variables=True)
dimensions = coord_sys.base_vectors()
scalars = coord_sys.base_scalars()
# Calculate scalar potential function
temp_function = integrate(field.dot(dimensions[0]), scalars[0])
for i, dim in enumerate(dimensions[1:]):
partial_diff = diff(temp_function, scalars[i + 1])
partial_diff = field.dot(dim) - partial_diff
temp_function += integrate(partial_diff, scalars[i + 1])
return temp_function
