Python Matrix.dot方法代码示例

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import dot [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 dot [as 别名]
def test_issue_15129_trigsimp_methods():
    t1 = Matrix([sin(Rational(1, 50)), cos(Rational(1, 50)), 0])
    t2 = Matrix([sin(Rational(1, 25)), cos(Rational(1, 25)), 0])
    t3 = Matrix([cos(Rational(1, 25)), sin(Rational(1, 25)), 0])
    r1 = t1.dot(t2)
    r2 = t1.dot(t3)
    assert trigsimp(r1) == cos(S(1)/50)
    assert trigsimp(r2) == sin(S(3)/50)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import dot [as 别名]
def get_interblock_dists(grid, direction, opposite=False):
    """Computes a list of indices of distributions that would be transferred
    to a node pointed to by the vector 'direction'.
    """
    d = Matrix((direction,))

    def process_dists(dists):
        if opposite:
            return [grid.idx_opposite[x] for x in dists]
        else:
            return dists

    ret = []
    for i, ei in enumerate(grid.basis):
        if ei.dot(d) >= d.dot(d):
            ret.append(i)
    return process_dists(ret)

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

#.........这里部分代码省略.........
    >>> refraction_angle(r1, 1, 1, n)
    Matrix([
    [ 1],
    [ 1],
    [-1]])
    >>> refraction_angle(r1, 1, 1, plane=P)
    Ray3D(Point3D(0, 0, 0), Point3D(1, 1, -1))

    With different index of refraction of the two media

    >>> n1, n2 = symbols('n1, n2')
    >>> refraction_angle(r1, n1, n2, n)
    Matrix([
    [                                n1/n2],
    [                                n1/n2],
    [-sqrt(3)*sqrt(-2*n1**2/(3*n2**2) + 1)]])
    >>> refraction_angle(r1, n1, n2, plane=P)
    Ray3D(Point3D(0, 0, 0), Point3D(n1/n2, n1/n2, -sqrt(3)*sqrt(-2*n1**2/(3*n2**2) + 1)))

    """
    # A flag to check whether to return Ray3D or not
    return_ray = False

    if plane is not None and normal is not None:
        raise ValueError("Either plane or normal is acceptable.")

    if not isinstance(incident, Matrix):
        if is_sequence(incident):
            _incident = Matrix(incident)
        elif isinstance(incident, Ray3D):
            _incident = Matrix(incident.direction_ratio)
        else:
            raise TypeError(
                "incident should be a Matrix, Ray3D, or sequence")
    else:
        _incident = incident

    # If plane is provided, get direction ratios of the normal
    # to the plane from the plane else go with `normal` param.
    if plane is not None:
        if not isinstance(plane, Plane):
            raise TypeError("plane should be an instance of geometry.plane.Plane")
        # If we have the plane, we can get the intersection
        # point of incident ray and the plane and thus return
        # an instance of Ray3D.
        if isinstance(incident, Ray3D):
            return_ray = True
            intersection_pt = plane.intersection(incident)[0]
        _normal = Matrix(plane.normal_vector)
    else:
        if not isinstance(normal, Matrix):
            if is_sequence(normal):
                _normal = Matrix(normal)
            elif isinstance(normal, Ray3D):
                _normal = Matrix(normal.direction_ratio)
                if isinstance(incident, Ray3D):
                    intersection_pt = intersection(incident, normal)
                    if len(intersection_pt) == 0:
                        raise ValueError(
                            "Normal isn't concurrent with the incident ray.")
                    else:
                        return_ray = True
                        intersection_pt = intersection_pt[0]
            else:
                raise TypeError(
                    "Normal should be a Matrix, Ray3D, or sequence")
        else:
            _normal = normal

    n1, n2 = None, None

    if isinstance(medium1, Medium):
        n1 = medium1.refractive_index
    else:
        n1 = sympify(medium1)

    if isinstance(medium2, Medium):
        n2 = medium2.refractive_index
    else:
        n2 = sympify(medium2)

    eta = n1/n2  # Relative index of refraction
    # Calculating magnitude of the vectors
    mag_incident = sqrt(sum([i**2 for i in _incident]))
    mag_normal = sqrt(sum([i**2 for i in _normal]))
    # Converting vectors to unit vectors by dividing
    # them with their magnitudes
    _incident /= mag_incident
    _normal /= mag_normal
    c1 = -_incident.dot(_normal)  # cos(angle_of_incidence)
    cs2 = 1 - eta**2*(1 - c1**2)  # cos(angle_of_refraction)**2
    if cs2.is_negative:  # This is the case of total internal reflection(TIR).
        return 0
    drs = eta*_incident + (eta*c1 - sqrt(cs2))*_normal
    # Multiplying unit vector by its magnitude
    drs = drs*mag_incident
    if not return_ray:
        return drs
    else:
        return Ray3D(intersection_pt, direction_ratio=drs)

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import dot [as 别名]
def deviation(incident, medium1, medium2, normal=None, plane=None):
    """
    This function calculates the angle of deviation of a ray
    due to refraction at planar surface.

    Parameters
    ==========

    incident : Matrix, Ray3D, or sequence
        Incident vector
    medium1 : sympy.physics.optics.medium.Medium or sympifiable
        Medium 1 or its refractive index
    medium2 : sympy.physics.optics.medium.Medium or sympifiable
        Medium 2 or its refractive index
    normal : Matrix, Ray3D, or sequence
        Normal vector
    plane : Plane
        Plane of separation of the two media.

    Examples
    ========

    >>> from sympy.physics.optics import deviation
    >>> from sympy.geometry import Point3D, Ray3D, Plane
    >>> from sympy.matrices import Matrix
    >>> from sympy import symbols
    >>> n1, n2 = symbols('n1, n2')
    >>> n = Matrix([0, 0, 1])
    >>> P = Plane(Point3D(0, 0, 0), normal_vector=[0, 0, 1])
    >>> r1 = Ray3D(Point3D(-1, -1, 1), Point3D(0, 0, 0))
    >>> deviation(r1, 1, 1, n)
    0
    >>> deviation(r1, n1, n2, plane=P)
    -acos(-sqrt(-2*n1**2/(3*n2**2) + 1)) + acos(-sqrt(3)/3)

    """
    refracted = refraction_angle(incident,
                                 medium1,
                                 medium2,
                                 normal=normal,
                                 plane=plane)
    if refracted != 0:
        if isinstance(refracted, Ray3D):
            refracted = Matrix(refracted.direction_ratio)

        if not isinstance(incident, Matrix):
            if is_sequence(incident):
                _incident = Matrix(incident)
            elif isinstance(incident, Ray3D):
                _incident = Matrix(incident.direction_ratio)
            else:
                raise TypeError(
                    "incident should be a Matrix, Ray3D, or sequence")
        else:
            _incident = incident

        if plane is None:
            if not isinstance(normal, Matrix):
                if is_sequence(normal):
                    _normal = Matrix(normal)
                elif isinstance(normal, Ray3D):
                    _normal = Matrix(normal.direction_ratio)
                else:
                    raise TypeError(
                        "normal should be a Matrix, Ray3D, or sequence")
            else:
                _normal = normal
        else:
            _normal = Matrix(plane.normal_vector)

        mag_incident = sqrt(sum([i**2 for i in _incident]))
        mag_normal = sqrt(sum([i**2 for i in _normal]))
        mag_refracted = sqrt(sum([i**2 for i in refracted]))
        _incident /= mag_incident
        _normal /= mag_normal
        refracted /= mag_refracted
        i = acos(_incident.dot(_normal))
        r = acos(refracted.dot(_normal))
        return i - r

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import dot [as 别名]
sigma1 = 0.53
sigma2 = 0.02
sigma12 = 0.3
P[0] = np.array([[sigma1, sigma12], [sigma12, sigma2]])
print("A posteriori initial error covariance: P_0(+) =\n{0}".format(P[0]))

for k in range(1, N):

    # Extrapolation to next time step
    x_est_k_pre = x[k - 1]  # state transition is identity
    Pk_pre = P[k - 1]  # no process noise

    # Measurement update
    zk = x[k, 1] + v[k]
    # Kalman gain
    Kk = Pk_pre.dot(H.transpose()) / ((H.dot(Pk_pre)).dot(H.transpose()) + r2)

    # A posteriori state estimate update
    x_est[k] = x_est_k_pre + Kk.dot(zk - H.dot(x_est_k_pre))

    # A posteriori covariance update
    P[k, 0, 0] = P[k - 1, 0, 0] - P[k - 1, 0, 1] ** 2.0 / (r2 + P[k - 1, 1, 1])
    P[k, 0, 1] = P[k - 1, 0, 1] - P[k - 1, 0, 1] * P[k - 1, 1, 1] / (r2 + P[k - 1, 1, 1])
    P[k, 1, 0] = P[k - 1, 0, 1] - P[k - 1, 0, 1] * P[k - 1, 1, 1] / (r2 + P[k - 1, 1, 1])
    P[k, 1, 1] = P[k - 1, 1, 1] - P[k - 1, 1, 1] ** 2.0 / (r2 + P[k - 1, 1, 1])

print("Final error covariance:\n{0}".format(P[-1]))


import matplotlib.pyplot as plt

# 需要导入模块: from sympy import Matrix [as 别名]
# 或者: from sympy.Matrix import dot [as 别名]
 def quad_casimir(self, highest):
     highest_weight = Matrix([highest]) * self.fund_weights
     return highest_weight.norm()**2 + 2 * highest_weight.dot(self.delta)