Python sympy 模块，diff() 实例源码

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 transmutagen_cram_error(degree, t0, prec=200):
    from ..partialfrac import (thetas_alphas, thetas_alphas_to_expr_complex,
        customre, t)
    from ..cram import get_CRAM_from_cache

    from sympy import re, exp, nsolve, diff

    expr = get_CRAM_from_cache(degree, prec)
    thetas, alphas, alpha0 = thetas_alphas(expr, prec)
    part_frac = thetas_alphas_to_expr_complex(thetas, alphas, alpha0)
    part_frac = part_frac.replace(customre, re)


    E = part_frac - exp(-t)
    E = E.evalf(20)

    return E.subs(t, nsolve(diff(E, t), t0, prec=prec))

def Simple_manifold_with_scalar_function_derivative():
    Print_Function()
    coords = (x, y, z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3', metric='[1,1,1]', coords=coords)
    # Define surface
    mfvar = (u, v) = symbols('u v')
    X = u*e1 + v*e2 + (u**2 + v**2)*e3
    print(X)
    MF = Manifold(X, mfvar)

    # Define field on the surface.
    g = (v + 1)*log(u)

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g, u) + (MF.rbasis[1]/MF.E_sq)*diff(g, v)
    print('Vector derivative =', VectorDerivative.subs({u: 1, v: 0}))

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print('Vector derivative =', dg.subs({u: 1, v: 0}))
    return

def Simple_manifold_with_scalar_function_derivative():
    Print_Function()
    coords = (x, y, z) = symbols('x y z')
    basis = (e1, e2, e3, grad) = MV.setup('e_1 e_2 e_3', metric='[1,1,1]', coords=coords)
    # Define surface
    mfvar = (u, v) = symbols('u v')
    X = u*e1 + v*e2 + (u**2 + v**2)*e3
    print('\\f{X}{u,v} =', X)
    MF = Manifold(X, mfvar)
    (eu, ev) = MF.Basis()
    # Define field on the surface.
    g = (v + 1)*log(u)

    print('\\f{g}{u,v} =', g)

    # Method 1: Using old Manifold routines.
    VectorDerivative = (MF.rbasis[0]/MF.E_sq)*diff(g, u) + (MF.rbasis[1]/MF.E_sq)*diff(g, v)
    print('\\eval{\\nabla g}{u=1,v=0} =', VectorDerivative.subs({u: 1, v: 0}))

    # Method 2: Using new Manifold routines.
    dg = MF.Grad(g)
    print('\\eval{\\f{Grad}{g}}{u=1,v=0} =', dg.subs({u: 1, v: 0}))
    dg = MF.grad*g
    print('\\eval{\\nabla g}{u=1,v=0} =', dg.subs({u: 1, v: 0}))
    return

def linearity_index_inverse_depth():
    """Linearity index of Inverse Depth Parameterization"""
    D, rho, rho_0, d_1, sigma_rho = sympy.symbols("D,rho,rho_0,d_1,sigma_rho")
    alpha = sympy.symbols("alpha")
    u = (rho * sympy.sin(alpha)) / (rho_0 * d_1 * (rho_0 - rho) + rho * sympy.cos(alpha))  # NOQA

    # first order derivative of u
    u_p = sympy.diff(u, rho)
    u_p = sympy.simplify(u_p)

    # second order derivative of u
    u_pp = sympy.diff(u_p, rho)
    u_pp = sympy.simplify(u_pp)

    # Linearity index
    L = (u_pp * 2 * sigma_rho) / (u_p)
    L = sympy.simplify(L)

    print()
    print("u: ", u)
    print("u': ", u_p)
    print("u'': ", u_pp)
    # print("L = ", L)
    print("L = ", L.subs(rho, 0))
    print()

def application():
    """makes list of the difference between the estimated derivative for h=0.01 and
    the analytical derivative"""
    application_list = []
    x = symbols('x')
    print
    print '%20s %20s' % ('Function', 'Error')
    print
    funcvar_tuples = [(lambda t: exp(t), 'exp(x)', float(0)), (lambda t : exp(-2*t**2), \
    'exp(-2*x**2)', float(0)), (lambda t: cos(t), 'cos(x)', 2*pi), (lambda t: log(t), \
    'log(x)', float(1))]
    for i in funcvar_tuples:
        prime = diff(i[1], x)
        exact = lambdify([x], prime)
        error = numdiff(i[0], i[2], 0.01) - exact(i[2])
        print '%20s %21.2E' % (i[1], error)
        print
        application_list.append(error)
    return application_list

def proj1(poly_set, x):
    p_out = []
    # F is a polynomial in A
    for F in poly_set:
        if degree(F, x) > 1:
            R = reducta(F, x)
            for G in R:
                if degree(G, x) > 0:
                    H = diff(G, x)
                    psc_1 = PSC(G, H, x)
                    if len(psc_1):
                        if degree(G) == 2:
                            p_out.append(poly(psc_1[0]))
                        else:
                            for aux in psc_1[0:degree(G) - 2]:
                                p_out.append(poly(aux))
    return p_out

def marginal_product_capital(self):
        r"""
        Symbolic expression for the marginal product of capital (per
        unit effective labor).

        :getter: Return the current marginal product of capital (per
        unit effective labor).
        :type: sym.Basic

        Notes
        -----
        The marginal product of capital is defined as follows:

        .. math::

            \frac{\partial F(K, AL)}{\partial K} \equiv f'(k)

        where :math:`k=K/AL` is capital stock (per unit effective labor)

        """
        return sym.diff(self.intensive_output, k)

def test_diff_and_applyfunc():
    from sympy.abc import x, y, z
    md = ImmutableDenseNDimArray([[x, y], [x*z, x*y*z]])
    assert md.diff(x) == ImmutableDenseNDimArray([[1, 0], [z, y*z]])
    assert diff(md, x) == ImmutableDenseNDimArray([[1, 0], [z, y*z]])

    sd = ImmutableSparseNDimArray(md)
    assert sd == ImmutableSparseNDimArray([x, y, x*z, x*y*z], (2, 2))
    assert sd.diff(x) == ImmutableSparseNDimArray([[1, 0], [z, y*z]])
    assert diff(sd, x) == ImmutableSparseNDimArray([[1, 0], [z, y*z]])

    mdn = md.applyfunc(lambda x: x*3)
    assert mdn == ImmutableDenseNDimArray([[3*x, 3*y], [3*x*z, 3*x*y*z]])
    assert md != mdn

    sdn = sd.applyfunc(lambda x: x/2)
    assert sdn == ImmutableSparseNDimArray([[x/2, y/2], [x*z/2, x*y*z/2]])
    assert sd != sdn

def testCalculusIntegrate(self):
    dataset_objects = algorithmic_math.math_dataset_init(
        8, digits=5, functions={"log": "L"})
    counter = 0
    for d in algorithmic_math.calculus_integrate(8, 0, 3, 10):
      counter += 1
      decoded_input = dataset_objects.int_decoder(d["inputs"])
      var, expression = decoded_input.split(":")
      target = dataset_objects.int_decoder(d["targets"])

      for fn_name, fn_char in six.iteritems(dataset_objects.functions):
        target = target.replace(fn_char, fn_name)

      # Take the derivative of the target.
      derivative = str(sympy.diff(target, var))

      # Check that the derivative of the integral equals the input.
      self.assertEqual(0, sympy.simplify("%s-(%s)" % (expression, derivative)))
    self.assertEqual(counter, 10)

def q_affine(expr, vars):
    """
    Return True if ``expr`` is (separately) affine in the variables ``vars``,
    False otherwise.

    Readapted from: https://stackoverflow.com/questions/36283548\
        /check-if-an-equation-is-linear-for-a-specific-set-of-variables/
    """
    # A function is (separately) affine in a given set of variables if all
    # non-mixed second order derivatives are identically zero.
    for x in as_tuple(vars):
        if x not in expr.atoms():
            return False
        try:
            if diff(expr, x) == nan or not Eq(diff(expr, x, x), 0):
                return False
        except TypeError:
            return False
    return True

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 differentiate(self, *, equation : str):
        '''
        Differentiate an equation
        with respect to x (dx)
        '''
        x = sympy.symbols('x')
        try:
            await self.bot.embed_reply("`{}`".format(sympy.diff(equation.strip('`'), x)), title = "Derivative of {}".format(equation))
        except Exception as e:
            await self.bot.embed_reply(py_code_block.format("{}: {}".format(type(e).__name__, e)), title = "Error")

def __call__(self, sat_pos, args=(), **kwds):
        """
        Return the definite integral of *self.fun*(pos, *args*[0],
        ..., *args*[-1]) for the line of site from *stn_pos* to
        *sat_pos* (a :class:`PyPosition`) where pos is a
        :class:`PyPosition` on the line of site (and note the
        integration bounds on h defined in __init__). The remaining
        *kwds* are passed to the quadrature routine (:py:func:`quad`).
        """
        diff = NP.array(sat_pos.xyz) - NP.array(self.stn_pos.xyz)
        S_los = NP.linalg.norm(diff) / 1e3
        def pos(s):
            """
            Return the ECEF vector a distance *s* along the line-of-site (in
            [km]).
            """
            return PyPosition(*(NP.array(self.stn_pos.xyz) + (s / S_los) * diff))
        # determine integration bounds
        # distance along of line of site at which the geodetic height
        # is self.height1
        s1 = minimize_scalar(lambda l: (pos(l).height / 1e3 - self.height1)**2,
                             bounds=[0, S_los],
                             method='Bounded').x
        # distance along of line of site at which the geodetic height
        # is self.height2
        s2 = minimize_scalar(lambda l: (pos(l).height / 1e3 - self.height2)**2,
                             bounds=[0, S_los],
                             method='Bounded').x
        def wrapper(s, *args):
            return self.fun(pos(s), *args)
        return quad(wrapper, s1, s2, args=args, **kwds)[0]

def __init__(self, stn_pos, **kwds):
        """
        Setup Chapman electron density profile F2 peak density derivative
        slant integrator.
        """
        z, Nm, Hm, H_O = SYM.symbols('z Nm Hm H_O')
        f_sym = chapman_sym_scaled(z, Nm, Hm, H_O)
        DNm_sym = SYM.diff(f_sym, Nm)
        f = SYM.lambdify((z, Hm, H_O),
                         DNm_sym,
                         modules='numexpr')
        wrapper = lambda pos, *args: f(pos.height / 1e3, *args)
        super(DNmChapmanSI, self).__init__(wrapper, stn_pos, **kwds)

def __init__(self, stn_pos, **kwds):
        """
        Setup Chapman electron density profile F2 peak height derivative
        slant integrator.
        """
        z, Nm, Hm, H_O = SYM.symbols('z Nm Hm H_O')
        f_sym = chapman_sym_scaled(z, Nm, Hm, H_O)
        DHm_sym = SYM.diff(f_sym, Hm)
        f = SYM.lambdify((z, Nm, Hm, H_O),
                         DHm_sym,
                         modules='numexpr')
        wrapper = lambda pos, *args: f(pos.height / 1e3, *args)
        super(DHmChapmanSI, self).__init__(wrapper, stn_pos, **kwds)

def __init__(self, stn_pos, **kwds):
        """
        Setup Chapman electron density profile plasma scale height
        derivative slant integrator.
        """
        z, Nm, Hm, H_O = SYM.symbols('z Nm Hm H_O')
        f_sym = chapman_sym_scaled(z, Nm, Hm, H_O)
        DH_O_sym = SYM.diff(f_sym, H_O)
        f = SYM.lambdify((z, Nm, Hm, H_O),
                         DH_O_sym,
                         modules='numexpr')
        wrapper = lambda pos, *args: f(pos.height / 1e3, *args)
        super(DH_OChapmanSI, self).__init__(wrapper, stn_pos, **kwds)

def test_tanh_sinh(f, a, b, exact):
    # test fine error estimate
    mp.dps = 50

    tol = 10**(-mp.dps)
    tol2 = 10**(-mp.dps+1)

    t = sympy.Symbol('t')
    f_derivatives = {
        1: sympy.lambdify(t, sympy.diff(f(t), t, 1), modules=['mpmath']),
        2: sympy.lambdify(t, sympy.diff(f(t), t, 2), modules=['mpmath']),
        }

    value, _ = quadpy.line_segment.tanh_sinh(
                f, a, b, tol,
                f_derivatives=f_derivatives
                )
    assert abs(value - exact) < tol2

    # test with crude estimate
    value, _ = quadpy.line_segment.tanh_sinh(f, a, b, tol)
    assert abs(value - exact) < tol2
    return


# Test functions with singularities at both ends.

def test_singularities_at_both_ends(f_left, f_right, b, exact):
    # test fine error estimate
    tol = 10**(-mp.dps)

    t = sympy.Symbol('t')
    fl = {
        0: f_left,
        1: sympy.lambdify(t, sympy.diff(f_left(t), t, 1), modules=['mpmath']),
        2: sympy.lambdify(t, sympy.diff(f_left(t), t, 2), modules=['mpmath']),
        }
    fr = {
        0: f_right,
        1: sympy.lambdify(t, sympy.diff(f_right(t), t, 1), modules=['mpmath']),
        2: sympy.lambdify(t, sympy.diff(f_right(t), t, 2), modules=['mpmath']),
        }

    value, _ = quadpy.line_segment.tanh_sinh_lr(fl, fr, b, tol)
    tol2 = 10**(-mp.dps+1)
    assert abs(value - exact) < tol2

    # # test with crude estimate
    # fl = {0: f_left}
    # fr = {0: f_right}
    # value, _ = quadpy.line_segment.tanh_sinh_lr(fl, fr, b, tol)
    # tol2 = 10**(-mp.dps+2)
    # assert abs(value - exact) < tol2
    return

def paper_cram_error(degree, t0, prec=200):
    from ..partialfrac import t, customre

    from sympy import re, exp, nsolve, diff

    paper_part_frac = get_paper_part_frac(degree).replace(customre, re)

    E = paper_part_frac - exp(-t)

    return E.subs(t, nsolve(diff(E, t), t0, prec=prec))

def diff(self, x):
        Dself = Vector()
        if isinstance(Vector.dtau_dict, dict):
            Dself.obj = linear_derivation(self.obj, Vector.Diff, x)
        else:
            Dself.obj = diff(self.obj, x)
        return Dself

def _coords(self, qind, qdep=[], coneqs=[]):
        """Supply all the generalized coordinates in a list.

        If some coordinates are dependent, supply them as part of qdep. Their
        dependent nature will only show up in the linearization process though.

        Parameters
        ==========

        qind : list
            A list of independent generalized coords
        qdep : list
            List of dependent coordinates
        coneq : list
            List of expressions which are equal to zero; these are the
            configuration constraint equations
        """

        if not isinstance(qind, (list, tuple)):
            raise TypeError('Generalized coords. must be supplied in a list.')
        self._q = qind + qdep
        self._qdot = [diff(i, dynamicsymbols._t) for i in self._q]

        if not isinstance(qdep, (list, tuple)):
            raise TypeError('Dependent coordinates and constraints must each be '
                            'provided in their own list.')
        if len(qdep) != len(coneqs):
            raise ValueError('There must be an equal number of dependent '
                             'coordinates and constraints.')
        coneqs = Matrix(coneqs)
        self._qdep = qdep
        self._f_h = coneqs

def mass_matrix_full(self):
        # Returns the mass matrix from above, augmented by kin diff's k_kqdot
        if (self._frstar is None) & (self._fr is None):
            raise ValueError('Need to compute Fr, Fr* first.')
        o = len(self._u)
        n = len(self._q)
        return ((self._k_kqdot).row_join(zeros(n, o))).col_join((zeros(o,
                n)).row_join(self.mass_matrix))

def gradient(scalar, frame):
    """
    Returns the vector gradient of a scalar field computed wrt the
    coordinate symbols of the given frame.

    Parameters
    ==========

    scalar : sympifiable
        The scalar field to take the gradient of

    frame : ReferenceFrame
        The frame to calculate the gradient in

    Examples
    ========

    >>> from sympy.physics.vector import ReferenceFrame
    >>> from sympy.physics.vector import gradient
    >>> R = ReferenceFrame('R')
    >>> s1 = R[0]*R[1]*R[2]
    >>> gradient(s1, R)
    R_y*R_z*R.x + R_x*R_z*R.y + R_x*R_y*R.z
    >>> s2 = 5*R[0]**2*R[2]
    >>> gradient(s2, R)
    10*R_x*R_z*R.x + 5*R_x**2*R.z

    """

    _check_frame(frame)
    outvec = Vector(0)
    scalar = express(scalar, frame, variables=True)
    for i, x in enumerate(frame):
        outvec += diff(scalar, frame[i]) * x
    return outvec

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)*sqrt(-cos(th)**2 + 1)*exp(-I*ph)/(4*sqrt(pi))
    assert Ynm(1, -1, th, ph).expand(func=True) == sqrt(6)*sqrt(-cos(th)**2 + 1)*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)*sqrt(-cos(th)**2 + 1)*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)*sqrt(-cos(th)**2 + 1)*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 grad(f, basis, for_numerical=True):
    """
    Compute the symbolic gradient of a vector-valued function with respect to a
    basis.

    Parameters
    ----------
    f : 1D array_like of sympy Expressions
        The vector-valued function to compute the gradient of.
    basis : 1D array_like of sympy symbols
        The basis symbols to compute the gradient with respect to.
    for_numerical : bool, optional
        A placeholder for the option of numerically computing the gradient.

    Returns
    -------
    grad : 2D array_like of sympy Expressions
        The symbolic gradient.
    """
    if hasattr(f, '__len__'):  # as of version 1.1.1, Array isn't supported
        f = sp.Matrix(f)

    return f.__class__([
        [
            sp.diff(f[x], basis[y])
            if not for_numerical or not f[x].has(sp.sign(basis[y])) else 0
            for y in range(len(basis))
        ] for x in range(len(f))
    ])

def linearity_index_xyz_parameterization():
    """Linearity index of XYZ Parameterization"""
    d, d1, alpha, sigma_d = sympy.symbols("d,d1,alpha,sigma_d")
    alpha = sympy.symbols("alpha")
    u = (d * sympy.sin(alpha)) / (d1 + d * sympy.cos(alpha))

    # first order derivative of u
    u_p = sympy.diff(u, d)
    u_p = sympy.simplify(u_p)

    # second order derivative of u
    u_pp = sympy.diff(u_p, d)
    u_pp = sympy.simplify(u_pp)

    # Linearity index
    L = (u_pp * 2 * sigma_d) / (u_p)
    L = sympy.simplify(L)

    print()
    print("u: ", u)
    print("u': ", u_p)
    print("u'': ", u_pp)
    # print("L = ", L)
    print("L = ", L.subs(d, 0))
    print()


# def derive_J_jacobian():
#     r_W_C = sympy.symbols("r_W_C")
#
#     y =

def test_Quantity_derivative():
    x = symbols("x")
    assert diff(x*meter, x) == meter
    assert diff(x**3*meter**2, x) == 3*x**2*meter**2
    assert diff(meter, meter) == 1
    assert diff(meter**2, meter) == 2*meter

def divergence(vect, frame):
    """
    Returns the divergence 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 divergence in

    Examples
    ========

    >>> from sympy.physics.vector import ReferenceFrame
    >>> from sympy.physics.vector import divergence
    >>> R = ReferenceFrame('R')
    >>> v1 = R[0]*R[1]*R[2] * (R.x+R.y+R.z)
    >>> divergence(v1, R)
    R_x*R_y + R_x*R_z + R_y*R_z
    >>> v2 = 2*R[1]*R[2]*R.y
    >>> divergence(v2, R)
    2*R_z

    """

    _check_vector(vect)
    if vect == 0:
        return S(0)
    vect = express(vect, frame, variables=True)
    vectx = vect.dot(frame.x)
    vecty = vect.dot(frame.y)
    vectz = vect.dot(frame.z)
    out = S(0)
    out += diff(vectx, frame[0])
    out += diff(vecty, frame[1])
    out += diff(vectz, frame[2])
    return out

def gradient(scalar, frame):
    """
    Returns the vector gradient of a scalar field computed wrt the
    coordinate symbols of the given frame.

    Parameters
    ==========

    scalar : sympifiable
        The scalar field to take the gradient of

    frame : ReferenceFrame
        The frame to calculate the gradient in

    Examples
    ========

    >>> from sympy.physics.vector import ReferenceFrame
    >>> from sympy.physics.vector import gradient
    >>> R = ReferenceFrame('R')
    >>> s1 = R[0]*R[1]*R[2]
    >>> gradient(s1, R)
    R_y*R_z*R.x + R_x*R_z*R.y + R_x*R_y*R.z
    >>> s2 = 5*R[0]**2*R[2]
    >>> gradient(s2, R)
    10*R_x*R_z*R.x + 5*R_x**2*R.z

    """

    _check_frame(frame)
    outvec = Vector(0)
    scalar = express(scalar, frame, variables=True)
    for i, x in enumerate(frame):
        outvec += diff(scalar, frame[i]) * x
    return outvec

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))'

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