vector.py 20 KB

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