vector.py 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623
  1. from __future__ import annotations
  2. from itertools import product
  3. from sympy.core.add import Add
  4. from sympy.core.assumptions import StdFactKB
  5. from sympy.core.expr import AtomicExpr, Expr
  6. from sympy.core.power import Pow
  7. from sympy.core.singleton import S
  8. from sympy.core.sorting import default_sort_key
  9. from sympy.core.sympify import sympify
  10. from sympy.functions.elementary.miscellaneous import sqrt
  11. from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
  12. from sympy.vector.basisdependent import (BasisDependentZero,
  13. BasisDependent, BasisDependentMul, BasisDependentAdd)
  14. from sympy.vector.coordsysrect import CoordSys3D
  15. from sympy.vector.dyadic import Dyadic, BaseDyadic, DyadicAdd
  16. class Vector(BasisDependent):
  17. """
  18. Super class for all Vector classes.
  19. Ideally, neither this class nor any of its subclasses should be
  20. instantiated by the user.
  21. """
  22. is_scalar = False
  23. is_Vector = True
  24. _op_priority = 12.0
  25. _expr_type: type[Vector]
  26. _mul_func: type[Vector]
  27. _add_func: type[Vector]
  28. _zero_func: type[Vector]
  29. _base_func: type[Vector]
  30. zero: VectorZero
  31. @property
  32. def components(self):
  33. """
  34. Returns the components of this vector in the form of a
  35. Python dictionary mapping BaseVector instances to the
  36. corresponding measure numbers.
  37. Examples
  38. ========
  39. >>> from sympy.vector import CoordSys3D
  40. >>> C = CoordSys3D('C')
  41. >>> v = 3*C.i + 4*C.j + 5*C.k
  42. >>> v.components
  43. {C.i: 3, C.j: 4, C.k: 5}
  44. """
  45. # The '_components' attribute is defined according to the
  46. # subclass of Vector the instance belongs to.
  47. return self._components
  48. def magnitude(self):
  49. """
  50. Returns the magnitude of this vector.
  51. """
  52. return sqrt(self & self)
  53. def normalize(self):
  54. """
  55. Returns the normalized version of this vector.
  56. """
  57. return self / self.magnitude()
  58. def dot(self, other):
  59. """
  60. Returns the dot product of this Vector, either with another
  61. Vector, or a Dyadic, or a Del operator.
  62. If 'other' is a Vector, returns the dot product scalar (SymPy
  63. expression).
  64. If 'other' is a Dyadic, the dot product is returned as a Vector.
  65. If 'other' is an instance of Del, returns the directional
  66. derivative operator as a Python function. If this function is
  67. applied to a scalar expression, it returns the directional
  68. derivative of the scalar field wrt this Vector.
  69. Parameters
  70. ==========
  71. other: Vector/Dyadic/Del
  72. The Vector or Dyadic we are dotting with, or a Del operator .
  73. Examples
  74. ========
  75. >>> from sympy.vector import CoordSys3D, Del
  76. >>> C = CoordSys3D('C')
  77. >>> delop = Del()
  78. >>> C.i.dot(C.j)
  79. 0
  80. >>> C.i & C.i
  81. 1
  82. >>> v = 3*C.i + 4*C.j + 5*C.k
  83. >>> v.dot(C.k)
  84. 5
  85. >>> (C.i & delop)(C.x*C.y*C.z)
  86. C.y*C.z
  87. >>> d = C.i.outer(C.i)
  88. >>> C.i.dot(d)
  89. C.i
  90. """
  91. # Check special cases
  92. if isinstance(other, Dyadic):
  93. if isinstance(self, VectorZero):
  94. return Vector.zero
  95. outvec = Vector.zero
  96. for k, v in other.components.items():
  97. vect_dot = k.args[0].dot(self)
  98. outvec += vect_dot * v * k.args[1]
  99. return outvec
  100. from sympy.vector.deloperator import Del
  101. if not isinstance(other, (Del, Vector)):
  102. raise TypeError(str(other) + " is not a vector, dyadic or " +
  103. "del operator")
  104. # Check if the other is a del operator
  105. if isinstance(other, Del):
  106. def directional_derivative(field):
  107. from sympy.vector.functions import directional_derivative
  108. return directional_derivative(field, self)
  109. return directional_derivative
  110. return dot(self, other)
  111. def __and__(self, other):
  112. return self.dot(other)
  113. __and__.__doc__ = dot.__doc__
  114. def cross(self, other):
  115. """
  116. Returns the cross product of this Vector with another Vector or
  117. Dyadic instance.
  118. The cross product is a Vector, if 'other' is a Vector. If 'other'
  119. is a Dyadic, this returns a Dyadic instance.
  120. Parameters
  121. ==========
  122. other: Vector/Dyadic
  123. The Vector or Dyadic we are crossing with.
  124. Examples
  125. ========
  126. >>> from sympy.vector import CoordSys3D
  127. >>> C = CoordSys3D('C')
  128. >>> C.i.cross(C.j)
  129. C.k
  130. >>> C.i ^ C.i
  131. 0
  132. >>> v = 3*C.i + 4*C.j + 5*C.k
  133. >>> v ^ C.i
  134. 5*C.j + (-4)*C.k
  135. >>> d = C.i.outer(C.i)
  136. >>> C.j.cross(d)
  137. (-1)*(C.k|C.i)
  138. """
  139. # Check special cases
  140. if isinstance(other, Dyadic):
  141. if isinstance(self, VectorZero):
  142. return Dyadic.zero
  143. outdyad = Dyadic.zero
  144. for k, v in other.components.items():
  145. cross_product = self.cross(k.args[0])
  146. outer = cross_product.outer(k.args[1])
  147. outdyad += v * outer
  148. return outdyad
  149. return cross(self, other)
  150. def __xor__(self, other):
  151. return self.cross(other)
  152. __xor__.__doc__ = cross.__doc__
  153. def outer(self, other):
  154. """
  155. Returns the outer product of this vector with another, in the
  156. form of a Dyadic instance.
  157. Parameters
  158. ==========
  159. other : Vector
  160. The Vector with respect to which the outer product is to
  161. be computed.
  162. Examples
  163. ========
  164. >>> from sympy.vector import CoordSys3D
  165. >>> N = CoordSys3D('N')
  166. >>> N.i.outer(N.j)
  167. (N.i|N.j)
  168. """
  169. # Handle the special cases
  170. if not isinstance(other, Vector):
  171. raise TypeError("Invalid operand for outer product")
  172. elif (isinstance(self, VectorZero) or
  173. isinstance(other, VectorZero)):
  174. return Dyadic.zero
  175. # Iterate over components of both the vectors to generate
  176. # the required Dyadic instance
  177. args = [(v1 * v2) * BaseDyadic(k1, k2) for (k1, v1), (k2, v2)
  178. in product(self.components.items(), other.components.items())]
  179. return DyadicAdd(*args)
  180. def projection(self, other, scalar=False):
  181. """
  182. Returns the vector or scalar projection of the 'other' on 'self'.
  183. Examples
  184. ========
  185. >>> from sympy.vector.coordsysrect import CoordSys3D
  186. >>> C = CoordSys3D('C')
  187. >>> i, j, k = C.base_vectors()
  188. >>> v1 = i + j + k
  189. >>> v2 = 3*i + 4*j
  190. >>> v1.projection(v2)
  191. 7/3*C.i + 7/3*C.j + 7/3*C.k
  192. >>> v1.projection(v2, scalar=True)
  193. 7/3
  194. """
  195. if self.equals(Vector.zero):
  196. return S.Zero if scalar else Vector.zero
  197. if scalar:
  198. return self.dot(other) / self.dot(self)
  199. else:
  200. return self.dot(other) / self.dot(self) * self
  201. @property
  202. def _projections(self):
  203. """
  204. Returns the components of this vector but the output includes
  205. also zero values components.
  206. Examples
  207. ========
  208. >>> from sympy.vector import CoordSys3D, Vector
  209. >>> C = CoordSys3D('C')
  210. >>> v1 = 3*C.i + 4*C.j + 5*C.k
  211. >>> v1._projections
  212. (3, 4, 5)
  213. >>> v2 = C.x*C.y*C.z*C.i
  214. >>> v2._projections
  215. (C.x*C.y*C.z, 0, 0)
  216. >>> v3 = Vector.zero
  217. >>> v3._projections
  218. (0, 0, 0)
  219. """
  220. from sympy.vector.operators import _get_coord_systems
  221. if isinstance(self, VectorZero):
  222. return (S.Zero, S.Zero, S.Zero)
  223. base_vec = next(iter(_get_coord_systems(self))).base_vectors()
  224. return tuple([self.dot(i) for i in base_vec])
  225. def __or__(self, other):
  226. return self.outer(other)
  227. __or__.__doc__ = outer.__doc__
  228. def to_matrix(self, system):
  229. """
  230. Returns the matrix form of this vector with respect to the
  231. specified coordinate system.
  232. Parameters
  233. ==========
  234. system : CoordSys3D
  235. The system wrt which the matrix form is to be computed
  236. Examples
  237. ========
  238. >>> from sympy.vector import CoordSys3D
  239. >>> C = CoordSys3D('C')
  240. >>> from sympy.abc import a, b, c
  241. >>> v = a*C.i + b*C.j + c*C.k
  242. >>> v.to_matrix(C)
  243. Matrix([
  244. [a],
  245. [b],
  246. [c]])
  247. """
  248. return Matrix([self.dot(unit_vec) for unit_vec in
  249. system.base_vectors()])
  250. def separate(self):
  251. """
  252. The constituents of this vector in different coordinate systems,
  253. as per its definition.
  254. Returns a dict mapping each CoordSys3D to the corresponding
  255. constituent Vector.
  256. Examples
  257. ========
  258. >>> from sympy.vector import CoordSys3D
  259. >>> R1 = CoordSys3D('R1')
  260. >>> R2 = CoordSys3D('R2')
  261. >>> v = R1.i + R2.i
  262. >>> v.separate() == {R1: R1.i, R2: R2.i}
  263. True
  264. """
  265. parts = {}
  266. for vect, measure in self.components.items():
  267. parts[vect.system] = (parts.get(vect.system, Vector.zero) +
  268. vect * measure)
  269. return parts
  270. def _div_helper(one, other):
  271. """ Helper for division involving vectors. """
  272. if isinstance(one, Vector) and isinstance(other, Vector):
  273. raise TypeError("Cannot divide two vectors")
  274. elif isinstance(one, Vector):
  275. if other == S.Zero:
  276. raise ValueError("Cannot divide a vector by zero")
  277. return VectorMul(one, Pow(other, S.NegativeOne))
  278. else:
  279. raise TypeError("Invalid division involving a vector")
  280. class BaseVector(Vector, AtomicExpr):
  281. """
  282. Class to denote a base vector.
  283. """
  284. def __new__(cls, index, system, pretty_str=None, latex_str=None):
  285. if pretty_str is None:
  286. pretty_str = "x{}".format(index)
  287. if latex_str is None:
  288. latex_str = "x_{}".format(index)
  289. pretty_str = str(pretty_str)
  290. latex_str = str(latex_str)
  291. # Verify arguments
  292. if index not in range(0, 3):
  293. raise ValueError("index must be 0, 1 or 2")
  294. if not isinstance(system, CoordSys3D):
  295. raise TypeError("system should be a CoordSys3D")
  296. name = system._vector_names[index]
  297. # Initialize an object
  298. obj = super().__new__(cls, S(index), system)
  299. # Assign important attributes
  300. obj._base_instance = obj
  301. obj._components = {obj: S.One}
  302. obj._measure_number = S.One
  303. obj._name = system._name + '.' + name
  304. obj._pretty_form = '' + pretty_str
  305. obj._latex_form = latex_str
  306. obj._system = system
  307. # The _id is used for printing purposes
  308. obj._id = (index, system)
  309. assumptions = {'commutative': True}
  310. obj._assumptions = StdFactKB(assumptions)
  311. # This attr is used for re-expression to one of the systems
  312. # involved in the definition of the Vector. Applies to
  313. # VectorMul and VectorAdd too.
  314. obj._sys = system
  315. return obj
  316. @property
  317. def system(self):
  318. return self._system
  319. def _sympystr(self, printer):
  320. return self._name
  321. def _sympyrepr(self, printer):
  322. index, system = self._id
  323. return printer._print(system) + '.' + system._vector_names[index]
  324. @property
  325. def free_symbols(self):
  326. return {self}
  327. class VectorAdd(BasisDependentAdd, Vector):
  328. """
  329. Class to denote sum of Vector instances.
  330. """
  331. def __new__(cls, *args, **options):
  332. obj = BasisDependentAdd.__new__(cls, *args, **options)
  333. return obj
  334. def _sympystr(self, printer):
  335. ret_str = ''
  336. items = list(self.separate().items())
  337. items.sort(key=lambda x: x[0].__str__())
  338. for system, vect in items:
  339. base_vects = system.base_vectors()
  340. for x in base_vects:
  341. if x in vect.components:
  342. temp_vect = self.components[x] * x
  343. ret_str += printer._print(temp_vect) + " + "
  344. return ret_str[:-3]
  345. class VectorMul(BasisDependentMul, Vector):
  346. """
  347. Class to denote products of scalars and BaseVectors.
  348. """
  349. def __new__(cls, *args, **options):
  350. obj = BasisDependentMul.__new__(cls, *args, **options)
  351. return obj
  352. @property
  353. def base_vector(self):
  354. """ The BaseVector involved in the product. """
  355. return self._base_instance
  356. @property
  357. def measure_number(self):
  358. """ The scalar expression involved in the definition of
  359. this VectorMul.
  360. """
  361. return self._measure_number
  362. class VectorZero(BasisDependentZero, Vector):
  363. """
  364. Class to denote a zero vector
  365. """
  366. _op_priority = 12.1
  367. _pretty_form = '0'
  368. _latex_form = r'\mathbf{\hat{0}}'
  369. def __new__(cls):
  370. obj = BasisDependentZero.__new__(cls)
  371. return obj
  372. class Cross(Vector):
  373. """
  374. Represents unevaluated Cross product.
  375. Examples
  376. ========
  377. >>> from sympy.vector import CoordSys3D, Cross
  378. >>> R = CoordSys3D('R')
  379. >>> v1 = R.i + R.j + R.k
  380. >>> v2 = R.x * R.i + R.y * R.j + R.z * R.k
  381. >>> Cross(v1, v2)
  382. Cross(R.i + R.j + R.k, R.x*R.i + R.y*R.j + R.z*R.k)
  383. >>> Cross(v1, v2).doit()
  384. (-R.y + R.z)*R.i + (R.x - R.z)*R.j + (-R.x + R.y)*R.k
  385. """
  386. def __new__(cls, expr1, expr2):
  387. expr1 = sympify(expr1)
  388. expr2 = sympify(expr2)
  389. if default_sort_key(expr1) > default_sort_key(expr2):
  390. return -Cross(expr2, expr1)
  391. obj = Expr.__new__(cls, expr1, expr2)
  392. obj._expr1 = expr1
  393. obj._expr2 = expr2
  394. return obj
  395. def doit(self, **hints):
  396. return cross(self._expr1, self._expr2)
  397. class Dot(Expr):
  398. """
  399. Represents unevaluated Dot product.
  400. Examples
  401. ========
  402. >>> from sympy.vector import CoordSys3D, Dot
  403. >>> from sympy import symbols
  404. >>> R = CoordSys3D('R')
  405. >>> a, b, c = symbols('a b c')
  406. >>> v1 = R.i + R.j + R.k
  407. >>> v2 = a * R.i + b * R.j + c * R.k
  408. >>> Dot(v1, v2)
  409. Dot(R.i + R.j + R.k, a*R.i + b*R.j + c*R.k)
  410. >>> Dot(v1, v2).doit()
  411. a + b + c
  412. """
  413. def __new__(cls, expr1, expr2):
  414. expr1 = sympify(expr1)
  415. expr2 = sympify(expr2)
  416. expr1, expr2 = sorted([expr1, expr2], key=default_sort_key)
  417. obj = Expr.__new__(cls, expr1, expr2)
  418. obj._expr1 = expr1
  419. obj._expr2 = expr2
  420. return obj
  421. def doit(self, **hints):
  422. return dot(self._expr1, self._expr2)
  423. def cross(vect1, vect2):
  424. """
  425. Returns cross product of two vectors.
  426. Examples
  427. ========
  428. >>> from sympy.vector import CoordSys3D
  429. >>> from sympy.vector.vector import cross
  430. >>> R = CoordSys3D('R')
  431. >>> v1 = R.i + R.j + R.k
  432. >>> v2 = R.x * R.i + R.y * R.j + R.z * R.k
  433. >>> cross(v1, v2)
  434. (-R.y + R.z)*R.i + (R.x - R.z)*R.j + (-R.x + R.y)*R.k
  435. """
  436. if isinstance(vect1, Add):
  437. return VectorAdd.fromiter(cross(i, vect2) for i in vect1.args)
  438. if isinstance(vect2, Add):
  439. return VectorAdd.fromiter(cross(vect1, i) for i in vect2.args)
  440. if isinstance(vect1, BaseVector) and isinstance(vect2, BaseVector):
  441. if vect1._sys == vect2._sys:
  442. n1 = vect1.args[0]
  443. n2 = vect2.args[0]
  444. if n1 == n2:
  445. return Vector.zero
  446. n3 = ({0,1,2}.difference({n1, n2})).pop()
  447. sign = 1 if ((n1 + 1) % 3 == n2) else -1
  448. return sign*vect1._sys.base_vectors()[n3]
  449. from .functions import express
  450. try:
  451. v = express(vect1, vect2._sys)
  452. except ValueError:
  453. return Cross(vect1, vect2)
  454. else:
  455. return cross(v, vect2)
  456. if isinstance(vect1, VectorZero) or isinstance(vect2, VectorZero):
  457. return Vector.zero
  458. if isinstance(vect1, VectorMul):
  459. v1, m1 = next(iter(vect1.components.items()))
  460. return m1*cross(v1, vect2)
  461. if isinstance(vect2, VectorMul):
  462. v2, m2 = next(iter(vect2.components.items()))
  463. return m2*cross(vect1, v2)
  464. return Cross(vect1, vect2)
  465. def dot(vect1, vect2):
  466. """
  467. Returns dot product of two vectors.
  468. Examples
  469. ========
  470. >>> from sympy.vector import CoordSys3D
  471. >>> from sympy.vector.vector import dot
  472. >>> R = CoordSys3D('R')
  473. >>> v1 = R.i + R.j + R.k
  474. >>> v2 = R.x * R.i + R.y * R.j + R.z * R.k
  475. >>> dot(v1, v2)
  476. R.x + R.y + R.z
  477. """
  478. if isinstance(vect1, Add):
  479. return Add.fromiter(dot(i, vect2) for i in vect1.args)
  480. if isinstance(vect2, Add):
  481. return Add.fromiter(dot(vect1, i) for i in vect2.args)
  482. if isinstance(vect1, BaseVector) and isinstance(vect2, BaseVector):
  483. if vect1._sys == vect2._sys:
  484. return S.One if vect1 == vect2 else S.Zero
  485. from .functions import express
  486. try:
  487. v = express(vect2, vect1._sys)
  488. except ValueError:
  489. return Dot(vect1, vect2)
  490. else:
  491. return dot(vect1, v)
  492. if isinstance(vect1, VectorZero) or isinstance(vect2, VectorZero):
  493. return S.Zero
  494. if isinstance(vect1, VectorMul):
  495. v1, m1 = next(iter(vect1.components.items()))
  496. return m1*dot(v1, vect2)
  497. if isinstance(vect2, VectorMul):
  498. v2, m2 = next(iter(vect2.components.items()))
  499. return m2*dot(vect1, v2)
  500. return Dot(vect1, vect2)
  501. Vector._expr_type = Vector
  502. Vector._mul_func = VectorMul
  503. Vector._add_func = VectorAdd
  504. Vector._zero_func = VectorZero
  505. Vector._base_func = BaseVector
  506. Vector.zero = VectorZero()