python类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