dyadic.py 18 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545
  1. from sympy import sympify, Add, ImmutableMatrix as Matrix
  2. from sympy.core.evalf import EvalfMixin
  3. from sympy.printing.defaults import Printable
  4. from mpmath.libmp.libmpf import prec_to_dps
  5. __all__ = ['Dyadic']
  6. class Dyadic(Printable, EvalfMixin):
  7. """A Dyadic object.
  8. See:
  9. https://en.wikipedia.org/wiki/Dyadic_tensor
  10. Kane, T., Levinson, D. Dynamics Theory and Applications. 1985 McGraw-Hill
  11. A more powerful way to represent a rigid body's inertia. While it is more
  12. complex, by choosing Dyadic components to be in body fixed basis vectors,
  13. the resulting matrix is equivalent to the inertia tensor.
  14. """
  15. is_number = False
  16. def __init__(self, inlist):
  17. """
  18. Just like Vector's init, you should not call this unless creating a
  19. zero dyadic.
  20. zd = Dyadic(0)
  21. Stores a Dyadic as a list of lists; the inner list has the measure
  22. number and the two unit vectors; the outerlist holds each unique
  23. unit vector pair.
  24. """
  25. self.args = []
  26. if inlist == 0:
  27. inlist = []
  28. while len(inlist) != 0:
  29. added = 0
  30. for i, v in enumerate(self.args):
  31. if ((str(inlist[0][1]) == str(self.args[i][1])) and
  32. (str(inlist[0][2]) == str(self.args[i][2]))):
  33. self.args[i] = (self.args[i][0] + inlist[0][0],
  34. inlist[0][1], inlist[0][2])
  35. inlist.remove(inlist[0])
  36. added = 1
  37. break
  38. if added != 1:
  39. self.args.append(inlist[0])
  40. inlist.remove(inlist[0])
  41. i = 0
  42. # This code is to remove empty parts from the list
  43. while i < len(self.args):
  44. if ((self.args[i][0] == 0) | (self.args[i][1] == 0) |
  45. (self.args[i][2] == 0)):
  46. self.args.remove(self.args[i])
  47. i -= 1
  48. i += 1
  49. @property
  50. def func(self):
  51. """Returns the class Dyadic. """
  52. return Dyadic
  53. def __add__(self, other):
  54. """The add operator for Dyadic. """
  55. other = _check_dyadic(other)
  56. return Dyadic(self.args + other.args)
  57. __radd__ = __add__
  58. def __mul__(self, other):
  59. """Multiplies the Dyadic by a sympifyable expression.
  60. Parameters
  61. ==========
  62. other : Sympafiable
  63. The scalar to multiply this Dyadic with
  64. Examples
  65. ========
  66. >>> from sympy.physics.vector import ReferenceFrame, outer
  67. >>> N = ReferenceFrame('N')
  68. >>> d = outer(N.x, N.x)
  69. >>> 5 * d
  70. 5*(N.x|N.x)
  71. """
  72. newlist = list(self.args)
  73. other = sympify(other)
  74. for i in range(len(newlist)):
  75. newlist[i] = (other * newlist[i][0], newlist[i][1],
  76. newlist[i][2])
  77. return Dyadic(newlist)
  78. __rmul__ = __mul__
  79. def dot(self, other):
  80. """The inner product operator for a Dyadic and a Dyadic or Vector.
  81. Parameters
  82. ==========
  83. other : Dyadic or Vector
  84. The other Dyadic or Vector to take the inner product with
  85. Examples
  86. ========
  87. >>> from sympy.physics.vector import ReferenceFrame, outer
  88. >>> N = ReferenceFrame('N')
  89. >>> D1 = outer(N.x, N.y)
  90. >>> D2 = outer(N.y, N.y)
  91. >>> D1.dot(D2)
  92. (N.x|N.y)
  93. >>> D1.dot(N.y)
  94. N.x
  95. """
  96. from sympy.physics.vector.vector import Vector, _check_vector
  97. if isinstance(other, Dyadic):
  98. other = _check_dyadic(other)
  99. ol = Dyadic(0)
  100. for v in self.args:
  101. for v2 in other.args:
  102. ol += v[0] * v2[0] * (v[2].dot(v2[1])) * (v[1].outer(v2[2]))
  103. else:
  104. other = _check_vector(other)
  105. ol = Vector(0)
  106. for v in self.args:
  107. ol += v[0] * v[1] * (v[2].dot(other))
  108. return ol
  109. # NOTE : supports non-advertised Dyadic & Dyadic, Dyadic & Vector notation
  110. __and__ = dot
  111. def __truediv__(self, other):
  112. """Divides the Dyadic by a sympifyable expression. """
  113. return self.__mul__(1 / other)
  114. def __eq__(self, other):
  115. """Tests for equality.
  116. Is currently weak; needs stronger comparison testing
  117. """
  118. if other == 0:
  119. other = Dyadic(0)
  120. other = _check_dyadic(other)
  121. if (self.args == []) and (other.args == []):
  122. return True
  123. elif (self.args == []) or (other.args == []):
  124. return False
  125. return set(self.args) == set(other.args)
  126. def __ne__(self, other):
  127. return not self == other
  128. def __neg__(self):
  129. return self * -1
  130. def _latex(self, printer):
  131. ar = self.args # just to shorten things
  132. if len(ar) == 0:
  133. return str(0)
  134. ol = [] # output list, to be concatenated to a string
  135. for v in ar:
  136. # if the coef of the dyadic is 1, we skip the 1
  137. if v[0] == 1:
  138. ol.append(' + ' + printer._print(v[1]) + r"\otimes " +
  139. printer._print(v[2]))
  140. # if the coef of the dyadic is -1, we skip the 1
  141. elif v[0] == -1:
  142. ol.append(' - ' +
  143. printer._print(v[1]) +
  144. r"\otimes " +
  145. printer._print(v[2]))
  146. # If the coefficient of the dyadic is not 1 or -1,
  147. # we might wrap it in parentheses, for readability.
  148. elif v[0] != 0:
  149. arg_str = printer._print(v[0])
  150. if isinstance(v[0], Add):
  151. arg_str = '(%s)' % arg_str
  152. if arg_str.startswith('-'):
  153. arg_str = arg_str[1:]
  154. str_start = ' - '
  155. else:
  156. str_start = ' + '
  157. ol.append(str_start + arg_str + printer._print(v[1]) +
  158. r"\otimes " + printer._print(v[2]))
  159. outstr = ''.join(ol)
  160. if outstr.startswith(' + '):
  161. outstr = outstr[3:]
  162. elif outstr.startswith(' '):
  163. outstr = outstr[1:]
  164. return outstr
  165. def _pretty(self, printer):
  166. e = self
  167. class Fake:
  168. baseline = 0
  169. def render(self, *args, **kwargs):
  170. ar = e.args # just to shorten things
  171. mpp = printer
  172. if len(ar) == 0:
  173. return str(0)
  174. bar = "\N{CIRCLED TIMES}" if printer._use_unicode else "|"
  175. ol = [] # output list, to be concatenated to a string
  176. for v in ar:
  177. # if the coef of the dyadic is 1, we skip the 1
  178. if v[0] == 1:
  179. ol.extend([" + ",
  180. mpp.doprint(v[1]),
  181. bar,
  182. mpp.doprint(v[2])])
  183. # if the coef of the dyadic is -1, we skip the 1
  184. elif v[0] == -1:
  185. ol.extend([" - ",
  186. mpp.doprint(v[1]),
  187. bar,
  188. mpp.doprint(v[2])])
  189. # If the coefficient of the dyadic is not 1 or -1,
  190. # we might wrap it in parentheses, for readability.
  191. elif v[0] != 0:
  192. if isinstance(v[0], Add):
  193. arg_str = mpp._print(
  194. v[0]).parens()[0]
  195. else:
  196. arg_str = mpp.doprint(v[0])
  197. if arg_str.startswith("-"):
  198. arg_str = arg_str[1:]
  199. str_start = " - "
  200. else:
  201. str_start = " + "
  202. ol.extend([str_start, arg_str, " ",
  203. mpp.doprint(v[1]),
  204. bar,
  205. mpp.doprint(v[2])])
  206. outstr = "".join(ol)
  207. if outstr.startswith(" + "):
  208. outstr = outstr[3:]
  209. elif outstr.startswith(" "):
  210. outstr = outstr[1:]
  211. return outstr
  212. return Fake()
  213. def __rsub__(self, other):
  214. return (-1 * self) + other
  215. def _sympystr(self, printer):
  216. """Printing method. """
  217. ar = self.args # just to shorten things
  218. if len(ar) == 0:
  219. return printer._print(0)
  220. ol = [] # output list, to be concatenated to a string
  221. for v in ar:
  222. # if the coef of the dyadic is 1, we skip the 1
  223. if v[0] == 1:
  224. ol.append(' + (' + printer._print(v[1]) + '|' +
  225. printer._print(v[2]) + ')')
  226. # if the coef of the dyadic is -1, we skip the 1
  227. elif v[0] == -1:
  228. ol.append(' - (' + printer._print(v[1]) + '|' +
  229. printer._print(v[2]) + ')')
  230. # If the coefficient of the dyadic is not 1 or -1,
  231. # we might wrap it in parentheses, for readability.
  232. elif v[0] != 0:
  233. arg_str = printer._print(v[0])
  234. if isinstance(v[0], Add):
  235. arg_str = "(%s)" % arg_str
  236. if arg_str[0] == '-':
  237. arg_str = arg_str[1:]
  238. str_start = ' - '
  239. else:
  240. str_start = ' + '
  241. ol.append(str_start + arg_str + '*(' +
  242. printer._print(v[1]) +
  243. '|' + printer._print(v[2]) + ')')
  244. outstr = ''.join(ol)
  245. if outstr.startswith(' + '):
  246. outstr = outstr[3:]
  247. elif outstr.startswith(' '):
  248. outstr = outstr[1:]
  249. return outstr
  250. def __sub__(self, other):
  251. """The subtraction operator. """
  252. return self.__add__(other * -1)
  253. def cross(self, other):
  254. """Returns the dyadic resulting from the dyadic vector cross product:
  255. Dyadic x Vector.
  256. Parameters
  257. ==========
  258. other : Vector
  259. Vector to cross with.
  260. Examples
  261. ========
  262. >>> from sympy.physics.vector import ReferenceFrame, outer, cross
  263. >>> N = ReferenceFrame('N')
  264. >>> d = outer(N.x, N.x)
  265. >>> cross(d, N.y)
  266. (N.x|N.z)
  267. """
  268. from sympy.physics.vector.vector import _check_vector
  269. other = _check_vector(other)
  270. ol = Dyadic(0)
  271. for v in self.args:
  272. ol += v[0] * (v[1].outer((v[2].cross(other))))
  273. return ol
  274. # NOTE : supports non-advertised Dyadic ^ Vector notation
  275. __xor__ = cross
  276. def express(self, frame1, frame2=None):
  277. """Expresses this Dyadic in alternate frame(s)
  278. The first frame is the list side expression, the second frame is the
  279. right side; if Dyadic is in form A.x|B.y, you can express it in two
  280. different frames. If no second frame is given, the Dyadic is
  281. expressed in only one frame.
  282. Calls the global express function
  283. Parameters
  284. ==========
  285. frame1 : ReferenceFrame
  286. The frame to express the left side of the Dyadic in
  287. frame2 : ReferenceFrame
  288. If provided, the frame to express the right side of the Dyadic in
  289. Examples
  290. ========
  291. >>> from sympy.physics.vector import ReferenceFrame, outer, dynamicsymbols
  292. >>> from sympy.physics.vector import init_vprinting
  293. >>> init_vprinting(pretty_print=False)
  294. >>> N = ReferenceFrame('N')
  295. >>> q = dynamicsymbols('q')
  296. >>> B = N.orientnew('B', 'Axis', [q, N.z])
  297. >>> d = outer(N.x, N.x)
  298. >>> d.express(B, N)
  299. cos(q)*(B.x|N.x) - sin(q)*(B.y|N.x)
  300. """
  301. from sympy.physics.vector.functions import express
  302. return express(self, frame1, frame2)
  303. def to_matrix(self, reference_frame, second_reference_frame=None):
  304. """Returns the matrix form of the dyadic with respect to one or two
  305. reference frames.
  306. Parameters
  307. ----------
  308. reference_frame : ReferenceFrame
  309. The reference frame that the rows and columns of the matrix
  310. correspond to. If a second reference frame is provided, this
  311. only corresponds to the rows of the matrix.
  312. second_reference_frame : ReferenceFrame, optional, default=None
  313. The reference frame that the columns of the matrix correspond
  314. to.
  315. Returns
  316. -------
  317. matrix : ImmutableMatrix, shape(3,3)
  318. The matrix that gives the 2D tensor form.
  319. Examples
  320. ========
  321. >>> from sympy import symbols, trigsimp
  322. >>> from sympy.physics.vector import ReferenceFrame
  323. >>> from sympy.physics.mechanics import inertia
  324. >>> Ixx, Iyy, Izz, Ixy, Iyz, Ixz = symbols('Ixx, Iyy, Izz, Ixy, Iyz, Ixz')
  325. >>> N = ReferenceFrame('N')
  326. >>> inertia_dyadic = inertia(N, Ixx, Iyy, Izz, Ixy, Iyz, Ixz)
  327. >>> inertia_dyadic.to_matrix(N)
  328. Matrix([
  329. [Ixx, Ixy, Ixz],
  330. [Ixy, Iyy, Iyz],
  331. [Ixz, Iyz, Izz]])
  332. >>> beta = symbols('beta')
  333. >>> A = N.orientnew('A', 'Axis', (beta, N.x))
  334. >>> trigsimp(inertia_dyadic.to_matrix(A))
  335. Matrix([
  336. [ Ixx, Ixy*cos(beta) + Ixz*sin(beta), -Ixy*sin(beta) + Ixz*cos(beta)],
  337. [ Ixy*cos(beta) + Ixz*sin(beta), Iyy*cos(2*beta)/2 + Iyy/2 + Iyz*sin(2*beta) - Izz*cos(2*beta)/2 + Izz/2, -Iyy*sin(2*beta)/2 + Iyz*cos(2*beta) + Izz*sin(2*beta)/2],
  338. [-Ixy*sin(beta) + Ixz*cos(beta), -Iyy*sin(2*beta)/2 + Iyz*cos(2*beta) + Izz*sin(2*beta)/2, -Iyy*cos(2*beta)/2 + Iyy/2 - Iyz*sin(2*beta) + Izz*cos(2*beta)/2 + Izz/2]])
  339. """
  340. if second_reference_frame is None:
  341. second_reference_frame = reference_frame
  342. return Matrix([i.dot(self).dot(j) for i in reference_frame for j in
  343. second_reference_frame]).reshape(3, 3)
  344. def doit(self, **hints):
  345. """Calls .doit() on each term in the Dyadic"""
  346. return sum([Dyadic([(v[0].doit(**hints), v[1], v[2])])
  347. for v in self.args], Dyadic(0))
  348. def dt(self, frame):
  349. """Take the time derivative of this Dyadic in a frame.
  350. This function calls the global time_derivative method
  351. Parameters
  352. ==========
  353. frame : ReferenceFrame
  354. The frame to take the time derivative in
  355. Examples
  356. ========
  357. >>> from sympy.physics.vector import ReferenceFrame, outer, dynamicsymbols
  358. >>> from sympy.physics.vector import init_vprinting
  359. >>> init_vprinting(pretty_print=False)
  360. >>> N = ReferenceFrame('N')
  361. >>> q = dynamicsymbols('q')
  362. >>> B = N.orientnew('B', 'Axis', [q, N.z])
  363. >>> d = outer(N.x, N.x)
  364. >>> d.dt(B)
  365. - q'*(N.y|N.x) - q'*(N.x|N.y)
  366. """
  367. from sympy.physics.vector.functions import time_derivative
  368. return time_derivative(self, frame)
  369. def simplify(self):
  370. """Returns a simplified Dyadic."""
  371. out = Dyadic(0)
  372. for v in self.args:
  373. out += Dyadic([(v[0].simplify(), v[1], v[2])])
  374. return out
  375. def subs(self, *args, **kwargs):
  376. """Substitution on the Dyadic.
  377. Examples
  378. ========
  379. >>> from sympy.physics.vector import ReferenceFrame
  380. >>> from sympy import Symbol
  381. >>> N = ReferenceFrame('N')
  382. >>> s = Symbol('s')
  383. >>> a = s*(N.x|N.x)
  384. >>> a.subs({s: 2})
  385. 2*(N.x|N.x)
  386. """
  387. return sum([Dyadic([(v[0].subs(*args, **kwargs), v[1], v[2])])
  388. for v in self.args], Dyadic(0))
  389. def applyfunc(self, f):
  390. """Apply a function to each component of a Dyadic."""
  391. if not callable(f):
  392. raise TypeError("`f` must be callable.")
  393. out = Dyadic(0)
  394. for a, b, c in self.args:
  395. out += f(a) * (b.outer(c))
  396. return out
  397. def _eval_evalf(self, prec):
  398. if not self.args:
  399. return self
  400. new_args = []
  401. dps = prec_to_dps(prec)
  402. for inlist in self.args:
  403. new_inlist = list(inlist)
  404. new_inlist[0] = inlist[0].evalf(n=dps)
  405. new_args.append(tuple(new_inlist))
  406. return Dyadic(new_args)
  407. def xreplace(self, rule):
  408. """
  409. Replace occurrences of objects within the measure numbers of the
  410. Dyadic.
  411. Parameters
  412. ==========
  413. rule : dict-like
  414. Expresses a replacement rule.
  415. Returns
  416. =======
  417. Dyadic
  418. Result of the replacement.
  419. Examples
  420. ========
  421. >>> from sympy import symbols, pi
  422. >>> from sympy.physics.vector import ReferenceFrame, outer
  423. >>> N = ReferenceFrame('N')
  424. >>> D = outer(N.x, N.x)
  425. >>> x, y, z = symbols('x y z')
  426. >>> ((1 + x*y) * D).xreplace({x: pi})
  427. (pi*y + 1)*(N.x|N.x)
  428. >>> ((1 + x*y) * D).xreplace({x: pi, y: 2})
  429. (1 + 2*pi)*(N.x|N.x)
  430. Replacements occur only if an entire node in the expression tree is
  431. matched:
  432. >>> ((x*y + z) * D).xreplace({x*y: pi})
  433. (z + pi)*(N.x|N.x)
  434. >>> ((x*y*z) * D).xreplace({x*y: pi})
  435. x*y*z*(N.x|N.x)
  436. """
  437. new_args = []
  438. for inlist in self.args:
  439. new_inlist = list(inlist)
  440. new_inlist[0] = new_inlist[0].xreplace(rule)
  441. new_args.append(tuple(new_inlist))
  442. return Dyadic(new_args)
  443. def _check_dyadic(other):
  444. if not isinstance(other, Dyadic):
  445. raise TypeError('A Dyadic must be supplied')
  446. return other