Python Matrix.jacobian方法代码示例

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def getj_0j_1():
    a1, a2 = dynamicsymbols('a1, a2')

    mass_matrix = kane.mass_matrix
    forcing = kane.forcing
    tf_0 = Matrix([[cos(theta1), -sin(theta1), -leg_length*sin(theta1)], [sin(theta1), cos(theta1), leg_length*cos(theta1)], [0,0,1]])

    tf_1 = Matrix([[cos(theta2), -sin(theta2), -body_length*sin(theta2)], [sin(theta2), cos(theta2), body_length*cos(theta2)], [0,0,1]])
        
    tf_01 = simplify(tf_0*tf_1)
    
    com_0 = Matrix([-leg_length*sin(theta1), leg_length*cos(theta1), 1])
    
    com_1 = Matrix([-body_length*sin(theta2), body_length*cos(theta2), 1])
    
    com_01 = simplify(tf_0*com_1)
    
    com = Matrix(com_0 + com_01)
    
    j_0 = com_0.jacobian([theta1, theta2])
    j_1 = simplify(com_01.jacobian([theta1, theta2]))
    j = com.jacobian([theta1, theta2])
    
    q = Matrix(coordinates)
    qdot = Matrix(speeds)
    qddot = Matrix([a1, a2])
    return [j_0, j_1]

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def test_jacobian_metrics():
    rho, phi = symbols("rho phi")
    X = Matrix([rho*cos(phi), rho*sin(phi)])
    Y = Matrix([rho, phi])
    J = X.jacobian(Y)
    assert J == X.jacobian(Y.T)
    assert J == (X.T).jacobian(Y)
    assert J == (X.T).jacobian(Y.T)
    g = J.T*eye(J.shape[0])*J
    g = g.applyfunc(trigsimp)
    assert g == Matrix([[1, 0], [0, rho**2]])

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def jacobian_dHdV(nK, y, cap, inds):
	# Set up equations of power flow (power at line from nodal voltage) as symbolic equations and
	# calculate the corresponding Jacobian matrix.
	from sympy import symbols, Matrix

	g = Matrix(np.real(y))
	b = Matrix(np.imag(y))

	e_J = symbols("e1:%d"%(nK+1))
	f_J = symbols("f1:%d"%(nK+1))
	hSluV = []

	if "Pl" in inds:
		nPl = np.asarray(inds["Pl"]).astype(int)
		for k in range(nPl.shape[0]):
			i,j = nPl[k,:]
			hSluV.append((e_J[i]**2+f_J[i]**2)*g[i,j]-(e_J[i]*e_J[j]+f_J[i]*f_J[j])*g[i,j]+(e_J[i]*f_J[j]-e_J[j]*f_J[i])*b[i,j])
	if "Ql" in inds:
		nQl = np.asarray(inds["Ql"]).astype(int)
		for k in range(nQl.shape[0]):
			i,j = nQl[k,:]
			hSluV.append(-(e_J[i]**2+f_J[i]**2)*b[i,j]+(e_J[i]*e_J[j]+f_J[i]*f_J[j])*b[i,j]+(e_J[i]*f_J[j]-e_J[j]*f_J[i])*g[i,j]-(e_J[i]**2+f_J[i]**2)*cap[i,j]/2)
	if "Vm" in inds:
		nVk = np.asarray(inds["Vm"]).astype(int)
		for k in range(nVk.shape[0]):
			i = nVk[k]
			hSluV.append((e_J[i]**2+f_J[i]**2)**0.5)
	if len(hSluV)>0:
		hSluV = Matrix(hSluV)
		ef_J = e_J[1:] + f_J[1:]
		JMatrix_dHdV = hSluV.jacobian(ef_J)
		return JMatrix_dHdV, hSluV
	else:
		return None, None

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def calculate_dependencies(equations, task_variables, variables_values):
    
    replacements = replace.get_replacements(task_variables)[1]
    transformed_equations = replace.get_transformed_expressions(equations, task_variables)
    
    #prepare data for Numpy manipulations - change variables names to the names of special format - 'x%sy', also 
    #convert some strings to Symbols so they could be processed by Sympy and Numpy libs
        
    symbolized_replaced_task_variables = sorted(symbols(replacements.keys()))        
    replaced_variables_values = {}
    for k,v in replacements.iteritems():
        replaced_variables_values[Symbol(k)] = variables_values[v]
        
    #prepare the system of equations for solving - find Jacobian, in Jacobian matrix that was found replace variables by their values
    
    F = Matrix(transformed_equations)
    J = F.jacobian(symbolized_replaced_task_variables).evalf(subs=replaced_variables_values)
    X = np.transpose(np.array(J))
    
    #solve the system
    n,m = X.shape
    if (n != m - 1): 
        # if the system of equations is overdetermined (the number of equations is bigger than number of variables 
        #then we have to make regular, square systems of equations out of it (we do that by grabbing only some of the equations)
        square_systems = suggest_square_systems_list(X, n, m)
        for sq in square_systems:
            try:
                solve_system(sq, equations)
            except Exception as e:
                print e.args
                print sq
                
                
    else:
        solve_system(X, equations)  

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def main():
    u, v, R = symbols('u v R', real=True)
    xi, eta = symbols(r'\xi \eta', cls=Function)

    numer = 4*R**2
    denom = u**2 + v**2 + numer

    # inverse of a stereographic projection from the south pole
    # onto the XY plane:
    pinv = Matrix([numer * u / denom,
                   numer * v / denom,
                   -(2 * R * (u**2 + v**2)) / denom]) # OK
    if False:
        # textbook style
        Dpinv = simplify(pinv.jacobian([u, v]))
        print_latex(Dpinv, mat_str='pmatrix', mat_delim=None) # OK?

        tDpinvDpinv = factor(Dpinv.transpose() @ Dpinv)
        print_latex(tDpinvDpinv, mat_str='pmatrix', mat_delim=None) # OK

        tDpinvDpinv = tDpinvDpinv.subs([(u, xi(t)), (v, eta(t))])
        dcdt = Matrix([xi(t).diff(), eta(t).diff()])
        print_latex(simplify(
            sqrt((dcdt.transpose() @ tDpinvDpinv).dot(dcdt))))
    else:
        # directly 
        dpinvc = pinv.subs([(u, xi(t)), (v, eta(t))]).diff(t, 1)
        print_latex(sqrt(factor(dpinvc.dot(dpinvc))))

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def test_jacobian_hessian():
    x = Symbol('x')
    y = Symbol('y')
    L = Matrix(1,2,[x**2*y, 2*y**2 + x*y])
    syms = [x,y]
    assert L.jacobian(syms) == Matrix([[2*x*y, x**2],[y, 4*y+x]])

    L = Matrix(1,2,[x, x**2*y**3])
    assert L.jacobian(syms) == Matrix([[1, 0], [2*x*y**3, x**2*3*y**2]])

    f = x**2*y
    syms = [x,y]
    assert hessian(f, syms) == Matrix([[2*y, 2*x], [2*x, 0]])

    f = x**2*y**3
    assert hessian(f, syms) == Matrix([[2*y**3, 6*x*y**2],[6*x*y**2, 6*x**2*y]])

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
    def compute_formal_matrices(self):
        model = self.model
        exo_eqs = [eq for eq in model.equations if eq.info.get('exogenous') == 'true']
        non_exo_eqs = [eq for eq in model.equations if not eq in exo_eqs]
        exo_vars = [eq.lhs for eq in exo_eqs]
        non_exo_vars = [v for v in model.variables if not v in exo_vars]
        model.info['exo_vars'] = exo_vars
        model.info['non_exo_vars'] = non_exo_vars

        mat_exo_vars_f = Matrix([v(+1) for v in exo_vars]).T
        mat_exo_vars_c = Matrix([v for v in exo_vars]).T
        mat_exo_vars_p = Matrix([v(-1) for v in exo_vars]).T

        mat_non_exo_vars_f = Matrix( [v(+1) for v in non_exo_vars] ).T
        mat_non_exo_vars_c = Matrix( [v for v in non_exo_vars] ).T
        mat_non_exo_vars_p = Matrix( [v(-1) for v in non_exo_vars] ).T

        # Compute matrix for exogenous equations
        mat_exo_rhs = Matrix([eq.rhs for eq in exo_eqs]).T
        N = mat_exo_rhs.jacobian(mat_exo_vars_p).T

        # Compute matrices for non exogenous equations

        mat_non_exo_eqs = Matrix( [ eq.gap for eq in non_exo_eqs ] ).T
        F = mat_non_exo_eqs.jacobian(mat_non_exo_vars_f).T
        G = mat_non_exo_eqs.jacobian(mat_non_exo_vars_c).T
        H = mat_non_exo_eqs.jacobian(mat_non_exo_vars_p).T
        L = mat_non_exo_eqs.jacobian(mat_exo_vars_f).T
        M = mat_non_exo_eqs.jacobian(mat_exo_vars_c).T

        def steady_state_ify(m):
            # replaces all variables in m by steady state value
            for v in model.variables + model.shocks: # slow and inefficient
                m = m.subs(v(+1),v.P)
                m = m.subs(v(-1),v.P)
            return m

        d = dict()
        d['F'] = steady_state_ify(F)
        d['G'] = steady_state_ify(G)
        d['H'] = steady_state_ify(H)
        d['L'] = steady_state_ify(L)
        d['M'] = steady_state_ify(M)
        d['N'] = steady_state_ify(N)
        return d

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def test_jacobian2():
    rho, phi = symbols("rho phi")
    X = Matrix([rho*cos(phi), rho*sin(phi), rho**2])
    Y = Matrix([rho, phi])
    J = Matrix([
            [cos(phi), -rho*sin(phi)],
            [sin(phi),  rho*cos(phi)],
            [   2*rho,             0],
        ])
    assert X.jacobian(Y) == J

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def jacobian_dHdV(nK, y, cap, inds):
	# Set up equations of power flow (power at line from nodal voltage) as symbolic equations and
	# calculate the corresponding Jacobian matrix.
	from sympy import symbols, Matrix

	g = Matrix(np.real(y))
	b = Matrix(np.imag(y))

	e_J = symbols("e1:%d"%(nK+1))
	f_J = symbols("f1:%d"%(nK+1))
	hSluV = []

	if isinstance(inds["Pl"],np.ndarray):
		nPl = inds["Pl"].astype(int)
	elif isinstance(inds["Pl"],list):
		nPl = inds["Pl"]
	else:
		raise ValueError
	for k in range(len(nPl)):
		i,j = nPl[k,:]
		hSluV.append((e_J[i]**2+f_J[i]**2)*g[i,j]-(e_J[i]*e_J[j]+f_J[i]*f_J[j])*g[i,j]+(e_J[i]*f_J[j]-e_J[j]*f_J[i])*b[i,j])

	if isinstance(inds["Ql"],np.ndarray):
		nQl = inds["Ql"].astype(int)
	elif isinstance(inds["Ql"],list):
		nQl = inds["Ql"]
	else:
		raise ValueError
	for k in range(len(nQl)):
		i,j = nQl[k,:]
		hSluV.append(-(e_J[i]**2+f_J[i]**2)*b[i,j]+(e_J[i]*e_J[j]+f_J[i]*f_J[j])*b[i,j]+(e_J[i]*f_J[j]-e_J[j]*f_J[i])*g[i,j]-(e_J[i]**2+f_J[i]**2)*cap[i,j]/2)

	if isinstance(inds["Vm"],np.ndarray):
		nVk = inds["Vm"].astype(int)
	elif isinstance(inds["Vm"],list):
		nVk = inds["Vm"]
	else:
		raise ValueError
	for k in range(len(nVk)):
		i = nVk[k]
		hSluV.append((e_J[i]**2+f_J[i]**2)**0.5)
	hSluV = Matrix(hSluV)
	ef_J = e_J[1:] + f_J[1:]
	JMatrix_dHdV = hSluV.jacobian(ef_J)

	return JMatrix_dHdV, hSluV

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def get_difeomorphisms01(n,s,which_difeo,ini_tabix,end_tabix,a_index,b_index,asserts):
	import time

	correct_difeomorphism = []
	R = []	
	s2 = []	
	Z = []
	B = []
	B,Z = n_grados001(n)
	N = matrix_gtdie(Z)
	#print ('trabajamos con: ' + str(size(s)))
	tabix = 0
	
	start_time = time.time()
	for s1 in s:
		if tabix < ini_tabix:
			tabix = tabix + 1
			continue
			
		am_i_good = 0
		Z0 = which_difeo(Z,B,s1)
		#D = compatibility_equation(Z0,B,simplify,N)
		h = Matrix(Z0)
		D = h.jacobian(B).det()
		
		
		for asserto in range(0,size(asserts)):
			am_i_good = am_i_good + asserts[asserto](D,n)
			
		if am_i_good == 0:
			correct_difeomorphism.append(Z0)
			R.append(D)
			s2.append(s1)
		
		if tabix >= end_tabix:
			break
		
		tabix = tabix + 1
	print("--- %s seconds ---" % (time.time() - start_time))		
	return correct_difeomorphism, R, s2

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def main():
    x_1, x_2 = symbols('x_1 x_2', real=True)

    # Phi(x_1, x_2) = a torus.
    Phi = Matrix([(2 + cos(x_2)) * cos(x_1),
                  (2 + cos(x_2)) * sin(x_1),
                  sin(x_2),])
    DPhi = Phi.jacobian([x_1, x_2])
    print_latex(DPhi, fold_func_brackets=True,
        mat_str='pmatrix', mat_delim=None)

    tDPhi = DPhi.transpose()
    print_latex(tDPhi, fold_func_brackets=True,
        mat_str='pmatrix', mat_delim=None)

    tDPhiDPhi = simplify(tDPhi * DPhi) # simplify \sin^2 + \cos^2.
    print_latex(tDPhiDPhi, fold_func_brackets=True,
        mat_str='pmatrix', mat_delim=None)
    dcdt = Matrix([Derivative(xi, t, 1), Derivative(eta, t, 1)])

    print_latex(sqrt((dcdt.transpose() @ tDPhiDPhi).dot(dcdt)),
       fold_func_brackets=True)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
def jacobian_dSdV(Y, nK):
	# Set up equations of power flow (bus power from nodal voltage) as symbolic equations and
	# calculate the corresponding Jacobian matrix.
	from sympy import symbols, Matrix

	G = Matrix(np.real(Y))
	B = Matrix(np.imag(Y))

	e_J = symbols("e1:%d"%(nK+1))
	f_J = symbols("f1:%d"%(nK+1))

	hSK = []
	for i in range(nK):
		hSK.append(e_J[i]*(G[i,:].dot(e_J) - B[i,:].dot(f_J)) + f_J[i]*(G[i,:].dot(f_J) + B[i,:].dot(e_J)))
	for i in range(nK):
		hSK.append(f_J[i]*(G[i,:].dot(e_J) - B[i,:].dot(f_J)) - e_J[i]*(G[i,:].dot(f_J) + B[i,:].dot(e_J)))
	hSK = Matrix(hSK)
	ef_J = e_J[1:] + f_J[1:]

	Jac_equ = hSK.jacobian(ef_J)

	return Jac_equ, hSK

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]

#.........这里部分代码省略.........
        """Supply the following for the initialization of LagrangesMethod

        Lagrangian : Sympifyable

        qs: array_like
            The generalized coordinates

        hol_coneqs: array_like, optional
            The holonomic constraint equations

        nonhol_coneqs: array_like, optional
            The nonholonomic constraint equations

        forcelist : iterable, optional
            Takes an iterable 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
            and/or moments.

        frame : ReferenceFrame, optional
            Supply the inertial frame. This is used to determine the
            generalized forces due to non-conservative forces.
        """

        self._L = Matrix([sympify(Lagrangian)])
        self.eom = None
        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()      # The coeffecients of the multipliers

        forcelist = forcelist if forcelist else []
        if not iterable(forcelist):
            raise TypeError('Force pairs must be supplied in an iterable.')
        self._forcelist = forcelist
        if frame and not isinstance(frame, ReferenceFrame):
            raise TypeError('frame must be a valid ReferenceFrame')
        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
        if not iterable(qs):
            raise TypeError('Generalized coordinates must be an iterable')
        self._q = Matrix(qs)
        self._qdots = self.q.diff(dynamicsymbols._t)
        self._qdoubledots = self._qdots.diff(dynamicsymbols._t)

        # Deal with constraint equations
        if coneqs:
            SymPyDeprecationWarning("The `coneqs` kwarg is deprecated in "
                    "favor of `hol_coneqs` and `nonhol_coneqs`. Please "
                    "update your code").warn()
            self.coneqs = coneqs
        else:
            mat_build = lambda x: Matrix(x) if x else Matrix()
            hol_coneqs = mat_build(hol_coneqs)
            nonhol_coneqs = mat_build(nonhol_coneqs)
            self.coneqs = Matrix([hol_coneqs.diff(dynamicsymbols._t),
                    nonhol_coneqs])
            self._hol_coneqs = hol_coneqs

    def form_lagranges_equations(self):
        """Method to form Lagrange's equations of motion.

        Returns a vector of equations of motion using Lagrange's equations of
        the second kind.
        """

        qds = self._qdots
        qdd_zero = dict((i, 0) for i in self._qdoubledots)
        n = len(self.q)

        # Internally we represent the EOM as four terms:
        # EOM = term1 - term2 - term3 - term4 = 0

        # First term
        self._term1 = self._L.jacobian(qds)
        self._term1 = self._term1.diff(dynamicsymbols._t).T

        # Second term
        self._term2 = self._L.jacobian(self.q).T

        # Third term
        if self.coneqs:
            coneqs = self.coneqs
            m = len(coneqs)
            # Creating the multipliers
            self.lam_vec = Matrix(dynamicsymbols('lam1:' + str(m + 1)))
            self.lam_coeffs = -coneqs.jacobian(qds)
            self._term3 = self.lam_coeffs.T * self.lam_vec
            # Extracting the coeffecients of the qdds from the diff coneqs
            diffconeqs = coneqs.diff(dynamicsymbols._t)
            self._m_cd = diffconeqs.jacobian(self._qdoubl

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
# DNT============ Analytical derivatives and sstate evaluation ================

        nx = len(x);
        ny = len(y);
        nxp = len(xp);
        nyp = len(yp);
        n = len(f);


        # Array of steady state value
        xss=array(lambdify(args,x,'numpy')(*SSval)).reshape(nx,1)
        yss=array(lambdify(args,y,'numpy')(*SSval)).reshape(ny,1)


        # 1st derivatives
        fx = f.jacobian(x)
        nfx=array(lambdify(args,fx,'numpy')(*SSval))

        fxp = f.jacobian(xp)
        nfxp=array(lambdify(args,fxp,'numpy')(*SSval))

        fy = f.jacobian(y)
        nfy=array(lambdify(args,fy,'numpy')(*SSval))

        fyp = f.jacobian(yp)
        nfyp=array(lambdify(args,fyp,'numpy')(*SSval))


        # 2nd derivatives
        nfypyp=array(lambdify(args,[fyp[ii,:].jacobian(yp) for ii in range(n)],'numpy')(*SSval)).reshape(n,nyp,nyp)
        nfypy=array(lambdify(args,[fyp[ii,:].jacobian(y) for ii in range(n)],'numpy')(*SSval)).reshape(n,nyp,ny)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import jacobian [as 别名]
    def form_lagranges_equations(self):
        """Method to form Lagrange's equations of motion.

        Returns a vector of equations of motion using Lagrange's equations of
        the second kind.

        """

        q = self._q
        qd = self._qdots
        qdd = self._qdoubledots
        n = len(q)

        #Putting the Lagrangian in a Matrix
        L = Matrix([self._L])

        #Determining the first term in Lagrange's EOM
        self._term1 = L.jacobian(qd)
        self._term1 = ((self._term1).diff(dynamicsymbols._t)).transpose()

        #Determining the second term in Lagrange's EOM
        self._term2 = (L.jacobian(q)).transpose()

        #term1 and term2 will be there no matter what so leave them as they are

        if self.coneqs is not None:
            coneqs = self.coneqs
            #If there are coneqs supplied then the following will be created
            coneqs = list(coneqs)
            if not isinstance(coneqs, list):
                raise TypeError('Enter the constraint equations in a list')

            o = len(coneqs)

            #Creating the multipliers
            self.lam_vec = Matrix(dynamicsymbols('lam1:' + str(o + 1)))

            #Extracting the coeffecients of the multipliers
            coneqs_mat = Matrix(coneqs)
            qd = self._qdots
            self.lam_coeffs = -coneqs_mat.jacobian(qd)

            #Determining the third term in Lagrange's EOM
            #term3 = ((self.lam_vec).transpose() * self.lam_coeffs).transpose()
            self._term3 = self.lam_coeffs.transpose() * self.lam_vec

            #Taking the time derivative of the constraint equations
            diffconeqs = [diff(i, dynamicsymbols._t) for i in coneqs]

            #Extracting the coeffecients of the qdds from the diff coneqs
            diffconeqs_mat = Matrix(diffconeqs)
            qdd = self._qdoubledots
            self._m_cd = diffconeqs_mat.jacobian(qdd)

            #The remaining terms i.e. the 'forcing' terms in diff coneqs
            qddzero = dict(zip(qdd, [0] * n))
            self._f_cd = -diffconeqs_mat.subs(qddzero)

        else:
            self._term3 = zeros(n, 1)

        if self.forcelist is not None:
            forcelist = self.forcelist
            N = self.inertial
            if not isinstance(N, ReferenceFrame):
                raise TypeError('Enter a valid ReferenceFrame')
            if not isinstance(forcelist, (list, tuple)):
                raise TypeError('Forces must be supplied in a list of: lists'
                        ' or tuples')
            self._term4 = zeros(n, 1)
            for i, v in enumerate(qd):
                for j, w in enumerate(forcelist):
                    if isinstance(w[0], ReferenceFrame):
                        speed = w[0].ang_vel_in(N)
                        self._term4[i] += speed.diff(v, N) & w[1]
                    if isinstance(w[0], Point):
                        speed = w[0].vel(N)
                        self._term4[i] += speed.diff(v, N) & w[1]
                    else:
                        raise TypeError('First entry in force pair is a point'
                                        ' or frame')

        else:
            self._term4 = zeros(n, 1)

        self.eom = self._term1 - self._term2 - self._term3 - self._term4

        return self.eom