Repository: sympy/sympy

Issue:
qdots in KanesMethod.forcing vector
When running the following example, which forms the equations of motion of a rolling disc, it can be seen that there are $\dot{\bar{q}}$'s in the forcing vector:
```py
from sympy import symbols
from sympy.physics.mechanics import (dynamicsymbols, cross, inertia, RigidBody,
                                     ReferenceFrame, KanesMethod,
                                     find_dynamicsymbols)

t = dynamicsymbols._t
q1, q2, q3, q4, q5, u1, u2, u3, u4, u5 = dynamicsymbols('q1:6 u1:6')
g, r = symbols('g r')
ground = RigidBody('ground')
disc = RigidBody('disk')
disc.inertia = (disc.mass * r ** 2 / 4 * inertia(disc.frame, 1, 2, 1),
                disc.masscenter)
int_frame = ReferenceFrame('int_frame')
int_frame.orient_body_fixed(ground.frame, (q1, q2, 0), 'zxy')
disc.frame.orient_axis(int_frame, int_frame.y, q3)
g_w_d = disc.frame.ang_vel_in(ground.frame)
disc.frame.set_ang_vel(ground.frame,
                       u1 * disc.x + u2 * disc.y + u3 * disc.z)
cp = ground.masscenter.locatenew('contact_point',
                                 q4 * ground.x + q5 * ground.y)
cp.set_vel(ground.frame, u4 * ground.x + u5 * ground.y)
disc.masscenter.set_pos(cp, r * int_frame.z)
loads = [(disc.masscenter, -disc.mass * g * ground.z)]
v_cp = cp.vel(ground.frame) + cross(
    disc.frame.ang_vel_in(ground.frame), cp.pos_from(disc.masscenter))
fnh = [v_cp.dot(ground.x), v_cp.dot(ground.y)]
kdes = [g_w_d.dot(disc.x) - u1, g_w_d.dot(disc.y) - u2,
        g_w_d.dot(disc.z) - u3, q4.diff(t) - u4, q5.diff(t) - u5]

kane = KanesMethod(ground.frame, [q1, q2, q3, q4, q5], [u1, u2, u3], kdes,
                   u_dependent=[u4, u5], velocity_constraints=fnh,
                   bodies=[disc], forcelist=loads, explicit_kinematics=False)
kane.kanes_equations()
print(find_dynamicsymbols(kane.forcing))
# {q1(t), u1(t), u3(t), Derivative(q2(t), t), u2(t), q3(t), Derivative(q3(t), t), q2(t)}
```
I tracked this problem down to the [`qddot_to_u_map`](https://github.com/sympy/sympy/blob/2bf714de7c8c9e8d62211e781a3139943f3f6546/sympy/physics/mechanics/kane.py#L378) having $\dot{\bar{q}}$'s in its values. This can be easily fixed by using `q_ddot_u_map = {k.diff(t): v.diff(t).xreplace(self._qdot_u_map) for (k, v) in self._qdot_u_map.items()}` instead.

Discussion:

Patch:
diff --git a/sympy/physics/mechanics/kane.py b/sympy/physics/mechanics/kane.py
index c0197553f8..3404244440 100644
--- a/sympy/physics/mechanics/kane.py
+++ b/sympy/physics/mechanics/kane.py
@@ -375,8 +375,8 @@ def _form_frstar(self, bl):
         uauxdot = [diff(i, t) for i in self._uaux]
         uauxdot_zero = {i: 0 for i in uauxdot}
         # Dictionary of q' and q'' to u and u'
-        q_ddot_u_map = {k.diff(t): v.diff(t) for (k, v) in
-                self._qdot_u_map.items()}
+        q_ddot_u_map = {k.diff(t): v.diff(t).xreplace(
+            self._qdot_u_map) for (k, v) in self._qdot_u_map.items()}
         q_ddot_u_map.update(self._qdot_u_map)
 
         # Fill up the list of partials: format is a list with num elements


Test patch:
diff --git a/sympy/physics/mechanics/tests/test_kane.py b/sympy/physics/mechanics/tests/test_kane.py
index ea8ed28894..b166e7fd6f 100644
--- a/sympy/physics/mechanics/tests/test_kane.py
+++ b/sympy/physics/mechanics/tests/test_kane.py
@@ -4,7 +4,7 @@
 from sympy.simplify.simplify import simplify
 from sympy.physics.mechanics import (dynamicsymbols, ReferenceFrame, Point,
                                      RigidBody, KanesMethod, inertia, Particle,
-                                     dot)
+                                     dot, find_dynamicsymbols)
 from sympy.testing.pytest import raises
 from sympy.core.backend import USE_SYMENGINE
 
@@ -530,3 +530,29 @@ def test_implicit_kinematics():
     lambda_0_sol = solve(config_cons[0], q_att_vec[0])[1]
     lhs_candidate = simplify(mass_matrix_kin_implicit * qdot_candidate).subs({q_att_vec[0]: lambda_0_sol})
     assert lhs_candidate == forcing_kin_implicit
+
+def test_issue_24887():
+    # Spherical pendulum
+    g, l, m, c = symbols('g l m c')
+    q1, q2, q3, u1, u2, u3 = dynamicsymbols('q1:4 u1:4')
+    N = ReferenceFrame('N')
+    A = ReferenceFrame('A')
+    A.orient_body_fixed(N, (q1, q2, q3), 'zxy')
+    N_w_A = A.ang_vel_in(N)
+    # A.set_ang_vel(N, u1 * A.x + u2 * A.y + u3 * A.z)
+    kdes = [N_w_A.dot(A.x) - u1, N_w_A.dot(A.y) - u2, N_w_A.dot(A.z) - u3]
+    O = Point('O')
+    O.set_vel(N, 0)
+    Po = O.locatenew('Po', -l * A.y)
+    Po.set_vel(A, 0)
+    P = Particle('P', Po, m)
+    kane = KanesMethod(N, [q1, q2, q3], [u1, u2, u3], kdes, bodies=[P],
+                       forcelist=[(Po, -m * g * N.y)])
+    kane.kanes_equations()
+    expected_md = m * l ** 2 * Matrix([[1, 0, 0], [0, 0, 0], [0, 0, 1]])
+    expected_fd = Matrix([
+        [l*m*(g*(sin(q1)*sin(q3) - sin(q2)*cos(q1)*cos(q3)) - l*u2*u3)],
+        [0], [l*m*(-g*(sin(q1)*cos(q3) + sin(q2)*sin(q3)*cos(q1)) + l*u1*u2)]])
+    assert find_dynamicsymbols(kane.forcing).issubset({q1, q2, q3, u1, u2, u3})
+    assert simplify(kane.mass_matrix - expected_md) == zeros(3, 3)
+    assert simplify(kane.forcing - expected_fd) == zeros(3, 1)