Repository: sympy/sympy

Issue:
test_pinjoint_arbitrary_axis fails if ran separately
When running `bin/test sympy/physics/mechanics` all tests pass, but if I run `bin/test sympy/physics/mechanics/tests/test_joint` the [test_pinjoint_arbitrary_axis](https://github.com/sympy/sympy/blob/b1f35f8daf5efa9f801f905cf1037000d69b0278/sympy/physics/mechanics/tests/test_joint.py#L153) fails.
```
sympy\physics\mechanics\tests\test_joint.py[8] ....F...                                      [FAIL]
sympy\physics\mechanics\tests\test_jointsmethod.py[7] .......                                  [OK]
___________________________________________________________________________________________________
___________ sympy\physics\mechanics\tests\test_joint.py:test_pinjoint_arbitrary_axis_______________
Traceback (most recent call last):
  File "...\sympy\physics\mechanics\tests\test_joint.py", line 223, in test_pinjoint_arbitrary_axis
    assert angle.xreplace({omega: 1}) == 0
AssertionError
````
This fail is due to the fact that when running this test separately the angle is not simplified (when done it is `0`). It is only simplified when testing the entire module.

The test fails for [line 223](https://github.com/sympy/sympy/blob/b1f35f8daf5efa9f801f905cf1037000d69b0278/sympy/physics/mechanics/tests/test_joint.py#L223), [line 250](https://github.com/sympy/sympy/blob/b1f35f8daf5efa9f801f905cf1037000d69b0278/sympy/physics/mechanics/tests/test_joint.py#L250) and [line 282](https://github.com/sympy/sympy/blob/b1f35f8daf5efa9f801f905cf1037000d69b0278/sympy/physics/mechanics/tests/test_joint.py#L282).

Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index 6f19ca644f..478444ff02 100644
--- a/.mailmap
+++ b/.mailmap
@@ -384,7 +384,6 @@ Bhavya Srivastava <bhavya17037@iiitd.ac.in>
 Bilal Ahmed <b.ahmed0918@gmail.com>
 Bilal Akhtar <bilalakhtar@ubuntu.com> <bilalakhtar96@yahoo.com>
 Bill Flynn <wflynny@gmail.com>
-Biswadeep Purkayastha <98874428+metabiswadeep@users.noreply.github.com>
 Björn Dahlgren <bjodah@gmail.com>
 Blair Azzopardi <bsdz@users.noreply.github.com>
 Boris Atamanovskiy <shaomoron@gmail.com>
diff --git a/doc/src/explanation/active-deprecations.md b/doc/src/explanation/active-deprecations.md
index 51dfaac6c7..4c6d753019 100644
--- a/doc/src/explanation/active-deprecations.md
+++ b/doc/src/explanation/active-deprecations.md
@@ -76,11 +76,22 @@ SymPy deprecation warnings.
 
 ## Version 1.12
 
+(deprecated-mechanics-joint-coordinate-format)=
+### New Joint coordinate format
+The format, i.e. type and auto generated name, of the generalized coordinates
+and generalized speeds of the joints in the ``sympy.physics.mechanics`` module
+has changed. The data type has changed from ``list`` to ``Matrix``, which is the
+same as the type for the generalized coordinates within the ``KanesMethod``.
+The auto naming of the generalized coordinates and generalized speeds of the
+``PinJoint`` and ``PrismaticJoint`` have also changed to ``q_<joint.name>`` and
+``u_<joint.name>``. Previously each of those joints had an unique template for
+auto generating these names.
+
 (deprecated-mechanics-joint-axis)=
 ### New Joint intermediate frames
 
-The definition of the joint axis in the ``sympy.physics.mechanics`` has changed.
-Instead of using the arguments ``parent_axis`` and ``child_axis`` to
+The definition of the joint axis in the ``sympy.physics.mechanics`` module has
+changed. Instead of using the arguments ``parent_axis`` and ``child_axis`` to
 automatically determine the joint axis and an intermediate reference frame, the
 joints now use an intermediate frame argument for both the parent and the child
 body, i.e. ``parent_interframe`` and ``child_interframe``. This means that you
@@ -102,9 +113,9 @@ joint was:
 ...                child_axis=-child.z)   # doctest: +SKIP
 >>> parent.dcm(child)   # doctest: +SKIP
 Matrix([
-[-cos(theta_pin(t)), -sin(theta_pin(t)),  0],
-[-sin(theta_pin(t)),  cos(theta_pin(t)),  0],
-[                 0,                  0, -1]])
+[-cos(q_pin(t)), -sin(q_pin(t)),  0],
+[-sin(q_pin(t)),  cos(q_pin(t)),  0],
+[             0,              0, -1]])
 ```
 
 When inspecting this matrix you will notice that for ``theta_pin = 0`` the child
@@ -122,9 +133,9 @@ this exact rotation:
 ...                child_interframe=int_frame)
 >>> parent.dcm(child)
 Matrix([
-[-cos(theta_pin(t)), -sin(theta_pin(t)),  0],
-[-sin(theta_pin(t)),  cos(theta_pin(t)),  0],
-[                 0,                  0, -1]])
+[-cos(q_pin(t)), -sin(q_pin(t)),  0],
+[-sin(q_pin(t)),  cos(q_pin(t)),  0],
+[             0,              0, -1]])
 ```
 
 However if you liked the fact that the deprecated arguments aligned the frames
@@ -140,9 +151,9 @@ given vector:
 ...                child_interframe=-child.z)
 >>> parent.dcm(child)
 Matrix([
-[-cos(theta_pin(t)), -sin(theta_pin(t)),  0],
-[-sin(theta_pin(t)),  cos(theta_pin(t)),  0],
-[                 0,                  0, -1]])
+[-cos(q_pin(t)), -sin(q_pin(t)),  0],
+[-sin(q_pin(t)),  cos(q_pin(t)),  0],
+[             0,              0, -1]])
 ```
 
 (deprecated-mechanics-joint-pos)=
diff --git a/doc/src/modules/physics/mechanics/api/PinJoint.svg b/doc/src/modules/physics/mechanics/api/PinJoint.svg
index 34e0fda824..18fd68db66 100644
--- a/doc/src/modules/physics/mechanics/api/PinJoint.svg
+++ b/doc/src/modules/physics/mechanics/api/PinJoint.svg
@@ -29,7 +29,7 @@
         <text fill="#FF0000"
               font-family="Calibri,Calibri_MSFontService,sans-serif"
               font-weight="400" font-size="83"
-              transform="matrix(1 0 0 1 3259.58 1576)">θ, ω<tspan font-size="83"
+              transform="matrix(1 0 0 1 3259.58 1576)">q, u<tspan font-size="83"
                                                                   x="-42.585"
                                                                   y="395">joint_axis</tspan>
             <tspan fill="#008000" font-size="83" x="-1704.28" y="-33">parent.masscenter</tspan></text>
diff --git a/doc/src/modules/physics/mechanics/api/PrismaticJoint.svg b/doc/src/modules/physics/mechanics/api/PrismaticJoint.svg
index 5db1394228..cc0e1f8ac6 100644
--- a/doc/src/modules/physics/mechanics/api/PrismaticJoint.svg
+++ b/doc/src/modules/physics/mechanics/api/PrismaticJoint.svg
@@ -47,7 +47,7 @@
         <text fill="#FF0000"
               font-family="Calibri,Calibri_MSFontService,sans-serif"
               font-weight="400" font-size="64"
-              transform="matrix(1 0 0 1 2447.49 1370)">q1, u1<tspan
+              transform="matrix(1 0 0 1 2447.49 1370)">q, u<tspan
                 fill="#008000" font-size="64" x="-1301.45" y="2">parent.masscenter</tspan>
             <tspan fill="#0070C0" font-size="64" x="-132.069" y="-529">child.masscenter</tspan>
             <tspan fill="#008000" font-size="64" x="-1085.66" y="-232">parent_point</tspan>
diff --git a/doc/src/modules/physics/mechanics/joints.rst b/doc/src/modules/physics/mechanics/joints.rst
index b6a9a350ac..fea8644a47 100644
--- a/doc/src/modules/physics/mechanics/joints.rst
+++ b/doc/src/modules/physics/mechanics/joints.rst
@@ -46,8 +46,8 @@ the joint axis is defined in the ``parent_interframe`` as $2\hat{p}_x +
 4\hat{p}_y + 3\hat{p}_z$, then this will also be $2\hat{c}_x + 4\hat{c}_y +
 3\hat{c}_z$ in the ``child_interframe``. Practically this means that in the case
 of the :class:`~.PinJoint`, also shown below, the ``joint_axis`` is the axis of
-rotation, with the generalized coordinate :math:`\theta` as the angle of
-rotation and the generalized speed :math:`\omega` as the angular velocity.
+rotation, with the generalized coordinate :math:`q` as the angle of
+rotation and the generalized speed :math:`u` as the angular velocity.
 
 .. image:: api/PinJoint.svg
    :align: center
@@ -67,16 +67,16 @@ body's frame. ::
 
    >>> from sympy.physics.mechanics import *
    >>> mechanics_printing(pretty_print=False)
-   >>> theta, omega = dynamicsymbols('theta, omega')
+   >>> q, u = dynamicsymbols('q, u')
    >>> parent = Body('parent')
    >>> child = Body('child')
    >>> joint = PinJoint(
-   ...     'hinge', parent, child, coordinates=theta, speeds=omega,
+   ...     'hinge', parent, child, coordinates=q, speeds=u,
    ...     parent_point=3 * parent.frame.x,
    ...     child_point=-3 * child.frame.x,
    ...     joint_axis=parent.frame.z)
    >>> joint.kdes
-   [omega - theta']
+   Matrix([[u - q']])
    >>> joint.parent_point.pos_from(parent.masscenter)
    3*parent_frame.x
    >>> joint.parent_interframe
@@ -86,7 +86,7 @@ body's frame. ::
    >>> child.masscenter.pos_from(parent.masscenter)
    3*parent_frame.x + 3*child_frame.x
    >>> child.masscenter.vel(parent.frame)
-   3*omega*child_frame.y
+   3*u*child_frame.y
 
 JointsMethod in Physics/Mechanics
 =================================
@@ -101,6 +101,6 @@ In the code below we form the equations of motion of the single
 
    >>> method = JointsMethod(parent, joint)
    >>> method.form_eoms()
-   Matrix([[-(child_izz + 9*child_mass)*omega']])
+   Matrix([[-(child_izz + 9*child_mass)*u']])
    >>> type(method.method)  # The method working in the backend
    <class 'sympy.physics.mechanics.kane.KanesMethod'>
diff --git a/doc/src/modules/solvers/inequalities.rst b/doc/src/modules/solvers/inequalities.rst
index 18bc0bbd7a..2498217e0a 100644
--- a/doc/src/modules/solvers/inequalities.rst
+++ b/doc/src/modules/solvers/inequalities.rst
@@ -3,9 +3,6 @@
 Inequality Solvers
 ==================
 
-For general cases :func:`~.reduce_inequalities` should be used.
-Other functions are the subcategories useful for special dedicated operations, and will be called internally as needed by ``reduce_inequalities``.
-
 .. module:: sympy.solvers.inequalities
 
 .. autofunction:: solve_rational_inequalities
diff --git a/sympy/functions/elementary/integers.py b/sympy/functions/elementary/integers.py
index 871f5ea95c..2f53c40a1a 100644
--- a/sympy/functions/elementary/integers.py
+++ b/sympy/functions/elementary/integers.py
@@ -150,7 +150,14 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             r = floor(arg0)
         if arg0.is_finite:
             if arg0 == r:
-                ndir = arg.dir(x, cdir=cdir)
+                if cdir == 0:
+                    ndirl = arg.dir(x, cdir=-1)
+                    ndir = arg.dir(x, cdir=1)
+                    if ndir != ndirl:
+                        raise ValueError("Two sided limit of %s around 0"
+                                    "does not exist" % self)
+                else:
+                    ndir = arg.dir(x, cdir=cdir)
                 return r - 1 if ndir.is_negative else r
             else:
                 return r
@@ -308,7 +315,14 @@ def _eval_as_leading_term(self, x, logx=None, cdir=0):
             r = ceiling(arg0)
         if arg0.is_finite:
             if arg0 == r:
-                ndir = arg.dir(x, cdir=cdir)
+                if cdir == 0:
+                    ndirl = arg.dir(x, cdir=-1)
+                    ndir = arg.dir(x, cdir=1)
+                    if ndir != ndirl:
+                        raise ValueError("Two sided limit of %s around 0"
+                                    "does not exist" % self)
+                else:
+                    ndir = arg.dir(x, cdir=cdir)
                 return r if ndir.is_negative else r + 1
             else:
                 return r
@@ -573,43 +587,6 @@ def _value_one_or_more(self, other):
             if other.is_integer and other.is_positive:
                 return S.true
 
-    def _eval_as_leading_term(self, x, logx=None, cdir=0):
-        from sympy.calculus.accumulationbounds import AccumBounds
-        arg = self.args[0]
-        arg0 = arg.subs(x, 0)
-        r = self.subs(x, 0)
-
-        if arg0.is_finite:
-            if r.is_zero:
-                ndir = arg.dir(x, cdir=cdir)
-                if ndir.is_negative:
-                    return S.One
-                return (arg - arg0).as_leading_term(x, logx=logx, cdir=cdir)
-            else:
-                return r
-        elif arg0 in (S.ComplexInfinity, S.Infinity, S.NegativeInfinity):
-            return AccumBounds(0, 1)
-        return arg.as_leading_term(x, logx=logx, cdir=cdir)
-
-    def _eval_nseries(self, x, n, logx, cdir=0):
-        from sympy.series.order import Order
-        arg = self.args[0]
-        arg0 = arg.subs(x, 0)
-        r = self.subs(x, 0)
-
-        if arg0.is_infinite:
-            from sympy.calculus.accumulationbounds import AccumBounds
-            o = Order(1, (x, 0)) if n <= 0 else AccumBounds(0, 1) + Order(x**n, (x, 0))
-            return o
-        else:
-            res = (arg - arg0)._eval_nseries(x, n, logx=logx, cdir=cdir)
-            if r.is_zero:
-                ndir = arg.dir(x, cdir=cdir)
-                res += S.One if ndir.is_negative else S.Zero
-            else:
-                res += r
-            return res
-
 
 @dispatch(frac, Basic)  # type:ignore
 def _eval_is_eq(lhs, rhs): # noqa:F811
diff --git a/sympy/physics/mechanics/joint.py b/sympy/physics/mechanics/joint.py
index 590c0e3f7f..2718c43f03 100644
--- a/sympy/physics/mechanics/joint.py
+++ b/sympy/physics/mechanics/joint.py
@@ -2,7 +2,7 @@
 
 from abc import ABC, abstractmethod
 
-from sympy.core.numbers import pi
+from sympy.core.backend import pi, AppliedUndef, Derivative, Matrix
 from sympy.physics.mechanics.body import Body
 from sympy.physics.vector import (Vector, dynamicsymbols, cross, Point,
                                   ReferenceFrame)
@@ -40,9 +40,9 @@ class Joint(ABC):
         The parent body of joint.
     child : Body
         The child body of joint.
-    coordinates: List of dynamicsymbols, optional
+    coordinates: iterable of dynamicsymbols, optional
         Generalized coordinates of the joint.
-    speeds : List of dynamicsymbols, optional
+    speeds : iterable of dynamicsymbols, optional
         Generalized speeds of joint.
     parent_point : Point or Vector, optional
         Attachment point where the joint is fixed to the parent body. If a
@@ -96,10 +96,10 @@ class Joint(ABC):
         The joint's parent body.
     child : Body
         The joint's child body.
-    coordinates : list
-        List of the joint's generalized coordinates.
-    speeds : list
-        List of the joint's generalized speeds.
+    coordinates : Matrix
+        Matrix of the joint's generalized coordinates.
+    speeds : Matrix
+        Matrix of the joint's generalized speeds.
     parent_point : Point
         Attachment point where the joint is fixed to the parent body.
     child_point : Point
@@ -114,7 +114,7 @@ class Joint(ABC):
     child_interframe : ReferenceFrame
         Intermediate frame of the child body with respect to which the joint
         transformation is formulated.
-    kdes : list
+    kdes : Matrix
         Kinematical differential equations of the joint.
 
     Notes
@@ -215,12 +215,12 @@ def child(self):
 
     @property
     def coordinates(self):
-        """List generalized coordinates of the joint."""
+        """Matrix of the joint's generalized coordinates."""
         return self._coordinates
 
     @property
     def speeds(self):
-        """List generalized coordinates of the joint.."""
+        """Matrix of the joint's generalized speeds."""
         return self._speeds
 
     @property
@@ -260,12 +260,12 @@ def child_interframe(self):
 
     @abstractmethod
     def _generate_coordinates(self, coordinates):
-        """Generate list generalized coordinates of the joint."""
+        """Generate Matrix of the joint's generalized coordinates."""
         pass
 
     @abstractmethod
     def _generate_speeds(self, speeds):
-        """Generate list generalized speeds of the joint."""
+        """Generate Matrix of the joint's generalized speeds."""
         pass
 
     @abstractmethod
@@ -442,7 +442,7 @@ def _generate_kdes(self):
         t = dynamicsymbols._t
         for i in range(len(self.coordinates)):
             kdes.append(-self.coordinates[i].diff(t) + self.speeds[i])
-        return kdes
+        return Matrix(kdes)
 
     def _locate_joint_pos(self, body, joint_pos):
         """Returns the attachment point of a body."""
@@ -474,6 +474,55 @@ def _locate_joint_frame(self, body, interframe):
         body.masscenter.set_vel(interframe, 0)  # Fixate interframe to body
         return interframe
 
+    def _fill_coordinate_list(self, coordinates, n_coords, label='q', offset=0,
+                              number_single=False):
+        """Helper method for _generate_coordinates and _generate_speeds.
+
+        Parameters
+        ==========
+
+        coordinates : iterable
+            Iterable of coordinates or speeds that have been provided.
+        n_coords : Integer
+            Number of coordinates that should be returned.
+        label : String, optional
+            Coordinate type either 'q' (coordinates) or 'u' (speeds). The
+            Default is 'q'.
+        offset : Integer
+            Count offset when creating new dynamicsymbols. The default is 0.
+        number_single : Boolean
+            Boolean whether if n_coords == 1, number should still be used. The
+            default is False.
+
+        """
+
+        def create_symbol(number):
+            if n_coords == 1 and not number_single:
+                return dynamicsymbols(f'{label}_{self.name}')
+            return dynamicsymbols(f'{label}{number}_{self.name}')
+
+        name = 'generalized coordinate' if label == 'q' else 'generalized speed'
+        generated_coordinates = []
+        if coordinates is None:
+            coordinates = []
+        elif isinstance(coordinates, (AppliedUndef, Derivative)):
+            coordinates = [coordinates]
+        # Supports more iterables, also Matrix
+        for i, coord in enumerate(coordinates):
+            if coord is None:
+                generated_coordinates.append(create_symbol(i + offset))
+            elif isinstance(coord, (AppliedUndef, Derivative)):
+                generated_coordinates.append(coord)
+            else:
+                raise TypeError(f'The {name} {coord} should have been a '
+                                f'dynamicsymbol.')
+        for i in range(len(coordinates) + offset, n_coords + offset):
+            generated_coordinates.append(create_symbol(i))
+        if len(generated_coordinates) != n_coords:
+            raise ValueError(f'{len(coordinates)} {name}s have been provided. '
+                             f'The maximum number of {name}s is {n_coords}.')
+        return Matrix(generated_coordinates)
+
 
 class _JointAxisMixin:
     """Mixin class, which provides joint axis support methods.
@@ -586,10 +635,12 @@ class PinJoint(Joint, _JointAxisMixin):
         The joint's parent body.
     child : Body
         The joint's child body.
-    coordinates : list
-        List of the joint's generalized coordinates.
-    speeds : list
-        List of the joint's generalized speeds.
+    coordinates : Matrix
+        Matrix of the joint's generalized coordinates. The default value is
+        ``dynamicsymbols(f'q_{joint.name}')``.
+    speeds : Matrix
+        Matrix of the joint's generalized speeds. The default value is
+        ``dynamicsymbols(f'u_{joint.name}')``.
     parent_point : Point
         Attachment point where the joint is fixed to the parent body.
     child_point : Point
@@ -607,7 +658,7 @@ class PinJoint(Joint, _JointAxisMixin):
     joint_axis : Vector
         The axis about which the rotation occurs. Note that the components of
         this axis are the same in the parent_interframe and child_interframe.
-    kdes : list
+    kdes : Matrix
         Kinematical differential equations of the joint.
 
     Examples
@@ -641,16 +692,16 @@ class PinJoint(Joint, _JointAxisMixin):
     >>> joint.child_axis
     C_frame.x
     >>> joint.coordinates
-    [theta_PC(t)]
+    Matrix([[q_PC(t)]])
     >>> joint.speeds
-    [omega_PC(t)]
+    Matrix([[u_PC(t)]])
     >>> joint.child.frame.ang_vel_in(joint.parent.frame)
-    omega_PC(t)*P_frame.x
+    u_PC(t)*P_frame.x
     >>> joint.child.frame.dcm(joint.parent.frame)
     Matrix([
-    [1,                 0,                0],
-    [0,  cos(theta_PC(t)), sin(theta_PC(t))],
-    [0, -sin(theta_PC(t)), cos(theta_PC(t))]])
+    [1,             0,            0],
+    [0,  cos(q_PC(t)), sin(q_PC(t))],
+    [0, -sin(q_PC(t)), cos(q_PC(t))]])
     >>> joint.child_point.pos_from(joint.parent_point)
     0
 
@@ -689,14 +740,14 @@ class PinJoint(Joint, _JointAxisMixin):
 
     >>> upper_bob.frame.dcm(ceiling.frame)
     Matrix([
-    [ cos(theta_P1(t)), sin(theta_P1(t)), 0],
-    [-sin(theta_P1(t)), cos(theta_P1(t)), 0],
-    [                0,                0, 1]])
+    [ cos(q_P1(t)), sin(q_P1(t)), 0],
+    [-sin(q_P1(t)), cos(q_P1(t)), 0],
+    [            0,            0, 1]])
     >>> trigsimp(lower_bob.frame.dcm(ceiling.frame))
     Matrix([
-    [ cos(theta_P1(t) + theta_P2(t)), sin(theta_P1(t) + theta_P2(t)), 0],
-    [-sin(theta_P1(t) + theta_P2(t)), cos(theta_P1(t) + theta_P2(t)), 0],
-    [                              0,                              0, 1]])
+    [ cos(q_P1(t) + q_P2(t)), sin(q_P1(t) + q_P2(t)), 0],
+    [-sin(q_P1(t) + q_P2(t)), cos(q_P1(t) + q_P2(t)), 0],
+    [                      0,                      0, 1]])
 
     The position of the lower bob's masscenter is found with:
 
@@ -707,17 +758,17 @@ class PinJoint(Joint, _JointAxisMixin):
     respect to the ceiling.
 
     >>> upper_bob.frame.ang_vel_in(ceiling.frame)
-    omega_P1(t)*C_frame.z
+    u_P1(t)*C_frame.z
     >>> lower_bob.frame.ang_vel_in(ceiling.frame)
-    omega_P1(t)*C_frame.z + omega_P2(t)*U_frame.z
+    u_P1(t)*C_frame.z + u_P2(t)*U_frame.z
 
     And finally, the linear velocities of the two pendulum bobs can be computed
     with respect to the ceiling.
 
     >>> upper_bob.masscenter.vel(ceiling.frame)
-    l1*omega_P1(t)*U_frame.y
+    l1*u_P1(t)*U_frame.y
     >>> lower_bob.masscenter.vel(ceiling.frame)
-    l1*omega_P1(t)*U_frame.y + l2*(omega_P1(t) + omega_P2(t))*L_frame.y
+    l1*u_P1(t)*U_frame.y + l2*(u_P1(t) + u_P2(t))*L_frame.y
 
     """
 
@@ -742,20 +793,10 @@ def joint_axis(self):
         return self._joint_axis
 
     def _generate_coordinates(self, coordinate):
-        coordinates = []
-        if coordinate is None:
-            theta = dynamicsymbols('theta' + '_' + self._name)
-            coordinate = theta
-        coordinates.append(coordinate)
-        return coordinates
+        return self._fill_coordinate_list(coordinate, 1, 'q')
 
     def _generate_speeds(self, speed):
-        speeds = []
-        if speed is None:
-            omega = dynamicsymbols('omega' + '_' + self._name)
-            speed = omega
-        speeds.append(speed)
-        return speeds
+        return self._fill_coordinate_list(speed, 1, 'u')
 
     def _orient_frames(self):
         self._joint_axis = self._axis(
@@ -803,9 +844,11 @@ class PrismaticJoint(Joint, _JointAxisMixin):
     child : Body
         The child body of joint.
     coordinates: dynamicsymbol, optional
-        Generalized coordinates of the joint.
+        Generalized coordinates of the joint. The default value is
+        ``dynamicsymbols(f'q_{joint.name}')``.
     speeds : dynamicsymbol, optional
-        Generalized speeds of joint.
+        Generalized speeds of joint. The default value is
+        ``dynamicsymbols(f'u_{joint.name}')``.
     parent_point : Point or Vector, optional
         Attachment point where the joint is fixed to the parent body. If a
         vector is provided, then the attachment point is computed by adding the
@@ -861,10 +904,10 @@ class PrismaticJoint(Joint, _JointAxisMixin):
         The joint's parent body.
     child : Body
         The joint's child body.
-    coordinates : list
-        List of the joint's generalized coordinates.
-    speeds : list
-        List of the joint's generalized speeds.
+    coordinates : Matrix
+        Matrix of the joint's generalized coordinates.
+    speeds : Matrix
+        Matrix of the joint's generalized speeds.
     parent_point : Point
         Attachment point where the joint is fixed to the parent body.
     child_point : Point
@@ -879,7 +922,7 @@ class PrismaticJoint(Joint, _JointAxisMixin):
     child_interframe : ReferenceFrame
         Intermediate frame of the child body with respect to which the joint
         transformation is formulated.
-    kdes : list
+    kdes : Matrix
         Kinematical differential equations of the joint.
 
     Examples
@@ -913,9 +956,9 @@ class PrismaticJoint(Joint, _JointAxisMixin):
     >>> joint.child_axis
     C_frame.x
     >>> joint.coordinates
-    [x_PC(t)]
+    Matrix([[q_PC(t)]])
     >>> joint.speeds
-    [v_PC(t)]
+    Matrix([[u_PC(t)]])
     >>> joint.child.frame.ang_vel_in(joint.parent.frame)
     0
     >>> joint.child.frame.dcm(joint.parent.frame)
@@ -924,7 +967,7 @@ class PrismaticJoint(Joint, _JointAxisMixin):
     [0, 1, 0],
     [0, 0, 1]])
     >>> joint.child_point.pos_from(joint.parent_point)
-    x_PC(t)*P_frame.x
+    q_PC(t)*P_frame.x
 
     To further demonstrate the use of the prismatic joint, the kinematics of two
     masses sliding, one moving relative to a fixed body and the other relative
@@ -969,10 +1012,10 @@ class PrismaticJoint(Joint, _JointAxisMixin):
     The position of the particles' masscenter is found with:
 
     >>> Part1.masscenter.pos_from(wall.masscenter)
-    x_J1(t)*W_frame.x
+    q_J1(t)*W_frame.x
 
     >>> Part2.masscenter.pos_from(wall.masscenter)
-    x_J1(t)*W_frame.x + x_J2(t)*P1_frame.x
+    q_J1(t)*W_frame.x + q_J2(t)*P1_frame.x
 
     The angular velocities of the two particle links can be computed with
     respect to the ceiling.
@@ -987,10 +1030,10 @@ class PrismaticJoint(Joint, _JointAxisMixin):
     with respect to the ceiling.
 
     >>> Part1.masscenter_vel(wall)
-    v_J1(t)*W_frame.x
+    u_J1(t)*W_frame.x
 
     >>> Part2.masscenter.vel(wall.frame)
-    v_J1(t)*W_frame.x + Derivative(x_J2(t), t)*P1_frame.x
+    u_J1(t)*W_frame.x + Derivative(q_J2(t), t)*P1_frame.x
 
     """
 
@@ -1015,20 +1058,10 @@ def joint_axis(self):
         return self._joint_axis
 
     def _generate_coordinates(self, coordinate):
-        coordinates = []
-        if coordinate is None:
-            x = dynamicsymbols('x' + '_' + self._name)
-            coordinate = x
-        coordinates.append(coordinate)
-        return coordinates
+        return self._fill_coordinate_list(coordinate, 1, 'q')
 
     def _generate_speeds(self, speed):
-        speeds = []
-        if speed is None:
-            y = dynamicsymbols('v' + '_' + self._name)
-            speed = y
-        speeds.append(speed)
-        return speeds
+        return self._fill_coordinate_list(speed, 1, 'u')
 
     def _orient_frames(self):
         self._joint_axis = self._axis(
diff --git a/sympy/physics/mechanics/jointsmethod.py b/sympy/physics/mechanics/jointsmethod.py
index f52470d9e6..7691a3c4ea 100644
--- a/sympy/physics/mechanics/jointsmethod.py
+++ b/sympy/physics/mechanics/jointsmethod.py
@@ -1,6 +1,7 @@
 from sympy.physics.mechanics import (Body, Lagrangian, KanesMethod, LagrangesMethod,
                                     RigidBody, Particle)
 from sympy.physics.mechanics.method import _Methods
+from sympy.core.backend import Matrix
 
 __all__ = ['JointsMethod']
 
@@ -161,7 +162,7 @@ def _generate_q(self):
                 if coordinate in q_ind:
                     raise ValueError('Coordinates of joints should be unique.')
                 q_ind.append(coordinate)
-        return q_ind
+        return Matrix(q_ind)
 
     def _generate_u(self):
         u_ind = []
@@ -170,12 +171,12 @@ def _generate_u(self):
                 if speed in u_ind:
                     raise ValueError('Speeds of joints should be unique.')
                 u_ind.append(speed)
-        return u_ind
+        return Matrix(u_ind)
 
     def _generate_kdes(self):
-        kd_ind = []
+        kd_ind = Matrix(1, 0, []).T
         for joint in self._joints:
-            kd_ind.extend(joint.kdes)
+            kd_ind = kd_ind.col_join(joint.kdes)
         return kd_ind
 
     def _convert_bodies(self):


Test patch:
diff --git a/sympy/functions/elementary/tests/test_integers.py b/sympy/functions/elementary/tests/test_integers.py
index 428b3ce8cf..73d97a593f 100644
--- a/sympy/functions/elementary/tests/test_integers.py
+++ b/sympy/functions/elementary/tests/test_integers.py
@@ -582,20 +582,6 @@ def test_issue_14355():
     assert ceiling((x**3 + x)/(x**2 - x)).series(x, 0, 100, cdir = -1) == 0
 
 
-def test_frac_leading_term():
-    assert frac(x).as_leading_term(x) == x
-    assert frac(x).as_leading_term(x, cdir = 1) == x
-    assert frac(x).as_leading_term(x, cdir = -1) == 1
-    assert frac(x + S.Half).as_leading_term(x, cdir = 1) == S.Half
-    assert frac(x + S.Half).as_leading_term(x, cdir = -1) == S.Half
-    assert frac(-2*x + 1).as_leading_term(x, cdir = 1) == S.One
-    assert frac(-2*x + 1).as_leading_term(x, cdir = -1) == -2*x
-    assert frac(sin(x) + 5).as_leading_term(x, cdir = 1) == x
-    assert frac(sin(x) + 5).as_leading_term(x, cdir = -1) == S.One
-    assert frac(sin(x**2) + 5).as_leading_term(x, cdir = 1) == x**2
-    assert frac(sin(x**2) + 5).as_leading_term(x, cdir = -1) == x**2
-
-
 @XFAIL
 def test_issue_4149():
     assert floor(3 + pi*I + y*I) == 3 + floor(pi + y)*I
diff --git a/sympy/physics/mechanics/tests/test_functions.py b/sympy/physics/mechanics/tests/test_functions.py
index 64adbdec08..05260375f0 100644
--- a/sympy/physics/mechanics/tests/test_functions.py
+++ b/sympy/physics/mechanics/tests/test_functions.py
@@ -8,10 +8,9 @@
                                      find_dynamicsymbols, Lagrangian)
 
 from sympy.physics.mechanics.functions import gravity, center_of_mass
-from sympy.physics.vector.vector import Vector
 from sympy.testing.pytest import raises
 
-Vector.simp = True
+
 q1, q2, q3, q4, q5 = symbols('q1 q2 q3 q4 q5')
 N = ReferenceFrame('N')
 A = N.orientnew('A', 'Axis', [q1, N.z])
diff --git a/sympy/physics/mechanics/tests/test_joint.py b/sympy/physics/mechanics/tests/test_joint.py
index 6d6ad1f4ef..6e372b2df7 100644
--- a/sympy/physics/mechanics/tests/test_joint.py
+++ b/sympy/physics/mechanics/tests/test_joint.py
@@ -3,14 +3,15 @@
 from sympy.core.singleton import S
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import (cos, sin)
-from sympy.matrices.dense import Matrix
-from sympy.core.backend import _simplify_matrix
+from sympy.core.backend import Matrix, _simplify_matrix
 from sympy.core.symbol import symbols
 from sympy.physics.mechanics import dynamicsymbols, Body, PinJoint, PrismaticJoint
 from sympy.physics.mechanics.joint import Joint
 from sympy.physics.vector import Vector, ReferenceFrame, Point
 from sympy.testing.pytest import raises, XFAIL, warns_deprecated_sympy
 
+
+Vector.simp = True
 t = dynamicsymbols._t # type: ignore
 
 
@@ -33,18 +34,58 @@ def test_Joint():
     raises(TypeError, lambda: Joint('J', parent, child))
 
 
+def test_coordinate_generation():
+    q, u, qj, uj = dynamicsymbols('q u q_J u_J')
+    q0j, q1j, q2j, q3j, u0j, u1j, u2j, u3j = dynamicsymbols('q0:4_J u0:4_J')
+    q0, q1, q2, q3, u0, u1, u2, u3 = dynamicsymbols('q0:4 u0:4')
+    _, _, P, C = _generate_body()
+    # Using PinJoint to access Joint's coordinate generation method
+    J = PinJoint('J', P, C)
+    # Test single given
+    assert J._fill_coordinate_list(q, 1) == Matrix([q])
+    assert J._fill_coordinate_list([u], 1) == Matrix([u])
+    assert J._fill_coordinate_list([u], 1, offset=2) == Matrix([u])
+    # Test None
+    assert J._fill_coordinate_list(None, 1) == Matrix([qj])
+    assert J._fill_coordinate_list([None], 1) == Matrix([qj])
+    assert J._fill_coordinate_list([q0, None], 3) == Matrix([q0, q1j, q2j])
+    # Test autofill
+    assert J._fill_coordinate_list(None, 3) == Matrix([q0j, q1j, q2j])
+    assert J._fill_coordinate_list([], 3) == Matrix([q0j, q1j, q2j])
+    # Test offset
+    assert J._fill_coordinate_list([], 3, offset=1) == Matrix([q1j, q2j, q3j])
+    assert J._fill_coordinate_list(q1, 3, offset=1) == Matrix([q1, q2j, q3j])
+    assert J._fill_coordinate_list([q1, None, q3], 3, offset=1) == Matrix(
+        [q1, q2j, q3])
+    assert J._fill_coordinate_list(None, 2, offset=2) == Matrix([q2j, q3j])
+    # Test label
+    assert J._fill_coordinate_list(None, 1, 'u') == Matrix([uj])
+    assert J._fill_coordinate_list([], 3, 'u') == Matrix([u0j, u1j, u2j])
+    assert J._fill_coordinate_list([u0], 3, 'u', 1) == Matrix([u0, u2j, u3j])
+    # Test single numbering
+    assert J._fill_coordinate_list(None, 1, number_single=True) == Matrix([q0j])
+    assert J._fill_coordinate_list([], 1, 'u', 2, True) == Matrix([u2j])
+    assert J._fill_coordinate_list([], 3, 'q') == Matrix([q0j, q1j, q2j])
+    # Test too many coordinates supplied
+    raises(ValueError, lambda: J._fill_coordinate_list([q0, q1], 1))
+    raises(ValueError, lambda: J._fill_coordinate_list([u0, u1, None], 2, 'u'))
+    # Test incorrect coordinate type
+    raises(TypeError, lambda: J._fill_coordinate_list([q0, symbols('q1')], 2))
+    raises(TypeError, lambda: J._fill_coordinate_list([q0 + q1, q1], 2))
+
+
 def test_pin_joint():
     P = Body('P')
     C = Body('C')
     l, m = symbols('l m')
-    theta, omega = dynamicsymbols('theta_J, omega_J')
+    q, u = dynamicsymbols('q_J, u_J')
     Pj = PinJoint('J', P, C)
     assert Pj.name == 'J'
     assert Pj.parent == P
     assert Pj.child == C
-    assert Pj.coordinates == [theta]
-    assert Pj.speeds == [omega]
-    assert Pj.kdes == [omega - theta.diff(t)]
+    assert Pj.coordinates == Matrix([q])
+    assert Pj.speeds == Matrix([u])
+    assert Pj.kdes == Matrix([u - q.diff(t)])
     assert Pj.joint_axis == P.frame.x
     assert Pj.child_point.pos_from(C.masscenter) == Vector(0)
     assert Pj.parent_point.pos_from(P.masscenter) == Vector(0)
@@ -121,8 +162,8 @@ def test_pin_joint_double_pendulum():
     assert B.ang_vel_in(N) == u1 * N.z + u2 * A.z
 
     # Check kde
-    assert J1.kdes == [u1 - q1.diff(t)]
-    assert J2.kdes == [u2 - q2.diff(t)]
+    assert J1.kdes == Matrix([u1 - q1.diff(t)])
+    assert J2.kdes == Matrix([u2 - q2.diff(t)])
 
     # Check Linear Velocity
     assert PartP.masscenter.vel(N) == l*u1*A.y
@@ -164,8 +205,8 @@ def test_pin_joint_chaos_pendulum():
     assert N.ang_vel_in(B) == -omega*N.y - alpha*A.z
 
     # Check kde
-    assert J1.kdes == [omega - theta.diff(t)]
-    assert J2.kdes == [alpha - phi.diff(t)]
+    assert J1.kdes == Matrix([omega - theta.diff(t)])
+    assert J2.kdes == Matrix([alpha - phi.diff(t)])
 
     # Check pos of masscenters
     assert C.masscenter.pos_from(rod.masscenter) == lA*A.z
@@ -326,19 +367,19 @@ def test_pin_joint_joint_axis():
 
 
 def test_pin_joint_arbitrary_axis():
-    theta, omega = dynamicsymbols('theta_J, omega_J')
+    q, u = dynamicsymbols('q_J, u_J')
 
-    # When the bodies are attached though masscenters but axess are opposite.
+    # When the bodies are attached though masscenters but axes are opposite.
     N, A, P, C = _generate_body()
     PinJoint('J', P, C, child_interframe=-A.x)
 
     assert (-A.x).angle_between(N.x) == 0
     assert -A.x.express(N) == N.x
     assert A.dcm(N) == Matrix([[-1, 0, 0],
-                            [0, -cos(theta), -sin(theta)],
-                            [0, -sin(theta), cos(theta)]])
-    assert A.ang_vel_in(N) == omega*N.x
-    assert A.ang_vel_in(N).magnitude() == sqrt(omega**2)
+                               [0, -cos(q), -sin(q)],
+                               [0, -sin(q), cos(q)]])
+    assert A.ang_vel_in(N) == u*N.x
+    assert A.ang_vel_in(N).magnitude() == sqrt(u**2)
     assert C.masscenter.pos_from(P.masscenter) == 0
     assert C.masscenter.pos_from(P.masscenter).express(N).simplify() == 0
     assert C.masscenter.vel(N) == 0
@@ -350,17 +391,17 @@ def test_pin_joint_arbitrary_axis():
 
     assert A.y.angle_between(N.x) == 0  # Axis are aligned
     assert A.y.express(N) == N.x
-    assert A.dcm(N) == Matrix([[0, -cos(theta), -sin(theta)],
+    assert A.dcm(N) == Matrix([[0, -cos(q), -sin(q)],
                                [1, 0, 0],
-                               [0, -sin(theta), cos(theta)]])
-    assert A.ang_vel_in(N) == omega*N.x
-    assert A.ang_vel_in(N).express(A) == omega * A.y
-    assert A.ang_vel_in(N).magnitude() == sqrt(omega**2)
+                               [0, -sin(q), cos(q)]])
+    assert A.ang_vel_in(N) == u*N.x
+    assert A.ang_vel_in(N).express(A) == u * A.y
+    assert A.ang_vel_in(N).magnitude() == sqrt(u**2)
     assert A.ang_vel_in(N).cross(A.y) == 0
-    assert C.masscenter.vel(N) == omega*A.z
+    assert C.masscenter.vel(N) == u*A.z
     assert C.masscenter.pos_from(P.masscenter) == -A.x
     assert (C.masscenter.pos_from(P.masscenter).express(N).simplify() ==
-            cos(theta)*N.y + sin(theta)*N.z)
+            cos(q)*N.y + sin(q)*N.z)
     assert C.masscenter.vel(N).angle_between(A.x) == pi/2
 
     # Similar to previous case but wrt parent body
@@ -370,13 +411,13 @@ def test_pin_joint_arbitrary_axis():
     assert N.y.angle_between(A.x) == 0  # Axis are aligned
     assert N.y.express(A) == A.x
     assert A.dcm(N) == Matrix([[0, 1, 0],
-                               [-cos(theta), 0, sin(theta)],
-                               [sin(theta), 0, cos(theta)]])
-    assert A.ang_vel_in(N) == omega*N.y
-    assert A.ang_vel_in(N).express(A) == omega*A.x
-    assert A.ang_vel_in(N).magnitude() == sqrt(omega**2)
+                               [-cos(q), 0, sin(q)],
+                               [sin(q), 0, cos(q)]])
+    assert A.ang_vel_in(N) == u*N.y
+    assert A.ang_vel_in(N).express(A) == u*A.x
+    assert A.ang_vel_in(N).magnitude() == sqrt(u**2)
     angle = A.ang_vel_in(N).angle_between(A.x)
-    assert angle.xreplace({omega: 1}) == 0
+    assert angle.xreplace({u: 1}) == 0
     assert C.masscenter.vel(N) == 0
     assert C.masscenter.pos_from(P.masscenter) == N.x
 
@@ -387,22 +428,22 @@ def test_pin_joint_arbitrary_axis():
     assert expand_mul(N.x.angle_between(A.x + A.y)) == 0  # Axis are aligned
     assert (A.x + A.y).express(N).simplify() == sqrt(2)*N.x
     assert _simplify_matrix(A.dcm(N)) == Matrix([
-        [sqrt(2)/2, -sqrt(2)*cos(theta)/2, -sqrt(2)*sin(theta)/2],
-        [sqrt(2)/2, sqrt(2)*cos(theta)/2, sqrt(2)*sin(theta)/2],
-        [0, -sin(theta), cos(theta)]])
-    assert A.ang_vel_in(N) == omega*N.x
+        [sqrt(2)/2, -sqrt(2)*cos(q)/2, -sqrt(2)*sin(q)/2],
+        [sqrt(2)/2, sqrt(2)*cos(q)/2, sqrt(2)*sin(q)/2],
+        [0, -sin(q), cos(q)]])
+    assert A.ang_vel_in(N) == u*N.x
     assert (A.ang_vel_in(N).express(A).simplify() ==
-            (omega*A.x + omega*A.y)/sqrt(2))
-    assert A.ang_vel_in(N).magnitude() == sqrt(omega**2)
+            (u*A.x + u*A.y)/sqrt(2))
+    assert A.ang_vel_in(N).magnitude() == sqrt(u**2)
     angle = A.ang_vel_in(N).angle_between(A.x + A.y)
-    assert angle.xreplace({omega: 1}) == 0
-    assert C.masscenter.vel(N).simplify() == (omega * A.z)/sqrt(2)
+    assert angle.xreplace({u: 1}) == 0
+    assert C.masscenter.vel(N).simplify() == (u * A.z)/sqrt(2)
     assert C.masscenter.pos_from(P.masscenter) == N.x - A.x
     assert (C.masscenter.pos_from(P.masscenter).express(N).simplify() ==
-            (1 - sqrt(2)/2)*N.x + sqrt(2)*cos(theta)/2*N.y +
-            sqrt(2)*sin(theta)/2*N.z)
+            (1 - sqrt(2)/2)*N.x + sqrt(2)*cos(q)/2*N.y +
+            sqrt(2)*sin(q)/2*N.z)
     assert (C.masscenter.vel(N).express(N).simplify() ==
-            -sqrt(2)*omega*sin(theta)/2*N.y + sqrt(2)*omega*cos(theta)/2*N.z)
+            -sqrt(2)*u*sin(q)/2*N.y + sqrt(2)*u*cos(q)/2*N.z)
     assert C.masscenter.vel(N).angle_between(A.x) == pi/2
 
     N, A, P, C = _generate_body()
@@ -411,26 +452,26 @@ def test_pin_joint_arbitrary_axis():
     assert expand_mul(N.x.angle_between(A.x + A.y - A.z)) == 0  # Axis aligned
     assert (A.x + A.y - A.z).express(N).simplify() == sqrt(3)*N.x
     assert _simplify_matrix(A.dcm(N)) == Matrix([
-        [sqrt(3)/3, -sqrt(6)*sin(theta + pi/4)/3,
-         sqrt(6)*cos(theta + pi/4)/3],
-        [sqrt(3)/3, sqrt(6)*cos(theta + pi/12)/3,
-         sqrt(6)*sin(theta + pi/12)/3],
-        [-sqrt(3)/3, sqrt(6)*cos(theta + 5*pi/12)/3,
-         sqrt(6)*sin(theta + 5*pi/12)/3]])
-    assert A.ang_vel_in(N) == omega*N.x
-    assert A.ang_vel_in(N).express(A).simplify() == (omega*A.x + omega*A.y -
-                                                     omega*A.z)/sqrt(3)
-    assert A.ang_vel_in(N).magnitude() == sqrt(omega**2)
+        [sqrt(3)/3, -sqrt(6)*sin(q + pi/4)/3,
+         sqrt(6)*cos(q + pi/4)/3],
+        [sqrt(3)/3, sqrt(6)*cos(q + pi/12)/3,
+         sqrt(6)*sin(q + pi/12)/3],
+        [-sqrt(3)/3, sqrt(6)*cos(q + 5*pi/12)/3,
+         sqrt(6)*sin(q + 5*pi/12)/3]])
+    assert A.ang_vel_in(N) == u*N.x
+    assert A.ang_vel_in(N).express(A).simplify() == (u*A.x + u*A.y -
+                                                     u*A.z)/sqrt(3)
+    assert A.ang_vel_in(N).magnitude() == sqrt(u**2)
     angle = A.ang_vel_in(N).angle_between(A.x + A.y-A.z)
-    assert angle.xreplace({omega: 1}) == 0
-    assert C.masscenter.vel(N).simplify() == (omega*A.y + omega*A.z)/sqrt(3)
+    assert angle.xreplace({u: 1}) == 0
+    assert C.masscenter.vel(N).simplify() == (u*A.y + u*A.z)/sqrt(3)
     assert C.masscenter.pos_from(P.masscenter) == N.x - A.x
     assert (C.masscenter.pos_from(P.masscenter).express(N).simplify() ==
-            (1 - sqrt(3)/3)*N.x + sqrt(6)*sin(theta + pi/4)/3*N.y -
-            sqrt(6)*cos(theta + pi/4)/3*N.z)
+            (1 - sqrt(3)/3)*N.x + sqrt(6)*sin(q + pi/4)/3*N.y -
+            sqrt(6)*cos(q + pi/4)/3*N.z)
     assert (C.masscenter.vel(N).express(N).simplify() ==
-            sqrt(6)*omega*cos(theta + pi/4)/3*N.y +
-            sqrt(6)*omega*sin(theta + pi/4)/3*N.z)
+            sqrt(6)*u*cos(q + pi/4)/3*N.y +
+            sqrt(6)*u*sin(q + pi/4)/3*N.z)
     assert C.masscenter.vel(N).angle_between(A.x) == pi/2
 
     N, A, P, C = _generate_body()
@@ -441,30 +482,30 @@ def test_pin_joint_arbitrary_axis():
     angle = (N.x - N.y + N.z).angle_between(A.x + A.y - A.z)
     assert expand_mul(angle) == 0  # Axis are aligned
     assert ((A.x-A.y+A.z).express(N).simplify() ==
-            (-4*cos(theta)/3 - S(1)/3)*N.x + (S(1)/3 - 4*sin(theta + pi/6)/3)*N.y +
-            (4*cos(theta + pi/3)/3 - S(1)/3)*N.z)
+            (-4*cos(q)/3 - S(1)/3)*N.x + (S(1)/3 - 4*sin(q + pi/6)/3)*N.y +
+            (4*cos(q + pi/3)/3 - S(1)/3)*N.z)
     assert _simplify_matrix(A.dcm(N)) == Matrix([
-        [S(1)/3 - 2*cos(theta)/3, -2*sin(theta + pi/6)/3 - S(1)/3,
-         2*cos(theta + pi/3)/3 + S(1)/3],
-        [2*cos(theta + pi/3)/3 + S(1)/3, 2*cos(theta)/3 - S(1)/3,
-         2*sin(theta + pi/6)/3 + S(1)/3],
-        [-2*sin(theta + pi/6)/3 - S(1)/3, 2*cos(theta + pi/3)/3 + S(1)/3,
-         2*cos(theta)/3 - S(1)/3]])
-    assert A.ang_vel_in(N) == (omega*N.x - omega*N.y + omega*N.z)/sqrt(3)
-    assert A.ang_vel_in(N).express(A).simplify() == (omega*A.x + omega*A.y -
-                                                     omega*A.z)/sqrt(3)
-    assert A.ang_vel_in(N).magnitude() == sqrt(omega**2)
+        [S(1)/3 - 2*cos(q)/3, -2*sin(q + pi/6)/3 - S(1)/3,
+         2*cos(q + pi/3)/3 + S(1)/3],
+        [2*cos(q + pi/3)/3 + S(1)/3, 2*cos(q)/3 - S(1)/3,
+         2*sin(q + pi/6)/3 + S(1)/3],
+        [-2*sin(q + pi/6)/3 - S(1)/3, 2*cos(q + pi/3)/3 + S(1)/3,
+         2*cos(q)/3 - S(1)/3]])
+    assert A.ang_vel_in(N) == (u*N.x - u*N.y + u*N.z)/sqrt(3)
+    assert A.ang_vel_in(N).express(A).simplify() == (u*A.x + u*A.y -
+                                                     u*A.z)/sqrt(3)
+    assert A.ang_vel_in(N).magnitude() == sqrt(u**2)
     angle = A.ang_vel_in(N).angle_between(A.x+A.y-A.z)
-    assert angle.xreplace({omega: 1}) == 0
+    assert angle.xreplace({u: 1}) == 0
     assert (C.masscenter.vel(N).simplify() ==
-            sqrt(3)*n*omega/3*A.y + sqrt(3)*n*omega/3*A.z)
+            sqrt(3)*n*u/3*A.y + sqrt(3)*n*u/3*A.z)
     assert C.masscenter.pos_from(P.masscenter) == m*N.x - n*A.x
     assert (C.masscenter.pos_from(P.masscenter).express(N).simplify() ==
-            (m + n*(2*cos(theta) - 1)/3)*N.x + n*(2*sin(theta + pi/6) +
-            1)/3*N.y - n*(2*cos(theta + pi/3) + 1)/3*N.z)
+            (m + n*(2*cos(q) - 1)/3)*N.x + n*(2*sin(q + pi/6) +
+            1)/3*N.y - n*(2*cos(q + pi/3) + 1)/3*N.z)
     assert (C.masscenter.vel(N).express(N).simplify() ==
-            - 2*n*omega*sin(theta)/3*N.x + 2*n*omega*cos(theta + pi/6)/3*N.y +
-            2*n*omega*sin(theta + pi/3)/3*N.z)
+            - 2*n*u*sin(q)/3*N.x + 2*n*u*cos(q + pi/6)/3*N.y +
+            2*n*u*sin(q + pi/3)/3*N.z)
     assert C.masscenter.vel(N).dot(N.x - N.y + N.z).simplify() == 0
 
 
@@ -573,20 +614,20 @@ def test_locate_joint_frame():
 
 def test_sliding_joint():
     _, _, P, C = _generate_body()
-    x, v = dynamicsymbols('x_S, v_S')
+    q, u = dynamicsymbols('q_S, u_S')
     S = PrismaticJoint('S', P, C)
     assert S.name == 'S'
     assert S.parent == P
     assert S.child == C
-    assert S.coordinates == [x]
-    assert S.speeds == [v]
-    assert S.kdes == [v - x.diff(t)]
+    assert S.coordinates == Matrix([q])
+    assert S.speeds == Matrix([u])
+    assert S.kdes == Matrix([u - q.diff(t)])
     assert S.joint_axis == P.frame.x
     assert S.child_point.pos_from(C.masscenter) == Vector(0)
     assert S.parent_point.pos_from(P.masscenter) == Vector(0)
-    assert S.parent_point.pos_from(S.child_point) == - x * P.frame.x
-    assert P.masscenter.pos_from(C.masscenter) == - x * P.frame.x
-    assert C.masscenter.vel(P.frame) == v * P.frame.x
+    assert S.parent_point.pos_from(S.child_point) == - q * P.frame.x
+    assert P.masscenter.pos_from(C.masscenter) == - q * P.frame.x
+    assert C.masscenter.vel(P.frame) == u * P.frame.x
     assert P.ang_vel_in(C) == 0
     assert C.ang_vel_in(P) == 0
     assert S.__str__() == 'PrismaticJoint: S  parent: P  child: C'
@@ -602,9 +643,9 @@ def test_sliding_joint():
     assert S.joint_axis == P.frame.z
     assert S.child_point.pos_from(C.masscenter) == m * C.frame.y
     assert S.parent_point.pos_from(P.masscenter) == l * P.frame.x
-    assert S.parent_point.pos_from(S.child_point) == - x * P.frame.z
-    assert P.masscenter.pos_from(C.masscenter) == - l*N.x - x*N.z + m*A.y
-    assert C.masscenter.vel(P.frame) == v * P.frame.z
+    assert S.parent_point.pos_from(S.child_point) == - q * P.frame.z
+    assert P.masscenter.pos_from(C.masscenter) == - l*N.x - q*N.z + m*A.y
+    assert C.masscenter.vel(P.frame) == u * P.frame.z
     assert P.masscenter.vel(Pint) == Vector(0)
     assert C.ang_vel_in(P) == 0
     assert P.ang_vel_in(C) == 0
@@ -618,15 +659,15 @@ def test_sliding_joint():
     assert S.joint_axis == P.frame.z
     assert S.child_point.pos_from(C.masscenter) == m * C.frame.x
     assert S.parent_point.pos_from(P.masscenter) == l * P.frame.z
-    assert S.parent_point.pos_from(S.child_point) == - x * P.frame.z
-    assert P.masscenter.pos_from(C.masscenter) == (-l - x)*P.frame.z + m*C.frame.x
-    assert C.masscenter.vel(P.frame) == v * P.frame.z
+    assert S.parent_point.pos_from(S.child_point) == - q * P.frame.z
+    assert P.masscenter.pos_from(C.masscenter) == (-l - q)*P.frame.z + m*C.frame.x
+    assert C.masscenter.vel(P.frame) == u * P.frame.z
     assert C.ang_vel_in(P) == 0
     assert P.ang_vel_in(C) == 0
 
 
 def test_sliding_joint_arbitrary_axis():
-    x, v = dynamicsymbols('x_S, v_S')
+    q, u = dynamicsymbols('q_S, u_S')
 
     N, A, P, C = _generate_body()
     PrismaticJoint('S', P, C, child_interframe=-A.x)
@@ -634,10 +675,10 @@ def test_sliding_joint_arbitrary_axis():
     assert (-A.x).angle_between(N.x) == 0
     assert -A.x.express(N) == N.x
     assert A.dcm(N) == Matrix([[-1, 0, 0], [0, -1, 0], [0, 0, 1]])
-    assert C.masscenter.pos_from(P.masscenter) == x * N.x
-    assert C.masscenter.pos_from(P.masscenter).express(A).simplify() == -x * A.x
-    assert C.masscenter.vel(N) == v * N.x
-    assert C.masscenter.vel(N).express(A) == -v * A.x
+    assert C.masscenter.pos_from(P.masscenter) == q * N.x
+    assert C.masscenter.pos_from(P.masscenter).express(A).simplify() == -q * A.x
+    assert C.masscenter.vel(N) == u * N.x
+    assert C.masscenter.vel(N).express(A) == -u * A.x
     assert A.ang_vel_in(N) == 0
     assert N.ang_vel_in(A) == 0
 
@@ -649,10 +690,10 @@ def test_sliding_joint_arbitrary_axis():
     assert A.y.angle_between(N.x) == 0 #Axis are aligned
     assert A.y.express(N) == N.x
     assert A.dcm(N) == Matrix([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
-    assert C.masscenter.vel(N) == v * N.x
-    assert C.masscenter.vel(N).express(A) == v * A.y
-    assert C.masscenter.pos_from(P.masscenter) == x*N.x - A.x
-    assert C.masscenter.pos_from(P.masscenter).express(N).simplify() == x*N.x + N.y
+    assert C.masscenter.vel(N) == u * N.x
+    assert C.masscenter.vel(N).express(A) == u * A.y
+    assert C.masscenter.pos_from(P.masscenter) == q*N.x - A.x
+    assert C.masscenter.pos_from(P.masscenter).express(N).simplify() == q*N.x + N.y
     assert A.ang_vel_in(N) == 0
     assert N.ang_vel_in(A) == 0
 
@@ -663,9 +704,9 @@ def test_sliding_joint_arbitrary_axis():
     assert N.y.angle_between(A.x) == 0 #Axis are aligned
     assert N.y.express(A) ==  A.x
     assert A.dcm(N) == Matrix([[0, 1, 0], [-1, 0, 0], [0, 0, 1]])
-    assert C.masscenter.vel(N) == v * N.y
-    assert C.masscenter.vel(N).express(A) == v * A.x
-    assert C.masscenter.pos_from(P.masscenter) == N.x + x*N.y
+    assert C.masscenter.vel(N) == u * N.y
+    assert C.masscenter.vel(N).express(A) == u * A.x
+    assert C.masscenter.pos_from(P.masscenter) == N.x + q*N.y
     assert A.ang_vel_in(N) == 0
     assert N.ang_vel_in(A) == 0
 
@@ -676,10 +717,10 @@ def test_sliding_joint_arbitrary_axis():
     assert N.x.angle_between(A.x + A.y) == 0 #Axis are aligned
     assert (A.x + A.y).express(N) == sqrt(2)*N.x
     assert A.dcm(N) == Matrix([[sqrt(2)/2, -sqrt(2)/2, 0], [sqrt(2)/2, sqrt(2)/2, 0], [0, 0, 1]])
-    assert C.masscenter.pos_from(P.masscenter) == (x + 1)*N.x - A.x
-    assert C.masscenter.pos_from(P.masscenter).express(N) == (x - sqrt(2)/2 + 1)*N.x + sqrt(2)/2*N.y
-    assert C.masscenter.vel(N).express(A) == v * (A.x + A.y)/sqrt(2)
-    assert C.masscenter.vel(N) == v*N.x
+    assert C.masscenter.pos_from(P.masscenter) == (q + 1)*N.x - A.x
+    assert C.masscenter.pos_from(P.masscenter).express(N) == (q - sqrt(2)/2 + 1)*N.x + sqrt(2)/2*N.y
+    assert C.masscenter.vel(N).express(A) == u * (A.x + A.y)/sqrt(2)
+    assert C.masscenter.vel(N) == u*N.x
     assert A.ang_vel_in(N) == 0
     assert N.ang_vel_in(A) == 0
 
@@ -691,11 +732,11 @@ def test_sliding_joint_arbitrary_axis():
     assert _simplify_matrix(A.dcm(N)) == Matrix([[sqrt(3)/3, -sqrt(3)/3, sqrt(3)/3],
                                                  [sqrt(3)/3, sqrt(3)/6 + S(1)/2, S(1)/2 - sqrt(3)/6],
                                                  [-sqrt(3)/3, S(1)/2 - sqrt(3)/6, sqrt(3)/6 + S(1)/2]])
-    assert C.masscenter.pos_from(P.masscenter) == (x + 1)*N.x - A.x
+    assert C.masscenter.pos_from(P.masscenter) == (q + 1)*N.x - A.x
     assert C.masscenter.pos_from(P.masscenter).express(N) == \
-        (x - sqrt(3)/3 + 1)*N.x + sqrt(3)/3*N.y - sqrt(3)/3*N.z
-    assert C.masscenter.vel(N) == v*N.x
-    assert C.masscenter.vel(N).express(A) == sqrt(3)*v/3*A.x + sqrt(3)*v/3*A.y - sqrt(3)*v/3*A.z
+        (q - sqrt(3)/3 + 1)*N.x + sqrt(3)/3*N.y - sqrt(3)/3*N.z
+    assert C.masscenter.vel(N) == u*N.x
+    assert C.masscenter.vel(N).express(A) == sqrt(3)*u/3*A.x + sqrt(3)*u/3*A.y - sqrt(3)*u/3*A.z
     assert A.ang_vel_in(N) == 0
     assert N.ang_vel_in(A) == 0
 
@@ -711,27 +752,27 @@ def test_sliding_joint_arbitrary_axis():
                                                  [S(2)/3, S(1)/3, S(2)/3],
                                                  [-S(2)/3, S(2)/3, S(1)/3]])
     assert C.masscenter.pos_from(P.masscenter) == \
-        (m + sqrt(3)*x/3)*N.x - sqrt(3)*x/3*N.y + sqrt(3)*x/3*N.z - n*A.x
+        (m + sqrt(3)*q/3)*N.x - sqrt(3)*q/3*N.y + sqrt(3)*q/3*N.z - n*A.x
     assert C.masscenter.pos_from(P.masscenter).express(N) == \
-        (m + n/3 + sqrt(3)*x/3)*N.x + (2*n/3 - sqrt(3)*x/3)*N.y + (-2*n/3 + sqrt(3)*x/3)*N.z
-    assert C.masscenter.vel(N) == sqrt(3)*v/3*N.x - sqrt(3)*v/3*N.y + sqrt(3)*v/3*N.z
-    assert C.masscenter.vel(N).express(A) == sqrt(3)*v/3*A.x + sqrt(3)*v/3*A.y - sqrt(3)*v/3*A.z
+        (m + n/3 + sqrt(3)*q/3)*N.x + (2*n/3 - sqrt(3)*q/3)*N.y + (-2*n/3 + sqrt(3)*q/3)*N.z
+    assert C.masscenter.vel(N) == sqrt(3)*u/3*N.x - sqrt(3)*u/3*N.y + sqrt(3)*u/3*N.z
+    assert C.masscenter.vel(N).express(A) == sqrt(3)*u/3*A.x + sqrt(3)*u/3*A.y - sqrt(3)*u/3*A.z
     assert A.ang_vel_in(N) == 0
     assert N.ang_vel_in(A) == 0
 
 
 def test_deprecated_parent_child_axis():
-    theta, omega = dynamicsymbols('theta_J, omega_J')
+    q, u = dynamicsymbols('q_J, u_J')
     N, A, P, C = _generate_body()
     with warns_deprecated_sympy():
         PinJoint('J', P, C, child_axis=-A.x)
     assert (-A.x).angle_between(N.x) == 0
     assert -A.x.express(N) == N.x
     assert A.dcm(N) == Matrix([[-1, 0, 0],
-                               [0, -cos(theta), -sin(theta)],
-                               [0, -sin(theta), cos(theta)]])
-    assert A.ang_vel_in(N) == omega * N.x
-    assert A.ang_vel_in(N).magnitude() == sqrt(omega ** 2)
+                               [0, -cos(q), -sin(q)],
+                               [0, -sin(q), cos(q)]])
+    assert A.ang_vel_in(N) == u * N.x
+    assert A.ang_vel_in(N).magnitude() == sqrt(u ** 2)
 
     N, A, P, C = _generate_body()
     with warns_deprecated_sympy():
diff --git a/sympy/physics/mechanics/tests/test_jointsmethod.py b/sympy/physics/mechanics/tests/test_jointsmethod.py
index fbcc77b029..f3d9a80f5a 100644
--- a/sympy/physics/mechanics/tests/test_jointsmethod.py
+++ b/sympy/physics/mechanics/tests/test_jointsmethod.py
@@ -20,18 +20,18 @@ def test_jointsmethod():
     C = Body('C')
     Pin = PinJoint('P1', P, C)
     C_ixx, g = symbols('C_ixx g')
-    theta, omega = dynamicsymbols('theta_P1, omega_P1')
+    q, u = dynamicsymbols('q_P1, u_P1')
     P.apply_force(g*P.y)
     method = JointsMethod(P, Pin)
     assert method.frame == P.frame
     assert method.bodies == [C, P]
     assert method.loads == [(P.masscenter, g*P.frame.y)]
-    assert method.q == [theta]
-    assert method.u == [omega]
-    assert method.kdes == [omega - theta.diff()]
+    assert method.q == Matrix([q])
+    assert method.u == Matrix([u])
+    assert method.kdes == Matrix([u - q.diff()])
     soln = method.form_eoms()
-    assert soln == Matrix([[-C_ixx*omega.diff()]])
-    assert method.forcing_full == Matrix([[omega], [0]])
+    assert soln == Matrix([[-C_ixx*u.diff()]])
+    assert method.forcing_full == Matrix([[u], [0]])
     assert method.mass_matrix_full == Matrix([[1, 0], [0, C_ixx]])
     assert isinstance(method.method, KanesMethod)
 
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index 17686ce2e1..baf78e7e61 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -244,22 +244,6 @@ def test_ceiling_requires_robust_assumptions():
     assert limit(ceiling(5 + cos(x)), x, 0, "-") == 6
 
 
-def test_frac():
-    assert limit(frac(x), x, oo) == AccumBounds(0, 1)
-    assert limit(frac(x)**(1/x), x, oo) == AccumBounds(0, 1)
-    assert limit(frac(x)**(1/x), x, -oo) == AccumBounds(1, oo)
-    assert limit(frac(x)**x, x, oo) == AccumBounds(0, oo)  # wolfram gives (0, 1)
-    assert limit(frac(sin(x)), x, 0, "+") == 0
-    assert limit(frac(sin(x)), x, 0, "-") == 1
-    assert limit(frac(cos(x)), x, 0, "+-") == 1
-    assert limit(frac(x**2), x, 0, "+-") == 0
-    raises(ValueError, lambda: limit(frac(x), x, 0, '+-'))
-    assert limit(frac(-2*x + 1), x, 0, "+") == 1
-    assert limit(frac(-2*x + 1), x, 0, "-") == 0
-    assert limit(frac(x + S.Half), x, 0, "+-") == 1/2
-    assert limit(frac(1/x), x, 0) == AccumBounds(0, 1)
-
-
 def test_issue_14355():
     assert limit(floor(sin(x)/x), x, 0, '+') == 0
     assert limit(floor(sin(x)/x), x, 0, '-') == 0
@@ -356,6 +340,7 @@ def test_series_AccumBounds():
     t2 = Mul(AccumBounds(-1 + sin(1)/2, sin(1)/2 + 1), 1/(1 - cos(1)))
     assert limit(simplify(Sum(sin(n).rewrite(exp), (n, 0, k)).doit().rewrite(sin)), k, oo) == t2
 
+    assert limit(frac(x)**x, x, oo) == AccumBounds(0, oo)  # wolfram gives (0, 1)
     assert limit(((sin(x) + 1)/2)**x, x, oo) == AccumBounds(0, oo)  # wolfram says 0
 
     # https://github.com/sympy/sympy/issues/12312
diff --git a/sympy/series/tests/test_nseries.py b/sympy/series/tests/test_nseries.py
index a2f20add82..87d4eece6e 100644
--- a/sympy/series/tests/test_nseries.py
+++ b/sympy/series/tests/test_nseries.py
@@ -1,4 +1,3 @@
-from sympy.calculus.util import AccumBounds
 from sympy.core.function import (Derivative, PoleError)
 from sympy.core.numbers import (E, I, Integer, Rational, pi)
 from sympy.core.singleton import S
@@ -6,7 +5,7 @@
 from sympy.functions.elementary.complexes import sign
 from sympy.functions.elementary.exponential import (exp, log)
 from sympy.functions.elementary.hyperbolic import (acosh, acoth, asinh, atanh, cosh, coth, sinh, tanh)
-from sympy.functions.elementary.integers import (ceiling, floor, frac)
+from sympy.functions.elementary.integers import (ceiling, floor)
 from sympy.functions.elementary.miscellaneous import (cbrt, sqrt)
 from sympy.functions.elementary.trigonometric import (asin, cos, cot, sin, tan)
 from sympy.series.limits import limit
@@ -445,20 +444,6 @@ def test_floor():
     assert floor(x + 1.5).series(x) == 1
 
 
-def test_frac():
-    assert frac(x).series(x, cdir=1) == x
-    assert frac(x).series(x, cdir=-1) == 1 + x
-    assert frac(2*x + 1).series(x, cdir=1) == 2*x
-    assert frac(2*x + 1).series(x, cdir=-1) == 1 + 2*x
-    assert frac(x**2).series(x, cdir=1) == x**2
-    assert frac(x**2).series(x, cdir=-1) == x**2
-    assert frac(sin(x) + 5).series(x, cdir=1) == x - x**3/6 + x**5/120 + O(x**6)
-    assert frac(sin(x) + 5).series(x, cdir=-1) == 1 + x - x**3/6 + x**5/120 + O(x**6)
-    assert frac(sin(x) + S.Half).series(x) == S.Half + x - x**3/6 + x**5/120 + O(x**6)
-    assert frac(x**8).series(x, cdir=1) == O(x**6)
-    assert frac(1/x).series(x) == AccumBounds(0, 1) + O(x**6)
-
-
 def test_ceiling():
     assert ceiling(x).series(x) == 1
     assert ceiling(-x).series(x) == 0