test_vector.py 9.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247
  1. from sympy.core.numbers import (Float, pi)
  2. from sympy.core.symbol import symbols
  3. from sympy.core.sorting import ordered
  4. from sympy.functions.elementary.trigonometric import (cos, sin)
  5. from sympy.matrices.immutable import ImmutableDenseMatrix as Matrix
  6. from sympy.physics.vector import ReferenceFrame, Vector, dynamicsymbols, dot
  7. from sympy.physics.vector.vector import VectorTypeError
  8. from sympy.abc import x, y, z
  9. from sympy.testing.pytest import raises
  10. Vector.simp = True
  11. A = ReferenceFrame('A')
  12. def test_free_dynamicsymbols():
  13. A, B, C, D = symbols('A, B, C, D', cls=ReferenceFrame)
  14. a, b, c, d, e, f = dynamicsymbols('a, b, c, d, e, f')
  15. B.orient_axis(A, a, A.x)
  16. C.orient_axis(B, b, B.y)
  17. D.orient_axis(C, c, C.x)
  18. v = d*D.x + e*D.y + f*D.z
  19. assert set(ordered(v.free_dynamicsymbols(A))) == {a, b, c, d, e, f}
  20. assert set(ordered(v.free_dynamicsymbols(B))) == {b, c, d, e, f}
  21. assert set(ordered(v.free_dynamicsymbols(C))) == {c, d, e, f}
  22. assert set(ordered(v.free_dynamicsymbols(D))) == {d, e, f}
  23. def test_Vector():
  24. assert A.x != A.y
  25. assert A.y != A.z
  26. assert A.z != A.x
  27. assert A.x + 0 == A.x
  28. v1 = x*A.x + y*A.y + z*A.z
  29. v2 = x**2*A.x + y**2*A.y + z**2*A.z
  30. v3 = v1 + v2
  31. v4 = v1 - v2
  32. assert isinstance(v1, Vector)
  33. assert dot(v1, A.x) == x
  34. assert dot(v1, A.y) == y
  35. assert dot(v1, A.z) == z
  36. assert isinstance(v2, Vector)
  37. assert dot(v2, A.x) == x**2
  38. assert dot(v2, A.y) == y**2
  39. assert dot(v2, A.z) == z**2
  40. assert isinstance(v3, Vector)
  41. # We probably shouldn't be using simplify in dot...
  42. assert dot(v3, A.x) == x**2 + x
  43. assert dot(v3, A.y) == y**2 + y
  44. assert dot(v3, A.z) == z**2 + z
  45. assert isinstance(v4, Vector)
  46. # We probably shouldn't be using simplify in dot...
  47. assert dot(v4, A.x) == x - x**2
  48. assert dot(v4, A.y) == y - y**2
  49. assert dot(v4, A.z) == z - z**2
  50. assert v1.to_matrix(A) == Matrix([[x], [y], [z]])
  51. q = symbols('q')
  52. B = A.orientnew('B', 'Axis', (q, A.x))
  53. assert v1.to_matrix(B) == Matrix([[x],
  54. [ y * cos(q) + z * sin(q)],
  55. [-y * sin(q) + z * cos(q)]])
  56. #Test the separate method
  57. B = ReferenceFrame('B')
  58. v5 = x*A.x + y*A.y + z*B.z
  59. assert Vector(0).separate() == {}
  60. assert v1.separate() == {A: v1}
  61. assert v5.separate() == {A: x*A.x + y*A.y, B: z*B.z}
  62. #Test the free_symbols property
  63. v6 = x*A.x + y*A.y + z*A.z
  64. assert v6.free_symbols(A) == {x,y,z}
  65. raises(TypeError, lambda: v3.applyfunc(v1))
  66. def test_Vector_diffs():
  67. q1, q2, q3, q4 = dynamicsymbols('q1 q2 q3 q4')
  68. q1d, q2d, q3d, q4d = dynamicsymbols('q1 q2 q3 q4', 1)
  69. q1dd, q2dd, q3dd, q4dd = dynamicsymbols('q1 q2 q3 q4', 2)
  70. N = ReferenceFrame('N')
  71. A = N.orientnew('A', 'Axis', [q3, N.z])
  72. B = A.orientnew('B', 'Axis', [q2, A.x])
  73. v1 = q2 * A.x + q3 * N.y
  74. v2 = q3 * B.x + v1
  75. v3 = v1.dt(B)
  76. v4 = v2.dt(B)
  77. v5 = q1*A.x + q2*A.y + q3*A.z
  78. assert v1.dt(N) == q2d * A.x + q2 * q3d * A.y + q3d * N.y
  79. assert v1.dt(A) == q2d * A.x + q3 * q3d * N.x + q3d * N.y
  80. assert v1.dt(B) == (q2d * A.x + q3 * q3d * N.x + q3d *
  81. N.y - q3 * cos(q3) * q2d * N.z)
  82. assert v2.dt(N) == (q2d * A.x + (q2 + q3) * q3d * A.y + q3d * B.x + q3d *
  83. N.y)
  84. assert v2.dt(A) == q2d * A.x + q3d * B.x + q3 * q3d * N.x + q3d * N.y
  85. assert v2.dt(B) == (q2d * A.x + q3d * B.x + q3 * q3d * N.x + q3d * N.y -
  86. q3 * cos(q3) * q2d * N.z)
  87. assert v3.dt(N) == (q2dd * A.x + q2d * q3d * A.y + (q3d**2 + q3 * q3dd) *
  88. N.x + q3dd * N.y + (q3 * sin(q3) * q2d * q3d -
  89. cos(q3) * q2d * q3d - q3 * cos(q3) * q2dd) * N.z)
  90. assert v3.dt(A) == (q2dd * A.x + (2 * q3d**2 + q3 * q3dd) * N.x + (q3dd -
  91. q3 * q3d**2) * N.y + (q3 * sin(q3) * q2d * q3d -
  92. cos(q3) * q2d * q3d - q3 * cos(q3) * q2dd) * N.z)
  93. assert v3.dt(B) == (q2dd * A.x - q3 * cos(q3) * q2d**2 * A.y + (2 *
  94. q3d**2 + q3 * q3dd) * N.x + (q3dd - q3 * q3d**2) *
  95. N.y + (2 * q3 * sin(q3) * q2d * q3d - 2 * cos(q3) *
  96. q2d * q3d - q3 * cos(q3) * q2dd) * N.z)
  97. assert v4.dt(N) == (q2dd * A.x + q3d * (q2d + q3d) * A.y + q3dd * B.x +
  98. (q3d**2 + q3 * q3dd) * N.x + q3dd * N.y + (q3 *
  99. sin(q3) * q2d * q3d - cos(q3) * q2d * q3d - q3 *
  100. cos(q3) * q2dd) * N.z)
  101. assert v4.dt(A) == (q2dd * A.x + q3dd * B.x + (2 * q3d**2 + q3 * q3dd) *
  102. N.x + (q3dd - q3 * q3d**2) * N.y + (q3 * sin(q3) *
  103. q2d * q3d - cos(q3) * q2d * q3d - q3 * cos(q3) *
  104. q2dd) * N.z)
  105. assert v4.dt(B) == (q2dd * A.x - q3 * cos(q3) * q2d**2 * A.y + q3dd * B.x +
  106. (2 * q3d**2 + q3 * q3dd) * N.x + (q3dd - q3 * q3d**2) *
  107. N.y + (2 * q3 * sin(q3) * q2d * q3d - 2 * cos(q3) *
  108. q2d * q3d - q3 * cos(q3) * q2dd) * N.z)
  109. assert v5.dt(B) == q1d*A.x + (q3*q2d + q2d)*A.y + (-q2*q2d + q3d)*A.z
  110. assert v5.dt(A) == q1d*A.x + q2d*A.y + q3d*A.z
  111. assert v5.dt(N) == (-q2*q3d + q1d)*A.x + (q1*q3d + q2d)*A.y + q3d*A.z
  112. assert v3.diff(q1d, N) == 0
  113. assert v3.diff(q2d, N) == A.x - q3 * cos(q3) * N.z
  114. assert v3.diff(q3d, N) == q3 * N.x + N.y
  115. assert v3.diff(q1d, A) == 0
  116. assert v3.diff(q2d, A) == A.x - q3 * cos(q3) * N.z
  117. assert v3.diff(q3d, A) == q3 * N.x + N.y
  118. assert v3.diff(q1d, B) == 0
  119. assert v3.diff(q2d, B) == A.x - q3 * cos(q3) * N.z
  120. assert v3.diff(q3d, B) == q3 * N.x + N.y
  121. assert v4.diff(q1d, N) == 0
  122. assert v4.diff(q2d, N) == A.x - q3 * cos(q3) * N.z
  123. assert v4.diff(q3d, N) == B.x + q3 * N.x + N.y
  124. assert v4.diff(q1d, A) == 0
  125. assert v4.diff(q2d, A) == A.x - q3 * cos(q3) * N.z
  126. assert v4.diff(q3d, A) == B.x + q3 * N.x + N.y
  127. assert v4.diff(q1d, B) == 0
  128. assert v4.diff(q2d, B) == A.x - q3 * cos(q3) * N.z
  129. assert v4.diff(q3d, B) == B.x + q3 * N.x + N.y
  130. # diff() should only express vector components in the derivative frame if
  131. # the orientation of the component's frame depends on the variable
  132. v6 = q2**2*N.y + q2**2*A.y + q2**2*B.y
  133. # already expressed in N
  134. n_measy = 2*q2
  135. # A_C_N does not depend on q2, so don't express in N
  136. a_measy = 2*q2
  137. # B_C_N depends on q2, so express in N
  138. b_measx = (q2**2*B.y).dot(N.x).diff(q2)
  139. b_measy = (q2**2*B.y).dot(N.y).diff(q2)
  140. b_measz = (q2**2*B.y).dot(N.z).diff(q2)
  141. n_comp, a_comp = v6.diff(q2, N).args
  142. assert len(v6.diff(q2, N).args) == 2 # only N and A parts
  143. assert n_comp[1] == N
  144. assert a_comp[1] == A
  145. assert n_comp[0] == Matrix([b_measx, b_measy + n_measy, b_measz])
  146. assert a_comp[0] == Matrix([0, a_measy, 0])
  147. def test_vector_var_in_dcm():
  148. N = ReferenceFrame('N')
  149. A = ReferenceFrame('A')
  150. B = ReferenceFrame('B')
  151. u1, u2, u3, u4 = dynamicsymbols('u1 u2 u3 u4')
  152. v = u1 * u2 * A.x + u3 * N.y + u4**2 * N.z
  153. assert v.diff(u1, N, var_in_dcm=False) == u2 * A.x
  154. assert v.diff(u1, A, var_in_dcm=False) == u2 * A.x
  155. assert v.diff(u3, N, var_in_dcm=False) == N.y
  156. assert v.diff(u3, A, var_in_dcm=False) == N.y
  157. assert v.diff(u3, B, var_in_dcm=False) == N.y
  158. assert v.diff(u4, N, var_in_dcm=False) == 2 * u4 * N.z
  159. raises(ValueError, lambda: v.diff(u1, N))
  160. def test_vector_simplify():
  161. x, y, z, k, n, m, w, f, s, A = symbols('x, y, z, k, n, m, w, f, s, A')
  162. N = ReferenceFrame('N')
  163. test1 = (1 / x + 1 / y) * N.x
  164. assert (test1 & N.x) != (x + y) / (x * y)
  165. test1 = test1.simplify()
  166. assert (test1 & N.x) == (x + y) / (x * y)
  167. test2 = (A**2 * s**4 / (4 * pi * k * m**3)) * N.x
  168. test2 = test2.simplify()
  169. assert (test2 & N.x) == (A**2 * s**4 / (4 * pi * k * m**3))
  170. test3 = ((4 + 4 * x - 2 * (2 + 2 * x)) / (2 + 2 * x)) * N.x
  171. test3 = test3.simplify()
  172. assert (test3 & N.x) == 0
  173. test4 = ((-4 * x * y**2 - 2 * y**3 - 2 * x**2 * y) / (x + y)**2) * N.x
  174. test4 = test4.simplify()
  175. assert (test4 & N.x) == -2 * y
  176. def test_vector_evalf():
  177. a, b = symbols('a b')
  178. v = pi * A.x
  179. assert v.evalf(2) == Float('3.1416', 2) * A.x
  180. v = pi * A.x + 5 * a * A.y - b * A.z
  181. assert v.evalf(3) == Float('3.1416', 3) * A.x + Float('5', 3) * a * A.y - b * A.z
  182. assert v.evalf(5, subs={a: 1.234, b:5.8973}) == Float('3.1415926536', 5) * A.x + Float('6.17', 5) * A.y - Float('5.8973', 5) * A.z
  183. def test_vector_angle():
  184. A = ReferenceFrame('A')
  185. v1 = A.x + A.y
  186. v2 = A.z
  187. assert v1.angle_between(v2) == pi/2
  188. B = ReferenceFrame('B')
  189. B.orient_axis(A, A.x, pi)
  190. v3 = A.x
  191. v4 = B.x
  192. assert v3.angle_between(v4) == 0
  193. def test_vector_xreplace():
  194. x, y, z = symbols('x y z')
  195. v = x**2 * A.x + x*y * A.y + x*y*z * A.z
  196. assert v.xreplace({x : cos(x)}) == cos(x)**2 * A.x + y*cos(x) * A.y + y*z*cos(x) * A.z
  197. assert v.xreplace({x*y : pi}) == x**2 * A.x + pi * A.y + x*y*z * A.z
  198. assert v.xreplace({x*y*z : 1}) == x**2*A.x + x*y*A.y + A.z
  199. assert v.xreplace({x:1, z:0}) == A.x + y * A.y
  200. raises(TypeError, lambda: v.xreplace())
  201. raises(TypeError, lambda: v.xreplace([x, y]))
  202. def test_issue_23366():
  203. u1 = dynamicsymbols('u1')
  204. N = ReferenceFrame('N')
  205. N_v_A = u1*N.x
  206. raises(VectorTypeError, lambda: N_v_A.diff(N, u1))