functions.py 25 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735
  1. from sympy.utilities import dict_merge
  2. from sympy.utilities.iterables import iterable
  3. from sympy.physics.vector import (Dyadic, Vector, ReferenceFrame,
  4. Point, dynamicsymbols)
  5. from sympy.physics.vector.printing import (vprint, vsprint, vpprint, vlatex,
  6. init_vprinting)
  7. from sympy.physics.mechanics.particle import Particle
  8. from sympy.physics.mechanics.rigidbody import RigidBody
  9. from sympy.simplify.simplify import simplify
  10. from sympy import Matrix, Mul, Derivative, sin, cos, tan, S
  11. from sympy.core.function import AppliedUndef
  12. from sympy.physics.mechanics.inertia import (inertia as _inertia,
  13. inertia_of_point_mass as _inertia_of_point_mass)
  14. from sympy.utilities.exceptions import sympy_deprecation_warning
  15. __all__ = ['linear_momentum',
  16. 'angular_momentum',
  17. 'kinetic_energy',
  18. 'potential_energy',
  19. 'Lagrangian',
  20. 'mechanics_printing',
  21. 'mprint',
  22. 'msprint',
  23. 'mpprint',
  24. 'mlatex',
  25. 'msubs',
  26. 'find_dynamicsymbols']
  27. # These are functions that we've moved and renamed during extracting the
  28. # basic vector calculus code from the mechanics packages.
  29. mprint = vprint
  30. msprint = vsprint
  31. mpprint = vpprint
  32. mlatex = vlatex
  33. def mechanics_printing(**kwargs):
  34. """
  35. Initializes time derivative printing for all SymPy objects in
  36. mechanics module.
  37. """
  38. init_vprinting(**kwargs)
  39. mechanics_printing.__doc__ = init_vprinting.__doc__
  40. def inertia(frame, ixx, iyy, izz, ixy=0, iyz=0, izx=0):
  41. sympy_deprecation_warning(
  42. """
  43. The inertia function has been moved.
  44. Import it from "sympy.physics.mechanics".
  45. """,
  46. deprecated_since_version="1.13",
  47. active_deprecations_target="moved-mechanics-functions"
  48. )
  49. return _inertia(frame, ixx, iyy, izz, ixy, iyz, izx)
  50. def inertia_of_point_mass(mass, pos_vec, frame):
  51. sympy_deprecation_warning(
  52. """
  53. The inertia_of_point_mass function has been moved.
  54. Import it from "sympy.physics.mechanics".
  55. """,
  56. deprecated_since_version="1.13",
  57. active_deprecations_target="moved-mechanics-functions"
  58. )
  59. return _inertia_of_point_mass(mass, pos_vec, frame)
  60. def linear_momentum(frame, *body):
  61. """Linear momentum of the system.
  62. Explanation
  63. ===========
  64. This function returns the linear momentum of a system of Particle's and/or
  65. RigidBody's. The linear momentum of a system is equal to the vector sum of
  66. the linear momentum of its constituents. Consider a system, S, comprised of
  67. a rigid body, A, and a particle, P. The linear momentum of the system, L,
  68. is equal to the vector sum of the linear momentum of the particle, L1, and
  69. the linear momentum of the rigid body, L2, i.e.
  70. L = L1 + L2
  71. Parameters
  72. ==========
  73. frame : ReferenceFrame
  74. The frame in which linear momentum is desired.
  75. body1, body2, body3... : Particle and/or RigidBody
  76. The body (or bodies) whose linear momentum is required.
  77. Examples
  78. ========
  79. >>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
  80. >>> from sympy.physics.mechanics import RigidBody, outer, linear_momentum
  81. >>> N = ReferenceFrame('N')
  82. >>> P = Point('P')
  83. >>> P.set_vel(N, 10 * N.x)
  84. >>> Pa = Particle('Pa', P, 1)
  85. >>> Ac = Point('Ac')
  86. >>> Ac.set_vel(N, 25 * N.y)
  87. >>> I = outer(N.x, N.x)
  88. >>> A = RigidBody('A', Ac, N, 20, (I, Ac))
  89. >>> linear_momentum(N, A, Pa)
  90. 10*N.x + 500*N.y
  91. """
  92. if not isinstance(frame, ReferenceFrame):
  93. raise TypeError('Please specify a valid ReferenceFrame')
  94. else:
  95. linear_momentum_sys = Vector(0)
  96. for e in body:
  97. if isinstance(e, (RigidBody, Particle)):
  98. linear_momentum_sys += e.linear_momentum(frame)
  99. else:
  100. raise TypeError('*body must have only Particle or RigidBody')
  101. return linear_momentum_sys
  102. def angular_momentum(point, frame, *body):
  103. """Angular momentum of a system.
  104. Explanation
  105. ===========
  106. This function returns the angular momentum of a system of Particle's and/or
  107. RigidBody's. The angular momentum of such a system is equal to the vector
  108. sum of the angular momentum of its constituents. Consider a system, S,
  109. comprised of a rigid body, A, and a particle, P. The angular momentum of
  110. the system, H, is equal to the vector sum of the angular momentum of the
  111. particle, H1, and the angular momentum of the rigid body, H2, i.e.
  112. H = H1 + H2
  113. Parameters
  114. ==========
  115. point : Point
  116. The point about which angular momentum of the system is desired.
  117. frame : ReferenceFrame
  118. The frame in which angular momentum is desired.
  119. body1, body2, body3... : Particle and/or RigidBody
  120. The body (or bodies) whose angular momentum is required.
  121. Examples
  122. ========
  123. >>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
  124. >>> from sympy.physics.mechanics import RigidBody, outer, angular_momentum
  125. >>> N = ReferenceFrame('N')
  126. >>> O = Point('O')
  127. >>> O.set_vel(N, 0 * N.x)
  128. >>> P = O.locatenew('P', 1 * N.x)
  129. >>> P.set_vel(N, 10 * N.x)
  130. >>> Pa = Particle('Pa', P, 1)
  131. >>> Ac = O.locatenew('Ac', 2 * N.y)
  132. >>> Ac.set_vel(N, 5 * N.y)
  133. >>> a = ReferenceFrame('a')
  134. >>> a.set_ang_vel(N, 10 * N.z)
  135. >>> I = outer(N.z, N.z)
  136. >>> A = RigidBody('A', Ac, a, 20, (I, Ac))
  137. >>> angular_momentum(O, N, Pa, A)
  138. 10*N.z
  139. """
  140. if not isinstance(frame, ReferenceFrame):
  141. raise TypeError('Please enter a valid ReferenceFrame')
  142. if not isinstance(point, Point):
  143. raise TypeError('Please specify a valid Point')
  144. else:
  145. angular_momentum_sys = Vector(0)
  146. for e in body:
  147. if isinstance(e, (RigidBody, Particle)):
  148. angular_momentum_sys += e.angular_momentum(point, frame)
  149. else:
  150. raise TypeError('*body must have only Particle or RigidBody')
  151. return angular_momentum_sys
  152. def kinetic_energy(frame, *body):
  153. """Kinetic energy of a multibody system.
  154. Explanation
  155. ===========
  156. This function returns the kinetic energy of a system of Particle's and/or
  157. RigidBody's. The kinetic energy of such a system is equal to the sum of
  158. the kinetic energies of its constituents. Consider a system, S, comprising
  159. a rigid body, A, and a particle, P. The kinetic energy of the system, T,
  160. is equal to the vector sum of the kinetic energy of the particle, T1, and
  161. the kinetic energy of the rigid body, T2, i.e.
  162. T = T1 + T2
  163. Kinetic energy is a scalar.
  164. Parameters
  165. ==========
  166. frame : ReferenceFrame
  167. The frame in which the velocity or angular velocity of the body is
  168. defined.
  169. body1, body2, body3... : Particle and/or RigidBody
  170. The body (or bodies) whose kinetic energy is required.
  171. Examples
  172. ========
  173. >>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
  174. >>> from sympy.physics.mechanics import RigidBody, outer, kinetic_energy
  175. >>> N = ReferenceFrame('N')
  176. >>> O = Point('O')
  177. >>> O.set_vel(N, 0 * N.x)
  178. >>> P = O.locatenew('P', 1 * N.x)
  179. >>> P.set_vel(N, 10 * N.x)
  180. >>> Pa = Particle('Pa', P, 1)
  181. >>> Ac = O.locatenew('Ac', 2 * N.y)
  182. >>> Ac.set_vel(N, 5 * N.y)
  183. >>> a = ReferenceFrame('a')
  184. >>> a.set_ang_vel(N, 10 * N.z)
  185. >>> I = outer(N.z, N.z)
  186. >>> A = RigidBody('A', Ac, a, 20, (I, Ac))
  187. >>> kinetic_energy(N, Pa, A)
  188. 350
  189. """
  190. if not isinstance(frame, ReferenceFrame):
  191. raise TypeError('Please enter a valid ReferenceFrame')
  192. ke_sys = S.Zero
  193. for e in body:
  194. if isinstance(e, (RigidBody, Particle)):
  195. ke_sys += e.kinetic_energy(frame)
  196. else:
  197. raise TypeError('*body must have only Particle or RigidBody')
  198. return ke_sys
  199. def potential_energy(*body):
  200. """Potential energy of a multibody system.
  201. Explanation
  202. ===========
  203. This function returns the potential energy of a system of Particle's and/or
  204. RigidBody's. The potential energy of such a system is equal to the sum of
  205. the potential energy of its constituents. Consider a system, S, comprising
  206. a rigid body, A, and a particle, P. The potential energy of the system, V,
  207. is equal to the vector sum of the potential energy of the particle, V1, and
  208. the potential energy of the rigid body, V2, i.e.
  209. V = V1 + V2
  210. Potential energy is a scalar.
  211. Parameters
  212. ==========
  213. body1, body2, body3... : Particle and/or RigidBody
  214. The body (or bodies) whose potential energy is required.
  215. Examples
  216. ========
  217. >>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
  218. >>> from sympy.physics.mechanics import RigidBody, outer, potential_energy
  219. >>> from sympy import symbols
  220. >>> M, m, g, h = symbols('M m g h')
  221. >>> N = ReferenceFrame('N')
  222. >>> O = Point('O')
  223. >>> O.set_vel(N, 0 * N.x)
  224. >>> P = O.locatenew('P', 1 * N.x)
  225. >>> Pa = Particle('Pa', P, m)
  226. >>> Ac = O.locatenew('Ac', 2 * N.y)
  227. >>> a = ReferenceFrame('a')
  228. >>> I = outer(N.z, N.z)
  229. >>> A = RigidBody('A', Ac, a, M, (I, Ac))
  230. >>> Pa.potential_energy = m * g * h
  231. >>> A.potential_energy = M * g * h
  232. >>> potential_energy(Pa, A)
  233. M*g*h + g*h*m
  234. """
  235. pe_sys = S.Zero
  236. for e in body:
  237. if isinstance(e, (RigidBody, Particle)):
  238. pe_sys += e.potential_energy
  239. else:
  240. raise TypeError('*body must have only Particle or RigidBody')
  241. return pe_sys
  242. def gravity(acceleration, *bodies):
  243. from sympy.physics.mechanics.loads import gravity as _gravity
  244. sympy_deprecation_warning(
  245. """
  246. The gravity function has been moved.
  247. Import it from "sympy.physics.mechanics.loads".
  248. """,
  249. deprecated_since_version="1.13",
  250. active_deprecations_target="moved-mechanics-functions"
  251. )
  252. return _gravity(acceleration, *bodies)
  253. def center_of_mass(point, *bodies):
  254. """
  255. Returns the position vector from the given point to the center of mass
  256. of the given bodies(particles or rigidbodies).
  257. Example
  258. =======
  259. >>> from sympy import symbols, S
  260. >>> from sympy.physics.vector import Point
  261. >>> from sympy.physics.mechanics import Particle, ReferenceFrame, RigidBody, outer
  262. >>> from sympy.physics.mechanics.functions import center_of_mass
  263. >>> a = ReferenceFrame('a')
  264. >>> m = symbols('m', real=True)
  265. >>> p1 = Particle('p1', Point('p1_pt'), S(1))
  266. >>> p2 = Particle('p2', Point('p2_pt'), S(2))
  267. >>> p3 = Particle('p3', Point('p3_pt'), S(3))
  268. >>> p4 = Particle('p4', Point('p4_pt'), m)
  269. >>> b_f = ReferenceFrame('b_f')
  270. >>> b_cm = Point('b_cm')
  271. >>> mb = symbols('mb')
  272. >>> b = RigidBody('b', b_cm, b_f, mb, (outer(b_f.x, b_f.x), b_cm))
  273. >>> p2.point.set_pos(p1.point, a.x)
  274. >>> p3.point.set_pos(p1.point, a.x + a.y)
  275. >>> p4.point.set_pos(p1.point, a.y)
  276. >>> b.masscenter.set_pos(p1.point, a.y + a.z)
  277. >>> point_o=Point('o')
  278. >>> point_o.set_pos(p1.point, center_of_mass(p1.point, p1, p2, p3, p4, b))
  279. >>> expr = 5/(m + mb + 6)*a.x + (m + mb + 3)/(m + mb + 6)*a.y + mb/(m + mb + 6)*a.z
  280. >>> point_o.pos_from(p1.point)
  281. 5/(m + mb + 6)*a.x + (m + mb + 3)/(m + mb + 6)*a.y + mb/(m + mb + 6)*a.z
  282. """
  283. if not bodies:
  284. raise TypeError("No bodies(instances of Particle or Rigidbody) were passed.")
  285. total_mass = 0
  286. vec = Vector(0)
  287. for i in bodies:
  288. total_mass += i.mass
  289. masscenter = getattr(i, 'masscenter', None)
  290. if masscenter is None:
  291. masscenter = i.point
  292. vec += i.mass*masscenter.pos_from(point)
  293. return vec/total_mass
  294. def Lagrangian(frame, *body):
  295. """Lagrangian of a multibody system.
  296. Explanation
  297. ===========
  298. This function returns the Lagrangian of a system of Particle's and/or
  299. RigidBody's. The Lagrangian of such a system is equal to the difference
  300. between the kinetic energies and potential energies of its constituents. If
  301. T and V are the kinetic and potential energies of a system then it's
  302. Lagrangian, L, is defined as
  303. L = T - V
  304. The Lagrangian is a scalar.
  305. Parameters
  306. ==========
  307. frame : ReferenceFrame
  308. The frame in which the velocity or angular velocity of the body is
  309. defined to determine the kinetic energy.
  310. body1, body2, body3... : Particle and/or RigidBody
  311. The body (or bodies) whose Lagrangian is required.
  312. Examples
  313. ========
  314. >>> from sympy.physics.mechanics import Point, Particle, ReferenceFrame
  315. >>> from sympy.physics.mechanics import RigidBody, outer, Lagrangian
  316. >>> from sympy import symbols
  317. >>> M, m, g, h = symbols('M m g h')
  318. >>> N = ReferenceFrame('N')
  319. >>> O = Point('O')
  320. >>> O.set_vel(N, 0 * N.x)
  321. >>> P = O.locatenew('P', 1 * N.x)
  322. >>> P.set_vel(N, 10 * N.x)
  323. >>> Pa = Particle('Pa', P, 1)
  324. >>> Ac = O.locatenew('Ac', 2 * N.y)
  325. >>> Ac.set_vel(N, 5 * N.y)
  326. >>> a = ReferenceFrame('a')
  327. >>> a.set_ang_vel(N, 10 * N.z)
  328. >>> I = outer(N.z, N.z)
  329. >>> A = RigidBody('A', Ac, a, 20, (I, Ac))
  330. >>> Pa.potential_energy = m * g * h
  331. >>> A.potential_energy = M * g * h
  332. >>> Lagrangian(N, Pa, A)
  333. -M*g*h - g*h*m + 350
  334. """
  335. if not isinstance(frame, ReferenceFrame):
  336. raise TypeError('Please supply a valid ReferenceFrame')
  337. for e in body:
  338. if not isinstance(e, (RigidBody, Particle)):
  339. raise TypeError('*body must have only Particle or RigidBody')
  340. return kinetic_energy(frame, *body) - potential_energy(*body)
  341. def find_dynamicsymbols(expression, exclude=None, reference_frame=None):
  342. """Find all dynamicsymbols in expression.
  343. Explanation
  344. ===========
  345. If the optional ``exclude`` kwarg is used, only dynamicsymbols
  346. not in the iterable ``exclude`` are returned.
  347. If we intend to apply this function on a vector, the optional
  348. ``reference_frame`` is also used to inform about the corresponding frame
  349. with respect to which the dynamic symbols of the given vector is to be
  350. determined.
  351. Parameters
  352. ==========
  353. expression : SymPy expression
  354. exclude : iterable of dynamicsymbols, optional
  355. reference_frame : ReferenceFrame, optional
  356. The frame with respect to which the dynamic symbols of the
  357. given vector is to be determined.
  358. Examples
  359. ========
  360. >>> from sympy.physics.mechanics import dynamicsymbols, find_dynamicsymbols
  361. >>> from sympy.physics.mechanics import ReferenceFrame
  362. >>> x, y = dynamicsymbols('x, y')
  363. >>> expr = x + x.diff()*y
  364. >>> find_dynamicsymbols(expr)
  365. {x(t), y(t), Derivative(x(t), t)}
  366. >>> find_dynamicsymbols(expr, exclude=[x, y])
  367. {Derivative(x(t), t)}
  368. >>> a, b, c = dynamicsymbols('a, b, c')
  369. >>> A = ReferenceFrame('A')
  370. >>> v = a * A.x + b * A.y + c * A.z
  371. >>> find_dynamicsymbols(v, reference_frame=A)
  372. {a(t), b(t), c(t)}
  373. """
  374. t_set = {dynamicsymbols._t}
  375. if exclude:
  376. if iterable(exclude):
  377. exclude_set = set(exclude)
  378. else:
  379. raise TypeError("exclude kwarg must be iterable")
  380. else:
  381. exclude_set = set()
  382. if isinstance(expression, Vector):
  383. if reference_frame is None:
  384. raise ValueError("You must provide reference_frame when passing a "
  385. "vector expression, got %s." % reference_frame)
  386. else:
  387. expression = expression.to_matrix(reference_frame)
  388. return {i for i in expression.atoms(AppliedUndef, Derivative) if
  389. i.free_symbols == t_set} - exclude_set
  390. def msubs(expr, *sub_dicts, smart=False, **kwargs):
  391. """A custom subs for use on expressions derived in physics.mechanics.
  392. Traverses the expression tree once, performing the subs found in sub_dicts.
  393. Terms inside ``Derivative`` expressions are ignored:
  394. Examples
  395. ========
  396. >>> from sympy.physics.mechanics import dynamicsymbols, msubs
  397. >>> x = dynamicsymbols('x')
  398. >>> msubs(x.diff() + x, {x: 1})
  399. Derivative(x(t), t) + 1
  400. Note that sub_dicts can be a single dictionary, or several dictionaries:
  401. >>> x, y, z = dynamicsymbols('x, y, z')
  402. >>> sub1 = {x: 1, y: 2}
  403. >>> sub2 = {z: 3, x.diff(): 4}
  404. >>> msubs(x.diff() + x + y + z, sub1, sub2)
  405. 10
  406. If smart=True (default False), also checks for conditions that may result
  407. in ``nan``, but if simplified would yield a valid expression. For example:
  408. >>> from sympy import sin, tan
  409. >>> (sin(x)/tan(x)).subs(x, 0)
  410. nan
  411. >>> msubs(sin(x)/tan(x), {x: 0}, smart=True)
  412. 1
  413. It does this by first replacing all ``tan`` with ``sin/cos``. Then each
  414. node is traversed. If the node is a fraction, subs is first evaluated on
  415. the denominator. If this results in 0, simplification of the entire
  416. fraction is attempted. Using this selective simplification, only
  417. subexpressions that result in 1/0 are targeted, resulting in faster
  418. performance.
  419. """
  420. sub_dict = dict_merge(*sub_dicts)
  421. if smart:
  422. func = _smart_subs
  423. elif hasattr(expr, 'msubs'):
  424. return expr.msubs(sub_dict)
  425. else:
  426. func = lambda expr, sub_dict: _crawl(expr, _sub_func, sub_dict)
  427. if isinstance(expr, (Matrix, Vector, Dyadic)):
  428. return expr.applyfunc(lambda x: func(x, sub_dict))
  429. else:
  430. return func(expr, sub_dict)
  431. def _crawl(expr, func, *args, **kwargs):
  432. """Crawl the expression tree, and apply func to every node."""
  433. val = func(expr, *args, **kwargs)
  434. if val is not None:
  435. return val
  436. new_args = (_crawl(arg, func, *args, **kwargs) for arg in expr.args)
  437. return expr.func(*new_args)
  438. def _sub_func(expr, sub_dict):
  439. """Perform direct matching substitution, ignoring derivatives."""
  440. if expr in sub_dict:
  441. return sub_dict[expr]
  442. elif not expr.args or expr.is_Derivative:
  443. return expr
  444. def _tan_repl_func(expr):
  445. """Replace tan with sin/cos."""
  446. if isinstance(expr, tan):
  447. return sin(*expr.args) / cos(*expr.args)
  448. elif not expr.args or expr.is_Derivative:
  449. return expr
  450. def _smart_subs(expr, sub_dict):
  451. """Performs subs, checking for conditions that may result in `nan` or
  452. `oo`, and attempts to simplify them out.
  453. The expression tree is traversed twice, and the following steps are
  454. performed on each expression node:
  455. - First traverse:
  456. Replace all `tan` with `sin/cos`.
  457. - Second traverse:
  458. If node is a fraction, check if the denominator evaluates to 0.
  459. If so, attempt to simplify it out. Then if node is in sub_dict,
  460. sub in the corresponding value.
  461. """
  462. expr = _crawl(expr, _tan_repl_func)
  463. def _recurser(expr, sub_dict):
  464. # Decompose the expression into num, den
  465. num, den = _fraction_decomp(expr)
  466. if den != 1:
  467. # If there is a non trivial denominator, we need to handle it
  468. denom_subbed = _recurser(den, sub_dict)
  469. if denom_subbed.evalf() == 0:
  470. # If denom is 0 after this, attempt to simplify the bad expr
  471. expr = simplify(expr)
  472. else:
  473. # Expression won't result in nan, find numerator
  474. num_subbed = _recurser(num, sub_dict)
  475. return num_subbed / denom_subbed
  476. # We have to crawl the tree manually, because `expr` may have been
  477. # modified in the simplify step. First, perform subs as normal:
  478. val = _sub_func(expr, sub_dict)
  479. if val is not None:
  480. return val
  481. new_args = (_recurser(arg, sub_dict) for arg in expr.args)
  482. return expr.func(*new_args)
  483. return _recurser(expr, sub_dict)
  484. def _fraction_decomp(expr):
  485. """Return num, den such that expr = num/den."""
  486. if not isinstance(expr, Mul):
  487. return expr, 1
  488. num = []
  489. den = []
  490. for a in expr.args:
  491. if a.is_Pow and a.args[1] < 0:
  492. den.append(1 / a)
  493. else:
  494. num.append(a)
  495. if not den:
  496. return expr, 1
  497. num = Mul(*num)
  498. den = Mul(*den)
  499. return num, den
  500. def _f_list_parser(fl, ref_frame):
  501. """Parses the provided forcelist composed of items
  502. of the form (obj, force).
  503. Returns a tuple containing:
  504. vel_list: The velocity (ang_vel for Frames, vel for Points) in
  505. the provided reference frame.
  506. f_list: The forces.
  507. Used internally in the KanesMethod and LagrangesMethod classes.
  508. """
  509. def flist_iter():
  510. for pair in fl:
  511. obj, force = pair
  512. if isinstance(obj, ReferenceFrame):
  513. yield obj.ang_vel_in(ref_frame), force
  514. elif isinstance(obj, Point):
  515. yield obj.vel(ref_frame), force
  516. else:
  517. raise TypeError('First entry in each forcelist pair must '
  518. 'be a point or frame.')
  519. if not fl:
  520. vel_list, f_list = (), ()
  521. else:
  522. unzip = lambda l: list(zip(*l)) if l[0] else [(), ()]
  523. vel_list, f_list = unzip(list(flist_iter()))
  524. return vel_list, f_list
  525. def _validate_coordinates(coordinates=None, speeds=None, check_duplicates=True,
  526. is_dynamicsymbols=True, u_auxiliary=None):
  527. """Validate the generalized coordinates and generalized speeds.
  528. Parameters
  529. ==========
  530. coordinates : iterable, optional
  531. Generalized coordinates to be validated.
  532. speeds : iterable, optional
  533. Generalized speeds to be validated.
  534. check_duplicates : bool, optional
  535. Checks if there are duplicates in the generalized coordinates and
  536. generalized speeds. If so it will raise a ValueError. The default is
  537. True.
  538. is_dynamicsymbols : iterable, optional
  539. Checks if all the generalized coordinates and generalized speeds are
  540. dynamicsymbols. If any is not a dynamicsymbol, a ValueError will be
  541. raised. The default is True.
  542. u_auxiliary : iterable, optional
  543. Auxiliary generalized speeds to be validated.
  544. """
  545. t_set = {dynamicsymbols._t}
  546. # Convert input to iterables
  547. if coordinates is None:
  548. coordinates = []
  549. elif not iterable(coordinates):
  550. coordinates = [coordinates]
  551. if speeds is None:
  552. speeds = []
  553. elif not iterable(speeds):
  554. speeds = [speeds]
  555. if u_auxiliary is None:
  556. u_auxiliary = []
  557. elif not iterable(u_auxiliary):
  558. u_auxiliary = [u_auxiliary]
  559. msgs = []
  560. if check_duplicates: # Check for duplicates
  561. seen = set()
  562. coord_duplicates = {x for x in coordinates if x in seen or seen.add(x)}
  563. seen = set()
  564. speed_duplicates = {x for x in speeds if x in seen or seen.add(x)}
  565. seen = set()
  566. aux_duplicates = {x for x in u_auxiliary if x in seen or seen.add(x)}
  567. overlap_coords = set(coordinates).intersection(speeds)
  568. overlap_aux = set(coordinates).union(speeds).intersection(u_auxiliary)
  569. if coord_duplicates:
  570. msgs.append(f'The generalized coordinates {coord_duplicates} are '
  571. f'duplicated, all generalized coordinates should be '
  572. f'unique.')
  573. if speed_duplicates:
  574. msgs.append(f'The generalized speeds {speed_duplicates} are '
  575. f'duplicated, all generalized speeds should be unique.')
  576. if aux_duplicates:
  577. msgs.append(f'The auxiliary speeds {aux_duplicates} are duplicated,'
  578. f' all auxiliary speeds should be unique.')
  579. if overlap_coords:
  580. msgs.append(f'{overlap_coords} are defined as both generalized '
  581. f'coordinates and generalized speeds.')
  582. if overlap_aux:
  583. msgs.append(f'The auxiliary speeds {overlap_aux} are also defined '
  584. f'as generalized coordinates or generalized speeds.')
  585. if is_dynamicsymbols: # Check whether all coordinates are dynamicsymbols
  586. for coordinate in coordinates:
  587. if not (isinstance(coordinate, (AppliedUndef, Derivative)) and
  588. coordinate.free_symbols == t_set):
  589. msgs.append(f'Generalized coordinate "{coordinate}" is not a '
  590. f'dynamicsymbol.')
  591. for speed in speeds:
  592. if not (isinstance(speed, (AppliedUndef, Derivative)) and
  593. speed.free_symbols == t_set):
  594. msgs.append(
  595. f'Generalized speed "{speed}" is not a dynamicsymbol.')
  596. for aux in u_auxiliary:
  597. if not (isinstance(aux, (AppliedUndef, Derivative)) and
  598. aux.free_symbols == t_set):
  599. msgs.append(
  600. f'Auxiliary speed "{aux}" is not a dynamicsymbol.')
  601. if msgs:
  602. raise ValueError('\n'.join(msgs))
  603. def _parse_linear_solver(linear_solver):
  604. """Helper function to retrieve a specified linear solver."""
  605. if callable(linear_solver):
  606. return linear_solver
  607. return lambda A, b: Matrix.solve(A, b, method=linear_solver)