欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

线性代数基础(python3)

程序员文章站 2022-07-12 14:01:17
...
from math import sqrt,acos,pi
from decimal import Decimal,getcontext
getcontext().prec =30 #保留30小数

class Vector(object):
    def __init__(self, coordinates):
        try:
            if not coordinates:
                raise ValueError
            self.coordinates = tuple([Decimal(x) for x in coordinates])
            self.dimension = len(self.coordinates)
        except ValueError:
            raise ValueError('The coordinates must be nonempty')

        except TypeError:
            raise TypeError('The coordinates must be an iterable')
    CANNOT_NORMALIZE_ZERO_VECTOR_MSG = 'Cannot normalize the zero vector'

    #向量相加
    def plus(self,v):
        new_coordinates = [x+y for x,y in zip(self.coordinates,v.coordinates)]
        return Vector(new_coordinates)

    #向量相减
    def minus(self,v):
        new_coordinates = [x-y for x,y in zip(self.coordinates,v.coordinates)]
        return Vector(new_coordinates)

    #数乘
    def times_scalar(self,c):
        new_coordinates = [Decimal(c)*x for x in self.coordinates]
        return Vector(new_coordinates)

    #向量的模
    def magnitude(self):
        coordinates_squared = [x*x for x in self.coordinates]
        return sqrt(sum(coordinates_squared))

    #方向函数(化为单位向量)
    def normalized(self):
        try:
            magnitude = self.magnitude()
            return self.times_scalar(Decimal(1.0)/Decimal(magnitude))
        except ZeroDivisionError:
            raise Exception(self.CANNOT_NORMALIZE_ZERO_VECTOR_MSG)

    #向量点乘(内积)
    def dot(self,v):
        return sum([x*y for x,y in zip(self.coordinates,v.coordinates)])
    #向量叉乘(外积、向量积)
    def cross(self,v):
        try:
            x_1,y_1,z_1=self.coordinates
            x_2,y_2,z_2=v.coordinates
            new_coordinates=[y_1*z_2-y_2*z_1,
                             -(x_1*z_2-x_2*z_1),
                             x_1*y_2-x_2*y_1]
            return Vector(new_coordinates)
        except ValueError as e:
            msg = str(e)
            if msg == 'need more than 2 values to unpack':
                self_embedded_in_R3 = Vector(self.coordinates+('0',))
                v_embedded_in_R3 = Vector(v.coordinates+('0',))
                return self_embedded_in_R3*v_embedded_in_R3
            elif (msg == 'too many values to unpack' or
                  msg == 'need more than 1 values to unpack'):
                  raise Exception(self.ONLY_DEFINED_IN_TWO_THREE_DIMS_MSG)
            else:
                raise e
            
    #用于计算两个向量形成的平行四边形的面积
    def area_of_parallelogram_with(self,v):
        cross_product = self.cross(v)
        return cross_product.magnitude()
    
    #用于计算两个向量形成的三角形的面积
    def area_of_triangle_with(self,v):
        return self.area_of_parallelogram_with(v)/2.0
    
    #计算向量夹角
    
    def angle_with(self,v,in_degrees=False):
        try:
            u1=self.normalized()
            u2=v.normalized()
            dots=u1.dot(u2)
            if abs(abs(dots) - 1) < 1e-10:
                if dots < 0:
                    dots = -1
                else:
                    dots = 1
            angle_in_radians =acos(dots)

            if in_degrees:
                degrees_per_radian =180./pi
                return (angle_in_radians*degrees_per_radian)
            else:
                return(angle_in_radians)
        
        except Exception as e:
            if str(e)==self.CANNOT_NORMALIZE_ZERO_VECTOR_MSG:
                raise Exception ("Cannot compute an angle with the zero vector")
            else:
                raise e
            
    #判断两向量之间平行
    def is_parallel_to(self,v):
        return (self.is_zero() or
                v.is_zero() or
                self.angle_with(v)==0 or
                self.angle_with(v)==pi)
    
    def is_zero(self,tolerance=1e-10):#1乘以10的-10次方
        return self.magnitude()<tolerance

    #判断两向量正交
    def is_orthogonal_to(self,v,tolerance=1e-10):
        return abs(self.dot(v))<tolerance

    #利用向量投影分解的水平向量
    def component_parallel_to(self,basis):
        try:
            u=basis.normalized()
            weight=self.dot(u)
            return u.times_scalar(weight)
        except Exception as e:
            if str(e)==self.CANNOT_NORMALIZE_ZERO_VECTOR_MSG:
                raise Exception(self.NO_UNIQUE_PARALLEL_COMPONENT_MSG)
            else:
                raise e

    #利用向量投影分解出的垂直向量
    def component_orthogonal_to(self,basis):
        try:
            projection=self.component_parallel_to(basis)
            return self.minus(projection)
        except Exception as e:
            if str(e)==self.NO_UNIQUE_PARALLEL_COMPONENT_MSG:
                raise Exception(self.NO_UNIQUE_ORTHOGNOAL_COMPONENT_MSG)
            else:
                raise e

    #
    
    
    def __str__(self):
        return 'Vector: {}'.format(self.coordinates)

    def __eq__(self, v):
        return self.coordinates == v.coordinates
#向量相加
vector1 = Vector(['8.218','-9.341'])
vector2 = Vector(['-1.129','2.111'])
print (vector1.plus(vector2))#Vector.plus(vector1,vector2)


#向量相减
vector3 = Vector(['7.119','8.215'])
vector4 = Vector(['-8.223','0.878'])
print (vector3.minus(vector4))#Vector.minus(vector3,vector4)

#数乘
vector5 = Vector(['1.671','-1.012','-0.318'])
c=7.41
print (vector5.times_scalar(c))#Vector.times_scalar(vector5,c)

#向量的模和方向函数(化为单位向量)
vector6 = Vector(['3','4'])
print(Vector.magnitude(vector6))
print(Vector.normalized(vector6))


#向量点乘(内积)
v = Vector(['7.887','4.138'])
w = Vector(['-8.802','6.776'])
print(v.dot(w))

v = Vector(['-5.955','-4.904','-1.874'])
w = Vector(['-4.496','-8.755','7.103'])
print(v.dot(w))

#计算向量夹角
v = Vector(['3.183','-7.627'])
w =Vector(['-2.668','5.319'])
print(v.angle_with(w))

v = Vector(['7.6','1','5.188'])
w = Vector(['2.751','8.259','3.985'])
print(v.angle_with(w,in_degrees=True))

#判断两向量之间平行或正交
print('first pair')
v = Vector(['-7.579','-7.88'])
w = Vector(['22.737','23.64'])
print ('is parallel:',v.is_parallel_to(w))
print ('is orthogonal:',v.is_orthogonal_to(w))

print('second pair')
v = Vector(['-2.029','9.97','4.172'])
w = Vector(['-9.231','-6.639','-7.245'])
print ('is parallel:',v.is_parallel_to(w))
print ('is orthogonal:',v.is_orthogonal_to(w))

print('third pair')
v = Vector(['-2.328','-7.284','-1.214'])
w = Vector(['-1.1821','1.072','-2.94'])
print ('is parallel:',v.is_parallel_to(w))
print ('is orthogonal:',v.is_orthogonal_to(w))

print('four pair')
v = Vector(['2.118','4.827'])
w = Vector(['0','0'])
print ('is parallel:',v.is_parallel_to(w))
print ('is orthogonal:',v.is_orthogonal_to(w))

#利用向量投影分解出一个水平向量、一个垂直向量
print('\n#1')
v = Vector(['3.039','1.879'])
w = Vector(['0.825','2.036'])
print(v.component_parallel_to(w))

print('\n#2')
v = Vector(['-9.88','-3.264','-8.159'])
w = Vector(['-2.155','-9.353','-9.473'])
print(v.component_orthogonal_to(w))

print('\n#3')
v = Vector(['3.009','-6.172','3.692','-2.51'])
w = Vector(['6.404','-9.144','2.759','8.718'])
vpar = v.component_parallel_to(w)
vort = v.component_orthogonal_to(w)
print('parallel component:',vpar)
print('orthogonal component:',vort)

未完待续