Python Matrix.subs方法代码示例

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
    def _kindiffeq(self, kdeqs):
        """Supply all the kinematic differential equations in a list.

        They should be in the form [Expr1, Expr2, ...] where Expri is equal to
        zero

        Parameters
        ==========

        kdeqs : list (of Expr)
            The listof kinematic differential equations

        """
        if len(self._q) != len(kdeqs):
            raise ValueError('There must be an equal number of kinematic '
                             'differential equations and coordinates.')

        uaux = self._uaux
        # dictionary of auxiliary speeds which are equal to zero
        uaz = dict(zip(uaux, [0] * len(uaux)))

        kdeqs = Matrix(kdeqs).subs(uaz)

        qdot = self._qdot
        qdotzero = dict(zip(qdot, [0] * len(qdot)))
        u = self._u
        uzero = dict(zip(u, [0] * len(u)))

        f_k = kdeqs.subs(uzero).subs(qdotzero)
        k_kqdot = (kdeqs.subs(uzero) - f_k).jacobian(Matrix(qdot))
        k_ku = (kdeqs.subs(qdotzero) - f_k).jacobian(Matrix(u))

        self._k_ku = self._mat_inv_mul(k_kqdot, k_ku)
        self._f_k = self._mat_inv_mul(k_kqdot, f_k)
        self._k_kqdot = eye(len(qdot))

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
    def _initialize_kindiffeq_matrices(self, kdeqs):
        """Initialize the kinematic differential equation matrices."""

        if kdeqs:
            if len(self.q) != len(kdeqs):
                raise ValueError('There must be an equal number of kinematic '
                                 'differential equations and coordinates.')
            kdeqs = Matrix(kdeqs)

            u = self.u
            qdot = self._qdot
            # Dictionaries setting things to zero
            u_zero = dict((i, 0) for i in u)
            uaux_zero = dict((i, 0) for i in self._uaux)
            qdot_zero = dict((i, 0) for i in qdot)

            f_k = kdeqs.subs(u_zero).subs(qdot_zero)
            k_ku = (kdeqs.subs(qdot_zero) - f_k).jacobian(u)
            k_kqdot = (kdeqs.subs(u_zero) - f_k).jacobian(qdot)

            f_k = k_kqdot.LUsolve(f_k)
            k_ku = k_kqdot.LUsolve(k_ku)
            k_kqdot = eye(len(qdot))

            self._qdot_u_map = solve_linear_system_LU(
                    Matrix([k_kqdot.T, -(k_ku * u + f_k).T]).T, qdot)

            self._f_k = f_k.subs(uaux_zero)
            self._k_ku = k_ku.subs(uaux_zero)
            self._k_kqdot = k_kqdot
        else:
            self._qdot_u_map = None
            self._f_k = Matrix()
            self._k_ku = Matrix()
            self._k_kqdot = Matrix()

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
    def _mat_inv_mul(self, A, B):
        """Internal Function

        Computes A^-1 * B symbolically w/ substitution, where B is not
        necessarily a vector, but can be a matrix.

        """

        # Note: investigate difficulty in only creating symbols for non-zero
        # entries; this could speed things up, perhaps?

        r1, c1 = A.shape
        r2, c2 = B.shape
        temp1 = Matrix(r1, c1, lambda i, j: Symbol('x' + str(j + r1 * i)))
        temp2 = Matrix(r2, c2, lambda i, j: Symbol('y' + str(j + r2 * i)))
        for i in range(len(temp1)):
            if A[i] == 0:
                temp1[i] = 0
        for i in range(len(temp2)):
            if B[i] == 0:
                temp2[i] = 0
        temp3 = []
        for i in range(c2):
            temp3.append(temp1.LUsolve(temp2.extract(range(r2), [i])))
        temp3 = Matrix([i.T for i in temp3]).T
        if Kane.simp == True:
            temp3.simplify()
        return temp3.subs(dict(zip(temp1, A))).subs(dict(zip(temp2, B)))

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [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 subs [as 别名]
def test_has():
    x, y, z = symbols('x,y,z')
    A = Matrix(((x,y),(2,3)))
    assert A.has(x)
    assert not A.has(z)
    assert A.has(Symbol)

    A = A.subs(x,2)
    assert not A.has(x)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
    def get_jacobian(self, eqs, at=None):
        jac = Matrix(eqs).jacobian(Matrix(self.phis))

        if not at is None:
            assert len(self.phis) == len(at)
            for p, v in zip(self.phis, at):
                jac = jac.subs(p, v)

        return jac

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
def mintegrate(m, var, l1, l2, mname, sufix, norm=False, do_simplify=False,
               conds='none'):
    print('\tstarting integration of {mname} over {var}'.format(
            mname=mname, var=var))
    m = Matrix(m)
    sympy.var('AAA')
    if norm:
        # changing variables
        subs = {var: (l2-l1)*AAA + l1}
        m = m.subs(subs)
        m *= (l2-l1)
    # integration
    for (i, j), ki in np.ndenumerate(m):
        tmp = '{mname}_{sufix}[{i}, {j}] over {var}'.format(mname=mname,
                  sufix=sufix, i=i, j=j, var=var)
        try:
            if norm:
                ki = integrate(ki, (AAA, 0, 1), conds=conds)
            else:
                ki = integrate(ki, (var, l1, l2), conds=conds)
            print('\tfinished integrate {tmp}'.format(tmp=tmp))
        except:
            print('\t\tintegrate() failed for {tmp}'.format(tmp=tmp))
        m[i, j] = ki
    for (i, j), ki in np.ndenumerate(m):
        tmp = '{mname}_{sufix}[{i}, {j}] over {var}'.format(mname=mname,
                  sufix=sufix, i=i, j=j, var=var)
        try:
            if do_simplify:
                ki = simplify(ki)
            else:
                ki = trigsimp(ki)
        except:
            print('\t\ttrigsimp failed {tmp}'.format(tmp=tmp))
            ki = simplify(ki)
        print('\tfinished simplify {tmp}'.format(tmp=tmp))
        m[i, j] = ki
    # printing
    filename = 'print_{mname}_{sufix}_over_{var}.txt'.format(mname=mname,
                   sufix=sufix, var=var)
    with open(filename, 'w') as f:
        def myprint(sth):
            lines.append(str(sth).strip() + '\n')
        myprint('matrix ' + mname + ' in file ' + filename)
        for (i, j), v in np.ndenumerate(m):
            if v:
                myprint(mname + '[{0},{1}] = {2}'.format(i,j,v))
    return m

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
def _mat_inv_mul(A, B):
    """
    Computes A^-1 * B symbolically w/ substitution, where B is not
    necessarily a vector, but can be a matrix.

    """

    r1, c1 = A.shape
    r2, c2 = B.shape
    temp1 = Matrix(r1, c1, lambda i, j: Symbol('x' + str(j) + str(r1 * i)))
    temp2 = Matrix(r2, c2, lambda i, j: Symbol('y' + str(j) + str(r2 * i)))
    for i in range(len(temp1)):
        if A[i] == 0:
            temp1[i] = 0
    for i in range(len(temp2)):
        if B[i] == 0:
            temp2[i] = 0
    temp3 = []
    for i in range(c2):
        temp3.append(temp1.LDLsolve(temp2[:, i]))
    temp3 = Matrix([i.T for i in temp3]).T
    return temp3.subs(dict(list(zip(temp1, A)))).subs(dict(list(zip(temp2, B))))

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

#.........这里部分代码省略.........
        self._u = Matrix([u_ind, u_dep])
        self._udot = self.u.diff(dynamicsymbols._t)
        self._uaux = none_handler(u_aux)

    def _initialize_constraint_matrices(self, config, vel, acc):
        """Initializes constraint matrices."""

        # Define vector dimensions
        o = len(self.u)
        m = len(self._udep)
        p = o - m
        none_handler = lambda x: Matrix(x) if x else Matrix()

        # Initialize configuration constraints
        config = none_handler(config)
        if len(self._qdep) != len(config):
            raise ValueError('There must be an equal number of dependent '
                             'coordinates and configuration constraints.')
        self._f_h = none_handler(config)

        # Initialize velocity and acceleration constraints
        vel = none_handler(vel)
        acc = none_handler(acc)
        if len(vel) != m:
            raise ValueError('There must be an equal number of dependent '
                             'speeds and velocity constraints.')
        if acc and (len(acc) != m):
            raise ValueError('There must be an equal number of dependent '
                             'speeds and acceleration constraints.')
        if vel:
            u_zero = dict((i, 0) for i in self.u)
            udot_zero = dict((i, 0) for i in self._udot)

            self._f_nh = vel.subs(u_zero)
            self._k_nh = (vel - self._f_nh).jacobian(self.u)
            # If no acceleration constraints given, calculate them.
            if not acc:
                self._f_dnh = (self._k_nh.diff(dynamicsymbols._t) * self.u +
                               self._f_nh.diff(dynamicsymbols._t))
                self._k_dnh = self._k_nh
            else:
                self._f_dnh = acc.subs(udot_zero)
                self._k_dnh = (acc - self._f_dnh).jacobian(self._udot)

            # Form of non-holonomic constraints is B*u + C = 0.
            # We partition B into independent and dependent columns:
            # Ars is then -B_dep.inv() * B_ind, and it relates dependent speeds
            # to independent speeds as: udep = Ars*uind, neglecting the C term.
            B_ind = self._k_nh[:, :p]
            B_dep = self._k_nh[:, p:o]
            self._Ars = -B_dep.LUsolve(B_ind)
        else:
            self._f_nh = Matrix()
            self._k_nh = Matrix()
            self._f_dnh = Matrix()
            self._k_dnh = Matrix()
            self._Ars = Matrix()

    def _initialize_kindiffeq_matrices(self, kdeqs):
        """Initialize the kinematic differential equation matrices."""

        if kdeqs:
            if len(self.q) != len(kdeqs):
                raise ValueError('There must be an equal number of kinematic '
                                 'differential equations and coordinates.')
            kdeqs = Matrix(kdeqs)

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

#.........这里部分代码省略.........
    PaperForkCgX                    =  0.9
    PaperForkCgZ                    =  0.7
    FrameLength                     =  evalf.N(PaperWb*sin(HTA)-(rake-(PaperRadFront-PaperRadRear)*cos(HTA)))
    FrameCGNorm                     =  evalf.N((PaperFrameCgZ - PaperRadRear-(PaperFrameCgX/sin(HTA))*cos(HTA))*sin(HTA))
    FrameCGPar                      =  evalf.N((PaperFrameCgX / sin(HTA) + (PaperFrameCgZ - PaperRadRear - PaperFrameCgX / sin(HTA) * cos(HTA)) * cos(HTA)))
    tempa                           =  evalf.N((PaperForkCgZ - PaperRadFront))
    tempb                           =  evalf.N((PaperWb-PaperForkCgX))
    tempc                           =  evalf.N(sqrt(tempa**2+tempb**2))
    PaperForkL                      =  evalf.N((PaperWb*cos(HTA)-(PaperRadFront-PaperRadRear)*sin(HTA)))
    ForkCGNorm                      =  evalf.N(rake+(tempc * sin(pi/2-HTA-acos(tempa/tempc))))
    ForkCGPar                       =  evalf.N(tempc * cos((pi/2-HTA)-acos(tempa/tempc))-PaperForkL)

    # Here is the final assembly of the numerical values. The symbol 'v' is the
    # forward speed of the bicycle (a concept which only makes sense in the
    # upright, static equilibrium case?). These are in a dictionary which will
    # later be substituted in. Again the sign on the *product* of inertia
    # values is flipped here, due to different orientations of coordinate
    # systems.
    v = symbols('v')
    val_dict = {WFrad: PaperRadFront,
                WRrad: PaperRadRear,
                htangle: HTA,
                forkoffset: rake,
                forklength: PaperForkL,
                framelength: FrameLength,
                forkcg1: ForkCGPar,
                forkcg3: ForkCGNorm,
                framecg1: FrameCGNorm,
                framecg3: FrameCGPar,
                Iwr11: 0.0603,
                Iwr22: 0.12,
                Iwf11: 0.1405,
                Iwf22: 0.28,
                Ifork11: 0.05892,
                Ifork22: 0.06,
                Ifork33: 0.00708,
                Ifork31: 0.00756,
                Iframe11: 9.2,
                Iframe22: 11,
                Iframe33: 2.8,
                Iframe31: -2.4,
                mfork: 4,
                mframe: 85,
                mwf: 3,
                mwr: 2,
                g: 9.81,
                q1: 0,
                q2: 0,
                q4: 0,
                q5: 0,
                u1: 0,
                u2: 0,
                u3: v / PaperRadRear,
                u4: 0,
                u5: 0,
                u6: v / PaperRadFront}

    # Linearizes the forcing vector; the equations are set up as MM udot =
    # forcing, where MM is the mass matrix, udot is the vector representing the
    # time derivatives of the generalized speeds, and forcing is a vector which
    # contains both external forcing terms and internal forcing terms, such as
    # centripital or coriolis forces.  This actually returns a matrix with as
    # many rows as *total* coordinates and speeds, but only as many columns as
    # independent coordinates and speeds.

    forcing_lin = KM.linearize()[0]

    # As mentioned above, the size of the linearized forcing terms is expanded
    # to include both q's and u's, so the mass matrix must have this done as
    # well.  This will likely be changed to be part of the linearized process,
    # for future reference.
    MM_full = KM.mass_matrix_full

    MM_full_s = MM_full.subs(val_dict)
    forcing_lin_s = forcing_lin.subs(KM.kindiffdict()).subs(val_dict)


    MM_full_s = MM_full_s.evalf()
    forcing_lin_s = forcing_lin_s.evalf()

    # Finally, we construct an "A" matrix for the form xdot = A x (x being the
    # state vector, although in this case, the sizes are a little off). The
    # following line extracts only the minimum entries required for eigenvalue
    # analysis, which correspond to rows and columns for lean, steer, lean
    # rate, and steer rate.
    Amat = MM_full_s.inv() * forcing_lin_s
    A = Amat.extract([1, 2, 4, 6], [1, 2, 3, 5])

    # Precomputed for comparison
    Res = Matrix([[               0,                                           0,                  1.0,                    0],
                  [               0,                                           0,                    0,                  1.0],
                  [9.48977444677355, -0.891197738059089*v**2 - 0.571523173729245, -0.105522449805691*v, -0.330515398992311*v],
                  [11.7194768719633,   -1.97171508499972*v**2 + 30.9087533932407,   3.67680523332152*v,  -3.08486552743311*v]])


    # Actual eigenvalue comparison
    eps = 1.e-12
    for i in range(6):
        error = Res.subs(v, i) - A.subs(v, i)
        assert all(abs(x) < eps for x in error)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
def test_Matrix2():
    m = Matrix([[x, x**2], [5, 2/x]])
    assert (matrix(m.subs(x, 2)) == matrix([[2, 4], [5, 1]])).all()
    m = Matrix([[sin(x), x**2], [5, 2/x]])
    assert (matrix(m.subs(x, 2)) == matrix([[sin(2), 4], [5, 1]])).all()

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
# Velocity level constraints                                (Table 1)
f_v = Matrix([dot(P.vel(N), uv) for uv in C])
f_v_dq = f_v.jacobian(q)
f_v_du = f_v.jacobian(u)

# Kinematic differential equations
kindiffs = Matrix([dot(w_c_n_qd - C.ang_vel_in(N), uv) for uv in B] +
                  [dot(v_co_n_qd - CO.vel(N), uv) for uv in N])
qdots = solve(kindiffs, qd)

B.set_ang_vel(N, w_b_n_qd.subs(qdots))
C.set_ang_acc(N, C.ang_vel_in(N).dt(B) + cross(B.ang_vel_in(N),
                                               C.ang_vel_in(N)))

# f_0 and f_1                                               (Table 1)
f_0 = kindiffs.subs(u_zero)
f_1 = kindiffs.subs(qd_zero)

# Acceleration level constraints                            (Table 1)
f_a = f_v.diff(t)
# Alternatively, f_a can be formed as
#v_co_n = cross(C.ang_vel_in(N), CO.pos_from(P))
#a_co_n = v_co_n.dt(B) + cross(B.ang_vel_in(N), v_co_n)
#f_a = Matrix([((a_co_n - CO.acc(N)) & uv).expand() for uv in B])

# Kane's dynamic equations via elbow grease
# Partial angular velocities and velocities
partial_w_C = [C.ang_vel_in(N).diff(ui, N) for ui in u]
partial_v_CO = [CO.vel(N).diff(ui, N) for ui in u]

# Active forces

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import subs [as 别名]
def dilute_solution_model(structure, e0, vac_defs, antisite_defs, T,
        trial_chem_pot = None, generate='plot'):

    """
    Compute the defect densities using dilute solution model.

    Args:
        structure: pymatgen.core.structure.Structure object representing the
            primitive or unitcell of the crystal.
        e0: The total energy of the undefected system.
            This is E0 from VASP calculation.
        vac_defs: List of vacancy defect parameters in the dictionary format.
            The keys of the dict associated with each vacancy defect are
            1) site_index, 2) site_specie, 3) site_multiplicity, and
            4) energy. 1-3 can be obtained from
            pymatgen.analysis.defects.point_defects.Vacancy class.
            Site index is expected to start with 1 (fortran index).
        antisite_defs: List of antisite defect parameters in the dictionary
            format. The keys of the dict associated with each antisite defect
            are 1) site_index, 2) site_specie, 3) site_multiplicity,
            4) substitution_specie, and 5) energy. 1-3 can be obtained
            from pymatgen.analysis.defects.point_defects.Vacancy class.
        T: Temperature in Kelvin
        trial_chem_pot (optional): Trial chemical potentials to speedup
            the plot generation. Format is {el1:mu1,...}
        generate (string): Options are plot or energy
            Chemical potentials are also returned with energy option.
            If energy option is not chosen, plot is generated.

    Returns:
        If generate=plot, the plot data is generated and returned in
        HighCharts format.
        If generate=energy, defect formation enthalpies and chemical
        potentials are returned.
    """

    if not check_input(vac_defs):
        raise ValueError('Vacancy energy is not defined')
    if not check_input(antisite_defs):
        raise ValueError('Antisite energy is not defined')

    formation_energies = {}
    formation_energies['vacancies'] = copy.deepcopy(vac_defs)
    formation_energies['antisites'] = copy.deepcopy(antisite_defs)
    for vac in formation_energies['vacancies']:
        del vac['energy']
    for asite in formation_energies['antisites']:
        del asite['energy']
    # Setup the system
    site_species = [vac_def['site_specie'] for vac_def in vac_defs]
    multiplicity = [vac_def['site_multiplicity'] for vac_def in vac_defs]
    m = len(set(site_species))      # distinct species
    n = len(vac_defs)           # inequivalent sites

    # Reduce the system and associated parameters such that only distinctive
    # atoms are retained
    comm_div = gcd(*tuple(multiplicity))
    multiplicity = [val/comm_div for val in multiplicity]
    e0 = e0/comm_div
    T = Integer(T)

    c0 = np.diag(multiplicity)
    #print ('c0',c0)
    mu = [Symbol('mu'+i.__str__()) for i in range(m)]

    # Generate maps for hashing
    # Generate specie->mu map and use it for site->mu map
    specie_order = []       # Contains hash for site->mu map    Eg: [Al, Ni]
    site_specie_set = set()             # Eg: {Ni, Al}
    for i in range(n):
        site_specie  = site_species[i]
        if site_specie not in site_specie_set:
            site_specie_set.add(site_specie)
            specie_order.append(site_specie)
    site_mu_map = []     # Eg: [mu0,mu0,mu0,mu1] where mu0->Al, and mu1->Ni
    for i in range(n):
        site_specie  = site_species[i]
        j = specie_order.index(site_specie)
        site_mu_map.append(j)
    specie_site_index_map = []      # Eg: [(0,3),(3,4)] for Al & Ni
    for i in range(m):
        low_ind = site_species.index(specie_order[i])
        if i < m-1:
            hgh_ind = site_species.index(specie_order[i+1])
        else:
            hgh_ind = n
        specie_site_index_map.append((low_ind,hgh_ind))


    """
    dC: delta concentration matrix:
    dC[i,j,k]: Concentration change of atom i, due to presence of atom
    j on lattice site k
    Special case is [i,i,i] which is considered as vacancy
    Few cases: dC[i,i,i] = -1 due to being vacancy special case
                dC[k,k,i] = +1 due to increment in k at i lattice if i
                               lattice type is of different element
                dC[i,k,i] = -1 due to decrement of ith type atom due to
                presence of kth type atom on ith sublattice and kth type
                atom specie is different from ith sublattice atom specie
#.........这里部分代码省略.........

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

# Maybe I'm calculating the Jacobian wrong? Sympy has a Jacobian function. I also just realized that those angles have nothing to do with the actual Euler angle of the hand. Since the most important component of the hand orientation is the direction the palm of the hand faces, we will take a local palm vector, which corresponds to a point on the palm, and rotate it into world space. This is nice because all the units will match, allowing us to use Euclidean distance.

# <codecell>

local_palm_vector = Matrix([Symbol("palm_x"), Symbol("palm_y"), Symbol("palm_z"), 1])
global_palm_vector = R_0[7] * local_palm_vector

target2 = Matrix([Symbol("target_%s" % s) for s in ["x", "y", "z", "palm_x", "palm_y", "palm_z"]])

hand_pose2 = Matrix([R_0[7][0,3],
                     R_0[7][1,3],
                     R_0[7][2,3]]).col_join(global_palm_vector[:3,0])

hand_pose2 = hand_pose2.subs([("shoulder_%s" % axis, dummy_shoulder_mat[row, 3])
                                  for row, axis in enumerate(["x","y","z"])])
hand_pose2 = subs_trig_exprs(hand_pose2)
J2 = hand_pose2.jacobian(
        [Symbol("theta_%d" % i) for i in range(1,8)])
hand_dist = (Matrix(hand_pose2)-target2).norm()
delta_pose2 = delta_s * (Matrix(hand_pose2)-target2)
delta_config2 = J2.T * delta_pose2
new_config = Matrix([Symbol("theta_%d" % i) for i in range(1,8)]) + delta_config2

# <codecell>

step_config_terms, step_config_list = cse(new_config, optimizations='basic')
step_config_terms_js = jsify_terms(step_config_terms)
step_config_list_js = jsify_list(step_config_list[0])
extract_parameters(new_config)

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

#.........这里部分代码省略.........
        self._f_nh = Matrix([])
        self._k_dnh = Matrix([])
        self._f_dnh = Matrix([])

        self._coords(q_ind, q_dependent, configuration_constraints)
        self._speeds(u_ind, u_dependent, velocity_constraints,
                acceleration_constraints, u_auxiliary)
        if kd_eqs is not None:
            self._kindiffeq(kd_eqs)

    def _find_dynamicsymbols(self, inlist, insyms=[]):
        """Finds all non-supplied dynamicsymbols in the expressions."""
        from sympy.core.function import AppliedUndef, Derivative
        t = dynamicsymbols._t
        return reduce(set.union, [set([i]) for j in inlist
            for i in j.atoms(AppliedUndef, Derivative)
            if i.atoms() == set([t])], set()) - insyms

        temp_f = set().union(*[i.atoms(AppliedUndef) for i in inlist])
        temp_d = set().union(*[i.atoms(Derivative) for i in inlist])
        set_f = set([a for a in temp_f if a.args == (t,)])
        set_d = set([a for a in temp_d if ((a.args[0] in set_f) and all([i == t
                     for i in a.variables]))])
        return list(set.union(set_f, set_d) - set(insyms))

    def _find_othersymbols(self, inlist, insyms=[]):
        """Finds all non-dynamic symbols in the expressions."""
        return list(reduce(set.union, [i.atoms(Symbol) for i in inlist]) -
                    set(insyms))

    def _mat_inv_mul(self, A, B):
        """Internal Function

        Computes A^-1 * B symbolically w/ substitution, where B is not
        necessarily a vector, but can be a matrix.

        """

        r1, c1 = A.shape
        r2, c2 = B.shape
        temp1 = Matrix(r1, c1, lambda i, j: Symbol('x' + str(j) + str(r1 * i)))
        temp2 = Matrix(r2, c2, lambda i, j: Symbol('y' + str(j) + str(r2 * i)))
        for i in range(len(temp1)):
            if A[i] == 0:
                temp1[i] = 0
        for i in range(len(temp2)):
            if B[i] == 0:
                temp2[i] = 0
        temp3 = []
        for i in range(c2):
            temp3.append(temp1.LDLsolve(temp2[:, i]))
        temp3 = Matrix([i.T for i in temp3]).T
        return temp3.subs(dict(list(zip(temp1, A)))).subs(dict(list(zip(temp2, B))))

    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