Repository: sympy/sympy

Issue:
Incorrect velocity calculation
As already pointed out in #22956 an incorrect assumption seems to be added in [21730](https://github.com/sympy/sympy/pull/21730).

Here is an example:
```py
import sympy as sm
import sympy.physics.mechanics as me

q1, q2 = me.dynamicsymbols('q1:3')
u1, u2 = sm.symbols('u1:3')
A = me.ReferenceFrame('A')
B = me.ReferenceFrame('B')
C = me.ReferenceFrame('C')
B.orient_axis(A, A.z, q1)
B.set_ang_vel(A, u1*A.z)
C.orient_axis(B, B.z, q2)
C.set_ang_vel(B, u2*B.z)

Am = me.Point('Am')
Am.set_vel(A, 0)
Bm = me.Point('Bm')
Bm.set_pos(Am, B.x)
Bm.set_vel(B, 0)
Bm.set_vel(C, 0)
Cm = me.Point('Cm')
Cm.set_pos(Bm, C.x)
Cm.set_vel(C, 0)
temp = Cm._vel_dict.copy()
print(Cm.vel(A))  # Takes the position derivative
Cm._vel_dict = temp
Cm.v2pt_theory(Bm, B, C)
print(Cm.vel(A))  # Uses the velocity from an intermediate frame
```
This outputs:
```py
u1*B.y + (u1 + u2)*C.y
u2*C.y
```
Which should be the same, but as shown it isn't. In fact with the code implemented in [21730](https://github.com/sympy/sympy/pull/21730) the incorrect assumption is made that:
<img src="https://render.githubusercontent.com/render/math?math={{}^Av^{Cm} = {}^Nv^{Am} %2B {}^Bv^{Cm}}#gh-light-mode-only">
<img src="https://render.githubusercontent.com/render/math?math={\color{white}{}^Av^{Cm} = {}^Nv^{Am} %2B {}^Bv^{Cm}}#gh-dark-mode-only">
Because one the of components of the position vector from `Am` to `Cm` is expressed in the `B` frame.

@sidhu1012 or @moorepants could you elaborate why this assumption could normally be made? Also the test that was added, namely [`test_vel_frame`](https://github.com/sympy/sympy/blob/4c969e4ded0ab76a1a22fa3c93bc8f22af854e22/sympy/physics/vector/tests/test_point.py#L302) is also not really correct seen from a mathematical/mechanical point of view. Since the `N` frame is equal to the `A` frame the distance between `P` and `P1` is constant, so the velocity of `P` should be the same as `P1`.

PS As always correct me if I'm wrong, I am still learning more about mechanics.

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v163). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* physics.vector
  * Fixed bug in the calculation of the velocity of a point. ([#23362](https://github.com/sympy/sympy/pull/23362) by [@TJStienstra](https://github.com/TJStienstra) and [@Tstienstra](https://github.com/Tstienstra))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.11.

<details><summary>Click here to see the pull request description that was parsed.</summary>

    This reverts commit 4c969e4ded0ab76a1a22fa3c93bc8f22af854e22.

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->


    #### Brief description of what is fixed or changed
    Fixes #23358 by removing an incorrect assumption in the velocity calculation of a point.

    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below between the BEGIN and END
    statements. The basic format is a bulleted list with the name of the subpackage
    and the release note for this PR. For example:

    * solvers
      * Added a new solver for logarithmic equations.

    * functions
      * Fixed a bug with log of integers.

    or if no release note(s) should be included use:

    NO ENTRY

    See https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more
    information on how to write release notes. The bot will check your release
    notes automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->
    * physics.vector
      * Fixed bug in the calculation of the velocity of a point.
    <!-- END RELEASE NOTES -->


</details><p>

moorepants: Are you sure it was correct before @sidhu1012's PR?

You should add a test case based on your example you showed here to prove it works.
moorepants: Yes

Patch:
diff --git a/.mailmap b/.mailmap
index 1b21777276..ea7b3be7b0 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1149,6 +1149,9 @@ Tiffany Zhu <bubble.wubble.303@gmail.com>
 Tim Gates <tim.gates@iress.com>
 Tim Lahey <tim.lahey@gmail.com>
 Tim Swast <tswast@gmail.com>
+Timo Stienstra <timostienstra00@gmail.com> TJStienstra <97806294+TJStienstra@users.noreply.github.com>
+Timo Stienstra <timostienstra00@gmail.com> TJStienstra <T.J.Stienstra@student.tudelft.nl>
+Timo Stienstra <timostienstra00@gmail.com> TJStienstra <timostienstra00@gmail.com>
 Timothy Cyrus <tcyrus@users.noreply.github.com>
 Timothy Reluga <treluga@math.psu.edu> <tcr-public2@psu.edu>
 Timothy Reluga <treluga@math.psu.edu> me <none@none>
diff --git a/doc/src/modules/physics/mechanics/examples/multi_degree_freedom_holonomic_system.rst b/doc/src/modules/physics/mechanics/examples/multi_degree_freedom_holonomic_system.rst
index f9a8b21673..eef16e5920 100644
--- a/doc/src/modules/physics/mechanics/examples/multi_degree_freedom_holonomic_system.rst
+++ b/doc/src/modules/physics/mechanics/examples/multi_degree_freedom_holonomic_system.rst
@@ -22,7 +22,7 @@ mass of the pendulum, as well as gravity and others. ::
 
     >>> from sympy import zeros, symbols
     >>> from sympy.physics.mechanics import Body, PinJoint, PrismaticJoint, JointsMethod, inertia
-    >>> from sympy.physics.mechanics import dynamicsymbols 
+    >>> from sympy.physics.mechanics import dynamicsymbols
     >>> q1, q2, q3, u1, u2, u3 = dynamicsymbols('q1, q2, q3, u1, u2, u3')
     >>> l, k, c, g, kT = symbols('l, k, c, g, kT')
     >>> ma, mb, mc, IBzz= symbols('ma, mb, mc, IBzz')
@@ -57,7 +57,7 @@ specified torque T acts on compound_pend and block, specified force F acts on bl
     >>> block.apply_force(-k*q1*block.x, reaction_body=wall)
     >>> block.apply_force(-c*u1*block.x, reaction_body=wall)
     >>> compound_pend.apply_torque(T*compound_pend.z, reaction_body=block)
-    >>> simple_pend.apply_torque(kT*q3*simple_pend.z, reaction_body=compound_pend)
+    >>> simple_pend.apply_torque(-kT*q3*simple_pend.z, reaction_body=compound_pend)
     >>> block.apply_force(-wall.y*block.mass*g)
     >>> compound_pend.apply_force(-wall.y*compound_pend.mass*g)
     >>> simple_pend.apply_force(-wall.y*simple_pend.mass*g)
@@ -68,24 +68,24 @@ With the problem setup, the equations of motion can be generated using the
     >>> method = JointsMethod(wall, slider, rev1, rev2)
     >>> method.form_eoms()
     Matrix([
-    [               -c*u1(t) - k*q1(t) + 2*l*mb*u2(t)**2*sin(q2(t))/3 - l*mc*(-sin(q2(t))*cos(q3(t)) - sin(q3(t))*cos(q2(t)))*(u2(t) + u3(t))*u3(t) - mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))/3)*Derivative(u3(t), t) + mc*(2*l*u2(t)/3 + l*u3(t)/3)*u2(t)*sin(q2(t)) - (2*l*mb*cos(q2(t))/3 + 2*l*mc*cos(q2(t))/3)*Derivative(u2(t), t) - (ma + mb + mc)*Derivative(u1(t), t) + F(t)],
-    [                                                                                                                  -2*g*l*mb*sin(q2(t))/3 - 2*g*l*mc*sin(q2(t))/3 + 2*l**2*mc*(u2(t) + u3(t))*u3(t)*sin(q3(t))/3 - mc*(2*l**2*cos(q3(t))/3 + 2*l**2/9)*Derivative(u3(t), t) - (2*l*mb*cos(q2(t))/3 + 2*l*mc*cos(q2(t))/3)*Derivative(u1(t), t) - (IBzz + 4*l**2*mb/9 + 4*l**2*mc/9)*Derivative(u2(t), t) + T(t)],
-    [-g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - g*l*mc*sin(q2(t))/3 + kT*q3(t) + l**2*mc*(u2(t) + u3(t))*u3(t)*sin(q3(t))/3 - l*mc*(2*l*u2(t)/3 + l*u3(t)/3)*u2(t)*sin(q3(t)) - mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))/3)*Derivative(u1(t), t) - mc*(2*l**2*cos(q3(t))/3 + 2*l**2/9)*Derivative(u2(t), t) - mc*(2*l**2*cos(q3(t))/3 + 10*l**2/9)*Derivative(u3(t), t)]])
+    [                                -c*u1(t) - k*q1(t) + 2*l*mb*u2(t)**2*sin(q2(t))/3 - l*mc*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t)))*Derivative(u3(t), t) - l*mc*(-sin(q2(t))*cos(q3(t)) - sin(q3(t))*cos(q2(t)))*(u2(t) + u3(t))**2 + l*mc*u2(t)**2*sin(q2(t)) - (2*l*mb*cos(q2(t))/3 + mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))))*Derivative(u2(t), t) - (ma + mb + mc)*Derivative(u1(t), t) + F(t)],
+    [-2*g*l*mb*sin(q2(t))/3 - g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - g*l*mc*sin(q2(t)) + l**2*mc*(u2(t) + u3(t))**2*sin(q3(t)) - l**2*mc*u2(t)**2*sin(q3(t)) - mc*(l**2*cos(q3(t)) + l**2)*Derivative(u3(t), t) - (2*l*mb*cos(q2(t))/3 + mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))))*Derivative(u1(t), t) - (IBzz + 4*l**2*mb/9 + mc*(2*l**2*cos(q3(t)) + 2*l**2))*Derivative(u2(t), t) + T(t)],
+    [                                                                                                                                                                        -g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - kT*q3(t) - l**2*mc*u2(t)**2*sin(q3(t)) - l**2*mc*Derivative(u3(t), t) - l*mc*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t)))*Derivative(u1(t), t) - mc*(l**2*cos(q3(t)) + l**2)*Derivative(u2(t), t)]])
 
     >>> method.mass_matrix_full
     Matrix([
-    [1, 0, 0,                                                                        0,                                         0,                                                                        0],
-    [0, 1, 0,                                                                        0,                                         0,                                                                        0],
-    [0, 0, 1,                                                                        0,                                         0,                                                                        0],
-    [0, 0, 0,                                                             ma + mb + mc, 2*l*mb*cos(q2(t))/3 + 2*l*mc*cos(q2(t))/3, mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))/3)],
-    [0, 0, 0,                                2*l*mb*cos(q2(t))/3 + 2*l*mc*cos(q2(t))/3,          IBzz + 4*l**2*mb/9 + 4*l**2*mc/9,                                      mc*(2*l**2*cos(q3(t))/3 + 2*l**2/9)],
-    [0, 0, 0, mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))/3),       mc*(2*l**2*cos(q3(t))/3 + 2*l**2/9),                                     mc*(2*l**2*cos(q3(t))/3 + 10*l**2/9)]])
+    [1, 0, 0,                                                                                            0,                                                                                            0,                                                     0],
+    [0, 1, 0,                                                                                            0,                                                                                            0,                                                     0],
+    [0, 0, 1,                                                                                            0,                                                                                            0,                                                     0],
+    [0, 0, 0,                                                                                 ma + mb + mc, 2*l*mb*cos(q2(t))/3 + mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))), l*mc*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t)))],
+    [0, 0, 0, 2*l*mb*cos(q2(t))/3 + mc*(l*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))) + l*cos(q2(t))),                                         IBzz + 4*l**2*mb/9 + mc*(2*l**2*cos(q3(t)) + 2*l**2),                           mc*(l**2*cos(q3(t)) + l**2)],
+    [0, 0, 0,                                        l*mc*(-sin(q2(t))*sin(q3(t)) + cos(q2(t))*cos(q3(t))),                                                                  mc*(l**2*cos(q3(t)) + l**2),                                               l**2*mc]])
 
     >>> method.forcing_full
     Matrix([
-    [                                                                                                                                                                                  u1(t)],
-    [                                                                                                                                                                                  u2(t)],
-    [                                                                                                                                                                                  u3(t)],
-    [ -c*u1(t) - k*q1(t) + 2*l*mb*u2(t)**2*sin(q2(t))/3 - l*mc*(-sin(q2(t))*cos(q3(t)) - sin(q3(t))*cos(q2(t)))*(u2(t) + u3(t))*u3(t) + mc*(2*l*u2(t)/3 + l*u3(t)/3)*u2(t)*sin(q2(t)) + F(t)],
-    [                                                                                   -2*g*l*mb*sin(q2(t))/3 - 2*g*l*mc*sin(q2(t))/3 + 2*l**2*mc*(u2(t) + u3(t))*u3(t)*sin(q3(t))/3 + T(t)],
-    [-g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - g*l*mc*sin(q2(t))/3 + kT*q3(t) + l**2*mc*(u2(t) + u3(t))*u3(t)*sin(q3(t))/3 - l*mc*(2*l*u2(t)/3 + l*u3(t)/3)*u2(t)*sin(q3(t))]])
+    [                                                                                                                                                                           u1(t)],
+    [                                                                                                                                                                           u2(t)],
+    [                                                                                                                                                                           u3(t)],
+    [                  -c*u1(t) - k*q1(t) + 2*l*mb*u2(t)**2*sin(q2(t))/3 - l*mc*(-sin(q2(t))*cos(q3(t)) - sin(q3(t))*cos(q2(t)))*(u2(t) + u3(t))**2 + l*mc*u2(t)**2*sin(q2(t)) + F(t)],
+    [-2*g*l*mb*sin(q2(t))/3 - g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - g*l*mc*sin(q2(t)) + l**2*mc*(u2(t) + u3(t))**2*sin(q3(t)) - l**2*mc*u2(t)**2*sin(q3(t)) + T(t)],
+    [                                                                                -g*l*mc*(sin(q2(t))*cos(q3(t)) + sin(q3(t))*cos(q2(t))) - kT*q3(t) - l**2*mc*u2(t)**2*sin(q3(t))]])
diff --git a/sympy/physics/mechanics/joint.py b/sympy/physics/mechanics/joint.py
index 007e74b06c..b9e774479a 100644
--- a/sympy/physics/mechanics/joint.py
+++ b/sympy/physics/mechanics/joint.py
@@ -702,7 +702,7 @@ class PrismaticJoint(Joint):
     v_J1(t)*W_frame.x
 
     >>> Part2.masscenter.vel(wall.frame)
-    v_J1(t)*W_frame.x + v_J2(t)*P1_frame.x
+    v_J1(t)*W_frame.x + Derivative(x_J2(t), t)*P1_frame.x
 
     """
 
diff --git a/sympy/physics/vector/point.py b/sympy/physics/vector/point.py
index 3cfb9adc17..c54dd0c266 100644
--- a/sympy/physics/vector/point.py
+++ b/sympy/physics/vector/point.py
@@ -564,15 +564,7 @@ def vel(self, frame):
                             continue
                         candidate_neighbor.append(neighbor)
                         if not valid_neighbor_found:
-                            vel = None
-                            for f in self.pos_from(neighbor).args:
-                                if f[1] in self._vel_dict.keys():
-                                    if self._vel_dict[f[1]] != 0:
-                                        vel = self._vel_dict[f[1]]
-                                        break
-                            if vel is None:
-                                vel = self.pos_from(neighbor).dt(frame)
-                            self.set_vel(frame, vel + neighbor_velocity)
+                            self.set_vel(frame, self.pos_from(neighbor).dt(frame) + neighbor_velocity)
                             valid_neighbor_found = True
             if is_cyclic:
                 warn('Kinematic loops are defined among the positions of '


Test patch:
diff --git a/sympy/physics/vector/tests/test_point.py b/sympy/physics/vector/tests/test_point.py
index ce38cbdd00..a978f33146 100644
--- a/sympy/physics/vector/tests/test_point.py
+++ b/sympy/physics/vector/tests/test_point.py
@@ -299,13 +299,28 @@ def test_auto_vel_multiple_path_warning_msg():
         assert 'Velocity automatically calculated based on point' in str(w[-1].message)
         assert 'Velocities from these points are not necessarily the same. This may cause errors in your calculations.' in str(w[-1].message)
 
-def test_vel_frame():
-    N = ReferenceFrame('N')
+def test_auto_vel_derivative():
+    q1, q2 = dynamicsymbols('q1:3')
+    u1, u2 = dynamicsymbols('u1:3', 1)
     A = ReferenceFrame('A')
-    A.orient_axis(N, N.x, 0)
-    P = Point('P')
-    P.set_vel(N, 5*N.x)
-    P1 = Point('P1')
-    P1.set_pos(P, A.x+N.x)
-    P1.set_vel(A, 10*A.x)
-    assert P1.vel(N) == 5*N.x + 10*A.x
+    B = ReferenceFrame('B')
+    C = ReferenceFrame('C')
+    B.orient_axis(A, A.z, q1)
+    B.set_ang_vel(A, u1 * A.z)
+    C.orient_axis(B, B.z, q2)
+    C.set_ang_vel(B, u2 * B.z)
+
+    Am = Point('Am')
+    Am.set_vel(A, 0)
+    Bm = Point('Bm')
+    Bm.set_pos(Am, B.x)
+    Bm.set_vel(B, 0)
+    Bm.set_vel(C, 0)
+    Cm = Point('Cm')
+    Cm.set_pos(Bm, C.x)
+    Cm.set_vel(C, 0)
+    temp = Cm._vel_dict.copy()
+    assert Cm.vel(A) == (u1 * B.y + (u1 + u2) * C.y)
+    Cm._vel_dict = temp
+    Cm.v2pt_theory(Bm, B, C)
+    assert Cm.vel(A) == (u1 * B.y + (u1 + u2) * C.y)