Python Matrix.det方法代码示例

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
    def cal(self):
        global X,Y,Z,Gauss_Curvature

        xu=Matrix([diff(f,u) for f in self.S])
        xv=Matrix([diff(f,v) for f in self.S])
        xuu=Matrix([diff(f,u) for f in xu])
        xvv=Matrix([diff(f,v) for f in xv])
        xuv=Matrix([diff(f,v) for f in xu])
        E=simplify(xu.dot(xu))
        G=simplify(xv.dot(xv))
        F=simplify(xu.dot(xv))


        H1=Matrix([xuu,xu,xv]).reshape(3,3)
        H2=Matrix([xvv,xu,xv]).reshape(3,3)
        H3=Matrix([xuv,xu,xv]).reshape(3,3)
        K=simplify(((H1.det()*H2.det() -(H3.det()**2)))/ (xu.norm()**2*xv.norm()**2 -F*F)**2)

        #Pass to numpy
        du=float(self.umax-self.umin)/100
        dv=float(self.vmax-self.vmin)/100
        [U,V] = np.mgrid[self.umin:self.umax+du:du,self.vmin:self.vmax+dv:dv]
        # convert Sympy formula to numpy lambda functions
        F1=lambdify((u,v), self.S[0], "numpy")
        F2=lambdify((u,v), self.S[1], "numpy")
        F3=lambdify((u,v), self.S[2], "numpy")
        F4=lambdify((u,v), K, "numpy")
        #Calculate numpy arrays 
        self.X=F1(U,V)
        self.Y=F2(U,V)
        self.Z=F3(U,V)
        self.Gauss_Curvature=F4(U,V)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
    def get_dixon_polynomial(self):
        r"""
        Returns
        =======

        dixon_polynomial: polynomial
            Dixon's polynomial is calculated as:

            delta = Delta(A) / ((x_1 - a_1) ... (x_n - a_n)) where,

            A =  |p_1(x_1,... x_n), ..., p_n(x_1,... x_n)|
                 |p_1(a_1,... x_n), ..., p_n(a_1,... x_n)|
                 |...             , ...,              ...|
                 |p_1(a_1,... a_n), ..., p_n(a_1,... a_n)|
        """
        if self.m != (self.n + 1):
            raise ValueError('Method invalid for given combination.')

        # First row
        rows = [self.polynomials]

        temp = list(self.variables)

        for idx in range(self.n):
            temp[idx] = self.dummy_variables[idx]
            substitution = {var: t for var, t in zip(self.variables, temp)}
            rows.append([f.subs(substitution) for f in self.polynomials])

        A = Matrix(rows)

        terms = zip(self.variables, self.dummy_variables)
        product_of_differences = Mul(*[a - b for a, b in terms])
        dixon_polynomial = (A.det() / product_of_differences).factor()

        return poly_from_expr(dixon_polynomial, self.dummy_variables)[0]

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
def Ixx(sections):
    r"""
    from http://en.wikipedia.org/wiki/Second_moment_of_area
    \f[ J_{xx} = \frac{1}{12} \sum_{i = 1}^{n-1} ( y_i^2 + y_i y_{i+1} + y_{i+1}^2 ) a_i \f]
    \f[ J_{yy} = \frac{1}{12} \sum_{i = 1}^{n-1} ( x_i^2 + x_i x_{i+1} + x_{i+1}^2 ) a_i \f]
    \f[ J_{xy} = \frac{1}{24} \sum_{i = 1}^{n-1} ( x_i y_{i+1} + 2 x_i y_i + 2 x_{i+1} y_{i+1} + x_{i+1} y_i ) a_i \f]
    """
    h = Symbol('h')
    b = Symbol('b')
    Ixx = 0
    Iyy = 0
    Ixy = 0
    half = Rational(1, 2)
    CG = [0, 0]
    Area = 0
    o3 = Rational(1, 3)
    for i, section in enumerate(sections):
        (sign, p1, p2) = section
        #print "p%s = %s" %(i,p1)
        #print "p%s = %s" %(i+1,p2)
        #Atri = half * abs(p1[0]*p2[1])
        AA = Matrix([p1, p2])
        #print "AA = \n",AA
        detAA = AA.det()
        msgAA = '%s' % (detAA)
        #print "msgAA = ",msgAA
        if '-' in msgAA:
            detAA *= Rational(-1, 1)
        Atri = detAA
        #Atri = half*b*h
        #Atri = (p1[0]*p2[1] - p2[0]*p1[1])*half
        ixxi = p1[1] * p1[1] + p1[1] * p2[1] + p2[1] * p2[1]
        iyyi = p1[0] * p1[0] + p1[0] * p2[0] + p2[0] * p2[0]
        ixyi = (p1[0] * p1[1] + 2 * p1[0] * p1[1] + 2 * p2[0] * p2[1] + p2[0] * p1[1])
        Ixxi = ixxi * Atri
        Iyyi = iyyi * Atri
        Ixyi = ixyi * Atri
        #print "Atri = ",Atri
        #print "Ixxi = ",Ixxi
        Ixx += Ixxi
        Iyy += Iyyi
        Ixy += Ixyi
        #print "Ixx = ",simplify(Ixx),'\n'
        CG = [CG[0] + p1[0] * o3 + p2[0] * o3, CG[1] + p1[1] * o3 + p2[1] * o3]
        Area += Atri
    CG = [CG[0] / Area, CG[1] / Area]
    #print "Ixx = ",simplify(Ixx)
    o12 = Rational(1, 12)
    o24 = Rational(1, 24)
    Ixx = Ixx * o12
    Iyy = Iyy * o12
    Ixy = Ixy * o24
    J = Ixx + Iyy
    print "Area = ", simplify(Area * half)
    print "CG   = ", simplify(CG)
    print "Ixx  = ", simplify(Ixx)
    print "Iyy  = ", simplify(Iyy)
    print "Ixy  = ", simplify(Ixy)
    print "J    = ", simplify(J), '\n'

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
 def Delta(self, idx=-1):
     """
     Construct the matrix corresponding to `idx`'th point, if `idx>0`
     Otherwise returns the discriminant.
     """
     if self.Env == 'sympy':
         from sympy import Matrix, Subs
     elif self.Env == 'symengine':
         from symengine import Matrix, Subs
     M = []
     for j in range(self.num_points):
         pnt = self.Points[j]
         row = []
         for mon in self.MonoBase:
             if j != idx:
                 chng = {self.Vars[i]: pnt[i] for i in range(self.num_vars)}
                 row.append(mon.subs(chng))
             else:
                 row.append(mon)
         M.append(row)
     Mtx = Matrix(M)
     return Mtx.det()

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
def test_determinant():
    x, y, z = Symbol('x'), Symbol('y'), Symbol('z')

    M = Matrix((1,))

    assert M.det(method="bareis") == 1
    assert M.det(method="berkowitz") == 1

    M = Matrix(( (-3,  2),
                 ( 8, -5) ))

    assert M.det(method="bareis") == -1
    assert M.det(method="berkowitz") == -1

    M = Matrix(( (x,   1),
                 (y, 2*y) ))

    assert M.det(method="bareis") == 2*x*y-y
    assert M.det(method="berkowitz") == 2*x*y-y

    M = Matrix(( (1, 1, 1),
                 (1, 2, 3),
                 (1, 3, 6) ))

    assert M.det(method="bareis") == 1
    assert M.det(method="berkowitz") == 1

    M = Matrix(( ( 3, -2,  0, 5),
                 (-2,  1, -2, 2),
                 ( 0, -2,  5, 0),
                 ( 5,  0,  3, 4) ))

    assert M.det(method="bareis") == -289
    assert M.det(method="berkowitz") == -289

    M = Matrix(( ( 1,  2,  3,  4),
                 ( 5,  6,  7,  8),
                 ( 9, 10, 11, 12),
                 (13, 14, 15, 16) ))

    assert M.det(method="bareis") == 0
    assert M.det(method="berkowitz") == 0

    M = Matrix(( (3, 2, 0, 0, 0),
                 (0, 3, 2, 0, 0),
                 (0, 0, 3, 2, 0),
                 (0, 0, 0, 3, 2),
                 (2, 0, 0, 0, 3) ))

    assert M.det(method="bareis") == 275
    assert M.det(method="berkowitz") == 275

    M = Matrix(( (1, 0,  1,  2, 12),
                 (2, 0,  1,  1,  4),
                 (2, 1,  1, -1,  3),
                 (3, 2, -1,  1,  8),
                 (1, 1,  1,  0,  6) ))

    assert M.det(method="bareis") == -55
    assert M.det(method="berkowitz") == -55

    M = Matrix(( (-5,  2,  3,  4,  5),
                 ( 1, -4,  3,  4,  5),
                 ( 1,  2, -3,  4,  5),
                 ( 1,  2,  3, -2,  5),
                 ( 1,  2,  3,  4, -1) ))

    assert M.det(method="bareis") == 11664
    assert M.det(method="berkowitz") == 11664

    M = Matrix(( ( 2,  7, -1, 3, 2),
                 ( 0,  0,  1, 0, 1),
                 (-2,  0,  7, 0, 2),
                 (-3, -2,  4, 5, 3),
                 ( 1,  0,  0, 0, 1) ))

    assert M.det(method="bareis") == 123
    assert M.det(method="berkowitz") == 123

    M = Matrix(( (x,y,z),
                 (1,0,0),
                 (y,z,x) ))

    assert M.det(method="bareis") == z**2 - x*y
    assert M.det(method="berkowitz") == z**2 - x*y

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
the mass center of B.
"""

A = ReferenceFrame('A')

# part a
print("a) mass center of surface area")
SA = (2 * pi * R) * 2
rho_a = m / SA

B = A.orientnew('B', 'Axis', [theta, A.x])
p = x * A.x  + r * B.y
y_ = dot(p, A.y)
z_ = dot(p, A.z)
J = Matrix([x, y_, z_]).jacobian([x, r, theta])
dJ = simplify(J.det())
print("dJ = {0}".format(dJ))

# calculate center of mass for the cut cylinder
# ranges from x = [1, 3], y = [-1, 1], z = [-1, 1] in Fig. P5.1
mass_cc_a = rho_a * 2*pi*R
cm_cc_a = (integrate_v(integrate_v((rho_a * p * dJ).subs(r, R),
                                           A, (theta, acos(1-x),
                                               2*pi - acos(1-x))),
                               A, (x, 0, 2)) / mass_cc_a +
           A.x)
print("cm = {0}".format(cm_cc_a))

mass_cyl_a = rho_a * 2*pi*R
cm_cyl_a = A.x/S(2)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
# -*- coding: utf-8 -*-
"""
Spyder Editor

This is a temporary script file.
"""

from sympy import symbols, Matrix, simplify

a10, a20, a30, a11, a21, a31, a12, a22, a32, a42, a03, a23, a43, a04, a44 = symbols("a10 a20 a30 a11 a21 a31 a12 a22 a32 a42 a03 a23 a43 a04 a44")

M = Matrix([[0, a10, a20, a30, 0], [0, a11, a21, a31, 0], [0, a12, a22, a32, a42], [a03, 0, a23, 0, a43], [a04, 0, 0, 0, a44]])

M_det = M.det()
M_inv = M.inverse_ADJ()

print(M_det)

print(simplify(M_inv * M_det))

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
    def __init__(self, domain, mapping, orientation='+'):
        '''
        Construct the set. Orientation is used to construct normal if relevant.
        Positive for 2d-3d means normal = cross(d(mapping)/ds, d(mapping)/dt)
        Positve for 1d-2d means normal = rotate tangent counter-clockwie
        '''
        assert isinstance(domain, ParameterDomain)
        # Topological dimension
        self._tdim = len(domain) 

        # Map is a an iterable with length = gdim that specifies [x, y, z] as
        # functions of parameters of domain
        if not hasattr(mapping, '__len__'):
            mapping = (mapping, )

        # Check that the component uses known params or numbers
        for comp in mapping:
            comp_extras = comp.atoms() - domain.parameters
            
            params_only = len(comp_extras) == 0 
            extras_are_numbers = all(isinstance(item, (float, int, Number))
                                     for item in comp_extras)
            
            assert params_only or extras_are_numbers, 'Invalid mapping'

        # Geometrical dimension
        self._gdim = len(mapping)

        # I am only interested in at most 3d gdim
        assert 0 < self._gdim < 4, 'Geometrical dimension %s not supported' % self._gdim
        assert 0 < self._tdim < 4, 'Topological dimension %s not supported' % self._tdim
        assert self._tdim <= self._gdim, 'Topolog. dim > geometric. dim not allowed' 

        # We rememeber mapping as dictionary x->mapping[0], y->mapping[1], ...
        xyz = symbols('x, y, z')
        self._mapping = dict((var, comp) for var, comp in zip(xyz, mapping))
        # Remeber the domain
        self._pdomain = domain

        params = domain.parameters
        # Every mapping has a Jacobian but not every has normal and tangent
        self._n = None
        self._tau = None
       
        # Volumes, Square matrix only Jacobian 
        if self._tdim == self._gdim:
            Jac = Matrix([[diff(comp, var) for var in params] for comp in mapping])
            self._J = abs(Jac.det())
        # Curves and surfaces have normals or tangents in addition to Jacobian
        else:
            # Curves
            if self._tdim == 1:
                # Tagent
                self._tau = Vector([diff(comp, list(params)[0]) for comp in mapping])
                # Jacobian is length of tangent
                self._J = sqrt(sum(v**2 for v in self._tau))
                
                # And in 2d we can define a normal
                if self._gdim == 2:
                    R = Tensor([[0, -1], [1, 0]])
                    R = R if orientation == '+' else -R
                    self._n = dot(R, self._tau)

            # Surface in 3d has normal
            elif self._tdim == 2 and self._gdim == 3:
                u0 = Vector([diff(comp, list(params)[0]) for comp in mapping])
                u1 = Vector([diff(comp, list(params)[1]) for comp in mapping])

                self._n = cross(u0, u1)
                self._n = self._n if orientation == '+' else -self._n
                self._J = sqrt(sum(v**2 for v in self._n))

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
rho = m / V

A = ReferenceFrame('A')                    # Torus fixed, x-y in symmetry plane
B = A.orientnew('B', 'Axis', [phi, A.z])   # Intermediate frame
C = B.orientnew('C', 'Axis', [-theta, B.y])# Intermediate frame

# Position vector from torus center to arbitrary point of torus
# R : torus major radius
# s : distance >= 0 from center of torus cross section to point in torus
p = R * B.x + s * C.x

# Determinant of the Jacobian of the mapping from a, b, c to x, y, z
# See Wikipedia for a lucid explanation of why we must comput this Jacobian:
# http://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant#Further_examples
J = Matrix([dot(p, uv) for uv in A]).transpose().jacobian([phi, theta, s])
dv = J.det().trigsimp()      # Need to ensure this is positive
print("dx*dy*dz = {0}*dphi*dtheta*ds".format(dv))

# We want to compute the inertia scalars of the torus relative to it's mass
# center, for the following six unit vector pairs
unit_vector_pairs = [(A.x, A.x), (A.y, A.y), (A.z, A.z),
                     (A.x, A.y), (A.y, A.z), (A.x, A.z)]

# Calculate the six unique inertia scalars using equation 3.3.9 of Kane &
# Levinson, 1985.
inertia_scalars = []
for n_a, n_b in unit_vector_pairs:
    # Integrand of Equation 3.3.9
    integrand = rho * dot(cross(p, n_a), cross(p, n_b)) * dv

    # Compute the integral by integrating over the whole volume of the tours

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
## --- define points D, S*, Q on frame A and their velocities ---
pP_star = Point('P*')
pP_star.set_vel(P, 0)
pP_star.set_vel(C, u2*C.x + u3*C.y)

pQ = pP_star.locatenew('Q', x*C.x + y*C.y + z*C.z)
pQ.set_vel(P, 0)
pQ.v2pt_theory(pP_star, C, P)

## --- map from cartesian to cylindrical coordinates ---
coord_pairs = [(x, x), (y, r*cos(theta)), (z, r*sin(theta))]
coord_map = dict([(x, x),
                  (y, r*cos(theta)),
                  (z, r*sin(theta))])
J = Matrix([coord_map.values()]).jacobian([x, theta, r])
dJ = trigsimp(J.det())

## --- define contact/distance forces ---
# force for a point on ring R1, R2, R3
n = alpha + beta*cos(theta/2) # contact pressure
t = u_prime*n # kinetic friction
tau = -pQ.vel(C).subs(coord_map).normalize() # direction of friction
v = -P.y # direction of surface
point_force = sum(simplify(dot(n*v + t*tau, b)) * b for b in P)

# want to find gen. active forces for motions where u3 = 0
forces = [(pP_star, E*C.x + M*g*C.y),
          (pQ, subs(point_force, u3, 0),
           lambda i: integrate(i.subs(coord_map) * dJ,
                               (theta, -pi, pi)).subs(r, R))]
# 3 rings so repeat the last element twice more

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import det [as 别名]
import numpy as np
from sympy import simplify, Matrix, var

xi = var('xi')
eta = var('eta')

print('Chebyshev Polynomials of the Second Kind')
u = [1, 2*xi]
for i in range(2, 15):
    ui = simplify(2*xi*u[i-1] - u[i-2])
    u.append(ui)
for i, ui in enumerate(u):
    print('u({0}) = {1}'.format(i, ui))

with open('cheb.txt', 'w') as f:
    f.write('Number of terms: {0}\n\n'.format(len(u)))
    f.write(','.join(map(str, u)).replace('**','^') + '\n\n')
    f.write(','.join(map(str, u)).replace('xi', 'eta').replace('**','^'))

print np.outer(u, [ui.subs({xi:eta}) for ui in u])

if False:
    m = Matrix([[2*xi, 1, 0, 0],
                [1, 2*xi, 1, 0],
                [0, 1, 2*xi, 1],
                [0, 0, 1, 2*xi]])

    print('m.det() {0}'.format(simplify(m.det())))