Python Matrix.zeros方法代码示例

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def generalized_active_forces_V(V, q, u, kde_map, vc_map=None):
    """Returns a list of the generalized active forces using equation 5.1.18
    from Kane 1985.

    'V' is a potential energy function for the system.
    'q' is a list of generalized coordinates.
    'u' is a list of the independent generalized speeds.
    'kde_map' is a dictionary with q dots as keys and the equivalent
        expressions in terms of q's and u's as values.
    'vc_map' is a dictionay with the dependent u's as keys and the expression
        in terms of independent u's as values.
    """
    n = len(q)
    p = len(u)
    m = n - p
    if vc_map is None:
        A_kr = Matrix.zeros(m, p)
    else:
        A_kr, _ = vc_matrix(u, vc_map)
        u_ = u + sorted(vc_map.keys(), cmp=lambda x, y: x.compare(y))
    W_sr, _ = kde_matrix(u_, kde_map)

    dV_dq = map(lambda x: diff(V, x), q)
    Fr = Matrix.zeros(1, p)
    for s in range(n):
        Fr -= dV_dq[s] * (W_sr[s, :p] + W_sr[s, p:]*A_kr[:, :p])
    return Fr[:]

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def generalized_inertia_forces_K(K, q, u, kde_map, vc_map=None):
    """Returns a list of the generalized inertia forces using equation 5.6.6
    from Kane 1985.

    'K' is a potential energy function for the system.
    'q' is a list of generalized coordinates.
    'u' is a list of the independent generalized speeds.
    'kde_map' is a dictionary with q dots as keys and the equivalent
        expressions in terms of q's and u's as values.
    'vc_map' is a dictionay with the dependent u's as keys and the expression
        in terms of independent u's as values.
    """
    n = len(q)
    p = len(u)
    m = n - p
    t = symbols('t')
    if vc_map is None:
        A_kr = Matrix.zeros(m, p)
    else:
        A_kr, _ = vc_matrix(u, vc_map)
        u_ = u + sorted(vc_map.keys(), cmp=lambda x, y: x.compare(y))
    W_sr, _ = kde_matrix(u_, kde_map)

    qd = map(lambda x: x.diff(t), q)
    K_partial_term = [K.diff(qd_s).diff(t) - K.diff(q_s)
                      for qd_s, q_s in zip(qd, q)]
    K_partial_term = subs(K_partial_term, kde_map)
    if vc_map is not None:
        K_partial_term = subs(K_partial_term, vc_map)
    Fr_star = Matrix.zeros(1, p)
    for s in range(n):
        Fr_star -= K_partial_term[s] * (W_sr[s, :p] + W_sr[s, p:]*A_kr[:, :p])
    return Fr_star[:]

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def LeastSquares(Problem, n):
    A = Matrix.zeros(n)
    base_f = [BasicFunction(i, Problem.a, Problem.b, Problem.alpha, Problem.beta, Problem.gamma, Problem.delta) for i in range(1, n + 1)]
    for i in range(n):
        for j in range(n):
            A[i,j] = ScalarMultiply(Problem.Operator, base_f[i], Problem.Operator, base_f[j], Problem.a, Problem.b)
    b = Matrix.zeros(n, 1)
    for i in range(n):
        b[i,0] = ScalarMultiply(0, Problem.RightPart, Problem.Operator, base_f[i], Problem.a, Problem.b)
    c = A.solve_least_squares(b)
    sol = Solution(Problem.a, Problem.b, Problem.alpha, Problem.beta, Problem.gamma, Problem.delta, *c)
    return sol

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def decompose(m,k):
    m1_const = math.factorial(k)
    m2_const = int(math.pow(2,k)-1)*math.factorial(k-1)
    m1 = Matrix.zeros((m.rows,m.cols))
    m2 = Matrix.zeros((m.rows,m.cols))
    for i in xrange(m.rows):
        for j in xrange(m.cols):
            if m[i,j] == m1_const:
                   m1[i,j] = m[i,j]
            if m[i,j] == m2_const:
                   m2[i,j] = m[i,j]
    return m1,m2

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def Collocation(Problem, n):
    A = Matrix.zeros(n)
    base_f = [BasicFunction(i, Problem.a, Problem.b, Problem.alpha, Problem.beta, Problem.gamma, Problem.delta) for i in range(1, n + 1)]
    x_arr = [Problem.a + i * (Problem.b - Problem.a) / (n + 1) for i in range(1, n + 1)]
    for i in range(n):
        for j in range(n):
            A[i,j] = Problem.Operator(base_f[j],x_arr[i])
    b = Matrix.zeros(n, 1)
    for i in range(n):
        b[i,0] = Problem.RightPart(x_arr[i])
    c = A.LUsolve(b)
    sol = Solution(Problem.a, Problem.b, Problem.alpha, Problem.beta, Problem.gamma, Problem.delta, *c)
    return sol

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def str_to_vec(_l, _symbols): #_sym_symbols are eg. myMBS.dynamics
    if not type(_symbols) == str:
        _str_symbols = str(_symbols).replace("[","").replace("(t)","").replace("]","").replace(", t", "")
    else:
        _str_symbols = _symbols.replace("[","").replace("(t)","").replace("]","").replace(", t", "")
    rows = len(_l)
    cols = 1
    m = Matrix.zeros(rows,cols)
    for i in range(len(_l)):
        m[i,0] = str_to_dynexpr(_l[i], _str_symbols)
    return m

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def str_to_matrix(_l, _symbols): #_sym_symbols are eg. myMBS.dynamics
    if not type(_symbols) == str:
        _str_symbols = str(_symbols).replace("[","").replace("(t)","").replace("]","").replace(", t", "")
    else:
        _str_symbols = _symbols.replace("[","").replace("(t)","").replace("]","").replace(", t", "")
    rows = cols = np_sqrt(len(_l))
    m = Matrix.zeros(rows,cols)
    for i in range(len(_l)):
        c = int(i%cols)
        r = int(i/cols)
        m[r,c] = str_to_dynexpr(_l[i], _str_symbols)
    return m

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def eqResistorMesh(filename, Rtest=1000.0):
	# Import mesh from file
	with open(filename, 'r') as f:
		# First line has all the node names
		V = f.readline().split()
		
		# Create neighbor lists for each node from the connections in the file
		adj = {v: {} for v in V}
		for line in f.readlines():
			edge = line.split()
			
			# If two nodes are connected by two (or more) parallel resistors, simplify
			if edge[1] in adj[edge[0]]:
				adj[edge[0]][edge[1]] = 1.0 / ( 1.0/adj[edge[0]][edge[1]] + 1.0/float(edge[2]) )
				adj[edge[1]][edge[0]] = adj[edge[0]][edge[1]]
			else:
				adj[edge[0]][edge[1]] = float(edge[2])
				adj[edge[1]][edge[0]] = float(edge[2])
	
	# Connect a test node to the first node with a known series resistance
	V.sort()
	adj['TEST'] = {V[0]: Rtest}
	adj[V[0]]['TEST'] = Rtest
	V.insert(0, 'TEST')
	
	# Generate a matrix that will be used to solve for node voltages
	N = len(adj)
	A = Matrix.zeros(N, N+1)
	
	# Test node has a known voltage (of 1V)
	A[0, 0] = 1.0
	A[0, N] = 1.0
	
	# Last node has a known voltage (of 0V)
	A[N-1, N-1] = 1.0
	
	# Populate matrix using KCL: sum of all currents in/out of a node = 0
	for i in range(1, N-1):
		for key in adj[V[i]]:
			A[i, V.index(key)] = -1.0/adj[V[i]][key]
			A[i, i] += 1.0/adj[V[i]][key]
	
	# Get matrix in reduced row echelon form, get the node voltages,
	# and calculate Req with Va=I*Req where the I is the current through the test resistor
	solved = A.rref()[0].col(N)
	return (Rtest*solved[1])/(solved[0]-solved[1])

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def getJeffreysPrior(rv):
    """
    Uses SymPy to determine the Jeffreys prior of a random variable analytically.

    Args:
        rv: SymPy RandomSymbol, corresponding to a probability distribution

    Returns:
        list: List containing Jeffreys prior in symbolic form and corresponding lambda function

    Example:
        rate = Symbol('rate', positive=True)
        rv = stats.Exponential('exponential', rate)
        print getJeffreysPrior(rv)

        >>> (1/rate, <function <lambda> at 0x0000000007F79AC8>)
    """

    # get support of random variable
    support = rv._sorted_args[0].distribution.set

    # get list of free parameters
    parameters = list(rv._sorted_args[0].distribution.free_symbols)
    x = abc.x

    # symbolic probability density function
    symPDF = density(rv)(x)

    # compute Fisher information matrix
    dim = len(parameters)
    G = Matrix.zeros(dim, dim)

    func = summation if support.is_iterable else integrate
    for i in range(0, dim):
        for j in range(0, dim):
            G[i, j] = func(simplify(symPDF *
                                    diff(ln(symPDF), parameters[i]) *
                                    diff(ln(symPDF), parameters[j])),
                           (x, support.inf, support.sup))

    # symbolic Jeffreys prior
    symJeff = simplify(sqrt(G.det()))

    # return symbolic Jeffreys prior and corresponding lambda function
    return symJeff, lambdify(parameters, symJeff, 'numpy')

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def inertia_coefficient_matrix(system, partials):
    """Returns the inertia coefficient matrix for a system of RigidBody's
    and Particle's. Each entry in the matrix, m_rs, is the inertia
    coefficient for generalized speeds r, s.

    'system' is a list where the elements are instances of RigidBody
        or Particle.
    'partials' is an instance of a PartialVelocity.

    Note: The order of the inertia coefficients is dependent on the order
    of the generalized speeds used when calculating partial velocities.
    """
    ulist = partials.ulist
    M = Matrix.zeros(len(ulist))

    for i, r in enumerate(ulist):
        for j, s in enumerate(ulist[i:], i):
            for p in system:
                m_rs = inertia_coefficient_contribution(p, partials, r, s)
                M[i, j] += m_rs
            if i != j:
                M[j, i] = M[i, j]
    return M

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
def potential_energy(Fr, q, u, kde_map, vc_map=None):
    """Returns a potential energy function using the method from Section 5.1
    from Kane 1985.

    'Fr' is a list of the generalized active forces for the system.
    'q' is a list of generalized coordinates.
    'u' is a list of the independent generalized speeds.
    'kde_map' is a dictionary with q dots as keys and the equivalent
        expressions in terms of q's and u's as values.
    'vc_map' is a dictionay with the dependent u's as keys and the expression
        in terms of independent u's as values.
    """
    n = len(q)
    p = len(u)
    m = n - p
    if vc_map is None:
        A_kr = Matrix.zeros(m, p)
    else:
        A_kr, _ = vc_matrix(u, vc_map)
        u_ = u + sorted(vc_map.keys(), cmp=lambda x, y: x.compare(y))
    W_sr, _ = kde_matrix(u_, kde_map)

    dV_dq = symbols('∂V/∂q1:{0}'.format(n + 1))
    dV_eq = Matrix(Fr).T

    for s in range(n):
        dV_eq += dV_dq[s] * (W_sr[s, :p] + W_sr[s, p:]*A_kr[:, :p])

    if vc_map is not None:
        f_arg, non_arg = _f_variables(Fr, q, dV_eq, dV_dq)
        f = map(lambda x: x(*f_arg),
                symbols('f1:{0}'.format(m + 1)))
        dV_eq = subs(dV_eq, dict(zip(dV_dq[-m:], f)))
        dV_dq = dV_dq[:-m]

    dV_dq_map = solve(dV_eq, dV_dq)
    dV_dq_list = map(lambda x: dV_dq_map[x], dV_dq)

    if vc_map is None:
        #print('Checking ∂/∂qr(∂V/∂qs) = ∂/∂qs(∂V/∂qr) for all r, s '
        #      '= 1, ..., n.')
        dV_eq = _equivalent_derivatives(dV_dq_list, q)
        if dV_eq != [0] * (n*(n - 1)//2):
            rs = [(r, s) for r in range(n) for s in range(r + 1, n)]
            for (r, s), x in zip(rs, dV_eq):
                if trigsimp(expand(x)) != 0:
                    print(('∂/∂q{0}(∂V/∂q{1}) != ∂/∂q{1}(∂V/∂q{0}). ' +
                           'V does NOT exist.').format(r + 1, s + 1))
                    print('∂/∂q{0}(∂V/∂q{1}) = {2}'.format(
                            r + 1, s + 1, dV_dq_list[r].diff(q[s])))
                    print('∂/∂q{1}(∂V/∂q{0}) = {2}'.format(
                            r + 1, s + 1, dV_dq_list[s].diff(q[r])))
                    return None
    else:
        dV_dq_list += f
        # Unable to take diff of 'fm.diff(qs)', replace with dummy symbols.
        dfdq = [Dummy('∂f{0}/∂q{1}'.format(i + 1, j + 1))
                for i in range(len(f)) for j in range(n)]
        dfdq_replace = lambda x: reduce(
                lambda y, z: y.replace(z[0], z[1]) if z[0] != 0 else y,
                zip([fm.diff(qs) for fm in f for qs in q], dfdq),
                x)
        dV_eq = map(dfdq_replace,
                    _equivalent_derivatives(dV_dq_list, q))

        X = Matrix(dfdq)
        Z = Matrix([map(lambda x: diff(dV_eqi, x), dfdq)
                    for dV_eqi in dV_eq])
        if Z.rank() == n * (n - 1) / 2:
            print('ρ == n(n - 1)/2')
            print('V may exist but cannot be found by this procedure.')
            return None

        Y = expand(Z*X - Matrix(dV_eq))
        ZI_rref, _ = Matrix.hstack(Z, Matrix.eye(Z.shape[0])).rref()
        # E is the matrix of elementary row operations that gives rref(Z).
        E = ZI_rref[:, Z.shape[1]:]
        f_eq = (E * Y)[Z.rank():]
        f_map = solve(f_eq, f)
        if sorted(f_map.keys(), cmp=lambda x, y: x.compare(y)) != f:
            print('Unable to solve for all f uniquely.')
            return None
        for k, v in f_map.iteritems():
            for qi in non_arg:
                if v.diff(qi) != 0:
                    print('{0} should not be a function of {1}'.format(k, qi))
                    return None
        dV_dq_list = map(trigsimp, (subs(dV_dq_list, f_map)))

    return function_from_partials(dV_dq_list, q)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
forces = [(pS_star, -M*g*F.x), (pQ, Q1*A.x)] # no friction at point Q
torques = [(A, -TB*A.z), (A, -TC*A.z), (B, TB*A.z), (C, TC*A.z)]
partials = partial_velocities(zip(*forces + torques)[0], [u1, u2],
                              F, kde_map, vc_map, express_frame=A)
Fr, _ = generalized_active_forces(partials, forces + torques)

q = [q1, q2, q3, q4, q5]
u = [u1, u2]
n = len(q)
p = len(u)
m = n - p

if vc_map is not None:
    u += sorted(vc_map.keys(), cmp=lambda x, y: x.compare(y))

dV_dq = symbols('∂V/∂q1:{0}'.format(n + 1))
dV_eq = Matrix(Fr).T

W_sr, _ = kde_matrix(u, kde_map)
if vc_map is not None:
    A_kr, _ = vc_matrix(u, vc_map)
else:
    A_kr = Matrix.zeros(m, p)

for s in range(W_sr.shape[0]):
    dV_eq += dV_dq[s] * (W_sr[s, :p] + W_sr[s, p:]*A_kr[:, :p])

for elem in dV_eq:
    print('{0} = 0'.format(msprint(elem)))

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
pS.set_acc(Q, pS.acc(S) + pQ.acc(S))
pS.set_acc(E, pS.acc(Q) + pE.acc(Q))
pS.set_acc(M, pS.acc(E) + pM.acc(E))

#print('acc in frame\t{0}\t{1}\t{2}\t{3}'.format(*frames))
#for p in points:
#    print('point {0}:\t{1:0.3g}\t{2:0.3g}\t{3:0.3g}\t{4:0.3g}'.format(
#    #print('point {0}:\t{1}\t{2}\t{3}\t{4}'.format(
#            p, *map(lambda x: float(p.acc(x).magnitude().subs(symbol_values)),
#                    frames)))


idx_Q = frames.index(Q)
print('acc in frame ratio\t{0}\t{1}\t{2}'.format(S, E, M))
acc_ratios = Matrix.zeros(4, 3)
for i, p in enumerate(points[2:]):
    acc_values = map(lambda x: p.acc(x).magnitude().subs(symbol_values),
                     frames)
    a_p_Q = acc_values[idx_Q]
    acc_values = [float(x / a_p_Q)
                  for x in acc_values[:idx_Q] + acc_values[idx_Q + 1:]]
    acc_ratios[i, :] = Matrix(acc_values).T
    print('object {0}:\t\t{1:0.3g}\t{2:0.3g}\t{3:0.3g}'.format(
            p, *acc_values))
print('Approximately Newtonian reference frames have a near 1 ratio.')

min_ratio = 0.9 # minimum value if frame is approximately Newtonian.
assert acc_ratios[0, 0] >= min_ratio
assert acc_ratios[0, 1] < min_ratio
assert acc_ratios[0, 2] < min_ratio

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
kde = [u1 - dot(A.ang_vel_in(F), A.x), u2 - dot(pD.vel(F), A.y), u3 - q3d, u4 - q4d, u5 - q5d]
kde_map = solve(kde, [q1d, q2d, q3d, q4d, q5d])
vc = [dot(p.vel(F), A.y) for p in [pB_hat, pC_hat]] + [dot(pD.vel(F), A.z)]
vc_map = solve(subs(vc, kde_map), [u3, u4, u5])

# inertias of bodies A, B, C
# IA22, IA23, IA33 are not specified in the problem statement, but are
# necessary to define an inertia object. Although the values of
# IA22, IA23, IA33 are not known in terms of the variables given in the
# problem statement, they do not appear in the general inertia terms.
inertia_A = inertia(A, IA, IA22, IA33, 0, IA23, 0)
K = mB * R ** 2 / 4
J = mB * R ** 2 / 2
inertia_B = inertia(B, K, K, J)
inertia_C = inertia(C, K, K, J)

# define the rigid bodies A, B, C
rbA = RigidBody("rbA", pA_star, A, mA, (inertia_A, pA_star))
rbB = RigidBody("rbB", pB_star, B, mB, (inertia_B, pB_star))
rbC = RigidBody("rbC", pC_star, C, mB, (inertia_C, pC_star))
bodies = [rbA, rbB, rbC]

system = [i.masscenter for i in bodies] + [i.frame for i in bodies]
partials = partial_velocities(system, [u1, u2], F, kde_map, vc_map)

M = trigsimp(inertia_coefficient_matrix(bodies, partials))
print("inertia_coefficient_matrix = {0}".format(msprint(M)))

M_expected = Matrix([[IA + mA * a ** 2 + mB * (R ** 2 / 2 + 3 * b ** 2), 0], [0, mA + 3 * mB]])
assert expand(M - M_expected) == Matrix.zeros(2)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import zeros [as 别名]
)
Y2 = (
    Z38
    + Z43
    + Z52
    + mB * LB * Z22
    + mC * Z9 * Z26
    + mD * (Z14 * Z32 - p1 * Z34)
    - T2_A_B
    + g * ((mB * LB + mC * Z9 + mD * Z14) * cos(q1) - mD * p1 * sin(q1))
)
Y3 = mC * Z28 + mD * Z34 - K3_B_C + g * (mC + mD) * sin(q1)
X_rs_expected = Matrix([[X11, X12, X13], [X21, X22, X23], [X31, X32, X33]])
Y_s_expected = Matrix([Y1, Y2, Y3])

assert trigsimp(expand(X_rs - X_rs_expected)) == Matrix.zeros(3)
assert trigsimp(expand(Y_s - Y_s_expected)) == Matrix.zeros(3, 1)

# 11.8
# If D* lies on line B*C*, then p1 = 0, p2 = 0.
# If each central principal axis of D is parallel to one of b1, b2, b3,
# then D12, D23, D31 = 0.
D_prop_map = {p1: 0, p2: 0, D12: 0, D23: 0, D31: 0}
DE_2 = Matrix(dyn_eq).subs(D_prop_map)
X_rs_2 = Matrix(map(lambda x: DE_2.T.diff(x), ud)).T
Y_s_2 = -expand(DE_2 - X_rs_2 * Matrix(ud))

# Solution will have the form ud_r = Y_r/X_rr (r = 1, 2, 3)
# if the X_rs matrix is diagonal.
X_rs_expected_2 = Matrix([[X11, 0, 0], [0, X22, 0], [0, 0, X33]]).subs(D_prop_map)
Y_s_expected_2 = Y_s_expected.subs(D_prop_map)