Repository: sympy/sympy

Issue:
test_kane3 fails in master
I've been doing some investigation on `sympy/physics/mechanics/tests/test_kane3.py`. This test implements the bicycle test for mechanics that is a comprehensive test that exercises all of the details of the KanesMethod and Linearizer classes.

The test was introduced in PR #2182 (commit 919853f) and it worked then as a slow test which takes about 3-4 minutes on my machine.

The test was then removed from the slow tests on Travis mistakenly. I believe the code was broken by a commit to master and it was misdiagnosed as a too slow test on Travis, but it actually fails. So I'm not completely sure where the break is, but a git bisect based on timings of the completion of the function found in commit 6800beb.

The current failure in master is such that it takes a long time to complete >1 hour on my machine then when the numerics are subbed the comparisons fail because of NaNs in the expressions.


Discussion:

Patch:
diff --git a/.ci/blacklisted.json b/.ci/blacklisted.json
index 3075f09d43..857e76864e 100644
--- a/.ci/blacklisted.json
+++ b/.ci/blacklisted.json
@@ -1,7 +1,4 @@
 {
-    "sympy/physics/mechanics/tests/test_kane3.py": [
-        "test_bicycle"
-    ],
     "sympy/utilities/tests/test_wester.py": [
         "test_W25"
     ]
diff --git a/.ci/generate_durations_log.sh b/.ci/generate_durations_log.sh
index 8d8c087b7d..26a91d59e0 100755
--- a/.ci/generate_durations_log.sh
+++ b/.ci/generate_durations_log.sh
@@ -2,9 +2,6 @@
 ABS_REPO_PATH=$(unset CDPATH && cd "$(dirname "$0")/.." && echo $PWD)
 cat <<EOF >${ABS_REPO_PATH}/.ci/blacklisted.json
 {
-    "sympy/physics/mechanics/tests/test_kane3.py": [
-        "test_bicycle"
-    ],
     "sympy/utilities/tests/test_wester.py": [
         "test_W25"
     ]
diff --git a/sympy/physics/mechanics/kane.py b/sympy/physics/mechanics/kane.py
index 0df1fde869..e1e126a6ec 100644
--- a/sympy/physics/mechanics/kane.py
+++ b/sympy/physics/mechanics/kane.py
@@ -5,9 +5,10 @@
 from sympy.physics.mechanics.method import _Methods
 from sympy.physics.mechanics.particle import Particle
 from sympy.physics.mechanics.rigidbody import RigidBody
-from sympy.physics.mechanics.functions import (
-    msubs, find_dynamicsymbols, _f_list_parser, _validate_coordinates,
-    _parse_linear_solver)
+from sympy.physics.mechanics.functions import (msubs, find_dynamicsymbols,
+                                               _f_list_parser,
+                                               _validate_coordinates,
+                                               _parse_linear_solver)
 from sympy.physics.mechanics.linearize import Linearizer
 from sympy.utilities.iterables import iterable
 
@@ -494,11 +495,32 @@ def get_partial_velocity(body):
         self._f_d = -(self._fr - nonMM)
         return fr_star
 
-    def to_linearizer(self):
+    def to_linearizer(self, linear_solver='LU'):
         """Returns an instance of the Linearizer class, initiated from the
         data in the KanesMethod class. This may be more desirable than using
         the linearize class method, as the Linearizer object will allow more
-        efficient recalculation (i.e. about varying operating points)."""
+        efficient recalculation (i.e. about varying operating points).
+
+        Parameters
+        ==========
+        linear_solver : str, callable
+            Method used to solve the several symbolic linear systems of the
+            form ``A*x=b`` in the linearization process. If a string is
+            supplied, it should be a valid method that can be used with the
+            :meth:`sympy.matrices.matrices.MatrixBase.solve`. If a callable is
+            supplied, it should have the format ``x = f(A, b)``, where it
+            solves the equations and returns the solution. The default is
+            ``'LU'`` which corresponds to SymPy's ``A.LUsolve(b)``.
+            ``LUsolve()`` is fast to compute but will often result in
+            divide-by-zero and thus ``nan`` results.
+
+        Returns
+        =======
+        Linearizer
+            An instantiated
+            :class:`sympy.physics.mechanics.linearize.Linearizer`.
+
+        """
 
         if (self._fr is None) or (self._frstar is None):
             raise ValueError('Need to compute Fr, Fr* first.')
@@ -565,12 +587,30 @@ def to_linearizer(self):
                 raise ValueError('Cannot have derivatives of specified \
                                  quantities when linearizing forcing terms.')
         return Linearizer(f_0, f_1, f_2, f_3, f_4, f_c, f_v, f_a, q, u, q_i,
-                q_d, u_i, u_d, r)
+                q_d, u_i, u_d, r, linear_solver=linear_solver)
 
     # TODO : Remove `new_method` after 1.1 has been released.
-    def linearize(self, *, new_method=None, **kwargs):
+    def linearize(self, *, new_method=None, linear_solver='LU', **kwargs):
         """ Linearize the equations of motion about a symbolic operating point.
 
+        Parameters
+        ==========
+        new_method
+            Deprecated, does nothing and will be removed.
+        linear_solver : str, callable
+            Method used to solve the several symbolic linear systems of the
+            form ``A*x=b`` in the linearization process. If a string is
+            supplied, it should be a valid method that can be used with the
+            :meth:`sympy.matrices.matrices.MatrixBase.solve`. If a callable is
+            supplied, it should have the format ``x = f(A, b)``, where it
+            solves the equations and returns the solution. The default is
+            ``'LU'`` which corresponds to SymPy's ``A.LUsolve(b)``.
+            ``LUsolve()`` is fast to compute but will often result in
+            divide-by-zero and thus ``nan`` results.
+        **kwargs
+            Extra keyword arguments are passed to
+            :meth:`sympy.physics.mechanics.linearize.Linearizer.linearize`.
+
         Explanation
         ===========
 
@@ -594,8 +634,11 @@ def linearize(self, *, new_method=None, **kwargs):
         dictionaries. The values may be numeric or symbolic. The more values
         you can specify beforehand, the faster this computation will run.
 
-        For more documentation, please see the ``Linearizer`` class."""
-        linearizer = self.to_linearizer()
+        For more documentation, please see the ``Linearizer`` class.
+
+        """
+
+        linearizer = self.to_linearizer(linear_solver=linear_solver)
         result = linearizer.linearize(**kwargs)
         return result + (linearizer.r,)
 
diff --git a/sympy/physics/mechanics/lagrange.py b/sympy/physics/mechanics/lagrange.py
index 10a111ee2e..531994d35c 100644
--- a/sympy/physics/mechanics/lagrange.py
+++ b/sympy/physics/mechanics/lagrange.py
@@ -289,10 +289,11 @@ def forcing_full(self):
         else:
             return self._qdots.col_join(self.forcing)
 
-    def to_linearizer(self, q_ind=None, qd_ind=None, q_dep=None, qd_dep=None):
-        """Returns an instance of the Linearizer class, initiated from the
-        data in the LagrangesMethod class. This may be more desirable than using
-        the linearize class method, as the Linearizer object will allow more
+    def to_linearizer(self, q_ind=None, qd_ind=None, q_dep=None, qd_dep=None,
+                      linear_solver='LU'):
+        """Returns an instance of the Linearizer class, initiated from the data
+        in the LagrangesMethod class. This may be more desirable than using the
+        linearize class method, as the Linearizer object will allow more
         efficient recalculation (i.e. about varying operating points).
 
         Parameters
@@ -302,6 +303,23 @@ def to_linearizer(self, q_ind=None, qd_ind=None, q_dep=None, qd_dep=None):
             The independent generalized coordinates and speeds.
         q_dep, qd_dep : array_like, optional
             The dependent generalized coordinates and speeds.
+        linear_solver : str, callable
+            Method used to solve the several symbolic linear systems of the
+            form ``A*x=b`` in the linearization process. If a string is
+            supplied, it should be a valid method that can be used with the
+            :meth:`sympy.matrices.matrices.MatrixBase.solve`. If a callable is
+            supplied, it should have the format ``x = f(A, b)``, where it
+            solves the equations and returns the solution. The default is
+            ``'LU'`` which corresponds to SymPy's ``A.LUsolve(b)``.
+            ``LUsolve()`` is fast to compute but will often result in
+            divide-by-zero and thus ``nan`` results.
+
+        Returns
+        =======
+        Linearizer
+            An instantiated
+            :class:`sympy.physics.mechanics.linearize.Linearizer`.
+
         """
 
         # Compose vectors
@@ -352,12 +370,28 @@ def to_linearizer(self, q_ind=None, qd_ind=None, q_dep=None, qd_dep=None):
                                  quantities when linearizing forcing terms.')
 
         return Linearizer(f_0, f_1, f_2, f_3, f_4, f_c, f_v, f_a, q, u, q_i,
-                q_d, u_i, u_d, r, lams)
+                          q_d, u_i, u_d, r, lams, linear_solver=linear_solver)
 
     def linearize(self, q_ind=None, qd_ind=None, q_dep=None, qd_dep=None,
-            **kwargs):
+                  linear_solver='LU', **kwargs):
         """Linearize the equations of motion about a symbolic operating point.
 
+        Parameters
+        ==========
+        linear_solver : str, callable
+            Method used to solve the several symbolic linear systems of the
+            form ``A*x=b`` in the linearization process. If a string is
+            supplied, it should be a valid method that can be used with the
+            :meth:`sympy.matrices.matrices.MatrixBase.solve`. If a callable is
+            supplied, it should have the format ``x = f(A, b)``, where it
+            solves the equations and returns the solution. The default is
+            ``'LU'`` which corresponds to SymPy's ``A.LUsolve(b)``.
+            ``LUsolve()`` is fast to compute but will often result in
+            divide-by-zero and thus ``nan`` results.
+        **kwargs
+            Extra keyword arguments are passed to
+            :meth:`sympy.physics.mechanics.linearize.Linearizer.linearize`.
+
         Explanation
         ===========
 
@@ -383,7 +417,8 @@ def linearize(self, q_ind=None, qd_ind=None, q_dep=None, qd_dep=None,
 
         For more documentation, please see the ``Linearizer`` class."""
 
-        linearizer = self.to_linearizer(q_ind, qd_ind, q_dep, qd_dep)
+        linearizer = self.to_linearizer(q_ind, qd_ind, q_dep, qd_dep,
+                                        linear_solver=linear_solver)
         result = linearizer.linearize(**kwargs)
         return result + (linearizer.r,)
 
diff --git a/sympy/physics/mechanics/linearize.py b/sympy/physics/mechanics/linearize.py
index ca1dad2e8b..9467cd60cb 100644
--- a/sympy/physics/mechanics/linearize.py
+++ b/sympy/physics/mechanics/linearize.py
@@ -4,11 +4,12 @@
 from sympy.core.symbol import Dummy
 from sympy.utilities.iterables import flatten
 from sympy.physics.vector import dynamicsymbols
-from sympy.physics.mechanics.functions import msubs
+from sympy.physics.mechanics.functions import msubs, _parse_linear_solver
 
 from collections import namedtuple
 from collections.abc import Iterable
 
+
 class Linearizer:
     """This object holds the general model form for a dynamic system.
     This model is used for computing the linearized form of the system,
@@ -31,8 +32,9 @@ class Linearizer:
         Permutation matrix such that [q_ind, u_ind]^T = perm_mat*[q, u]^T
     """
 
-    def __init__(self, f_0, f_1, f_2, f_3, f_4, f_c, f_v, f_a, q, u,
-            q_i=None, q_d=None, u_i=None, u_d=None, r=None, lams=None):
+    def __init__(self, f_0, f_1, f_2, f_3, f_4, f_c, f_v, f_a, q, u, q_i=None,
+                 q_d=None, u_i=None, u_d=None, r=None, lams=None,
+                 linear_solver='LU'):
         """
         Parameters
         ==========
@@ -53,7 +55,19 @@ def __init__(self, f_0, f_1, f_2, f_3, f_4, f_c, f_v, f_a, q, u,
             The input variables.
         lams : array_like, optional
             The lagrange multipliers
+        linear_solver : str, callable
+            Method used to solve the several symbolic linear systems of the
+            form ``A*x=b`` in the linearization process. If a string is
+            supplied, it should be a valid method that can be used with the
+            :meth:`sympy.matrices.matrices.MatrixBase.solve`. If a callable is
+            supplied, it should have the format ``x = f(A, b)``, where it
+            solves the equations and returns the solution. The default is
+            ``'LU'`` which corresponds to SymPy's ``A.LUsolve(b)``.
+            ``LUsolve()`` is fast to compute but will often result in
+            divide-by-zero and thus ``nan`` results.
+
         """
+        self.linear_solver = _parse_linear_solver(linear_solver)
 
         # Generalized equation form
         self.f_0 = Matrix(f_0)
@@ -160,8 +174,9 @@ def _form_coefficient_matrices(self):
         # If not, C_0 is I_(nxn). Note that this works even if n=0
         if l > 0:
             f_c_jac_q = self.f_c.jacobian(self.q)
-            self._C_0 = (eye(n) - self._Pqd * (f_c_jac_q *
-                    self._Pqd).LUsolve(f_c_jac_q)) * self._Pqi
+            self._C_0 = (eye(n) - self._Pqd*
+                         self.linear_solver(f_c_jac_q*self._Pqd,
+                                            f_c_jac_q))*self._Pqi
         else:
             self._C_0 = eye(n)
         # If there are motion constraints (m > 0), form C_1 and C_2 as normal.
@@ -172,11 +187,11 @@ def _form_coefficient_matrices(self):
             temp = f_v_jac_u * self._Pud
             if n != 0:
                 f_v_jac_q = self.f_v.jacobian(self.q)
-                self._C_1 = -self._Pud * temp.LUsolve(f_v_jac_q)
+                self._C_1 = -self._Pud * self.linear_solver(temp, f_v_jac_q)
             else:
                 self._C_1 = zeros(o, n)
-            self._C_2 = (eye(o) - self._Pud *
-                    temp.LUsolve(f_v_jac_u)) * self._Pui
+            self._C_2 = (eye(o) - self._Pud*
+                         self.linear_solver(temp, f_v_jac_u))*self._Pui
         else:
             self._C_1 = zeros(o, n)
             self._C_2 = eye(o)
@@ -391,9 +406,9 @@ def linearize(self, op_point=None, A_and_B=False, simplify=False):
         # kwarg A_and_B indicates to return  A, B for forming the equation
         # dx = [A]x + [B]r, where x = [q_indnd, u_indnd]^T,
         if A_and_B:
-            A_cont = self.perm_mat.T * M_eq.LUsolve(Amat_eq)
+            A_cont = self.perm_mat.T * self.linear_solver(M_eq, Amat_eq)
             if Bmat_eq:
-                B_cont = self.perm_mat.T * M_eq.LUsolve(Bmat_eq)
+                B_cont = self.perm_mat.T * self.linear_solver(M_eq, Bmat_eq)
             else:
                 # Bmat = Matrix([]), so no need to sub
                 B_cont = Bmat_eq


Test patch:
diff --git a/sympy/physics/mechanics/tests/test_kane3.py b/sympy/physics/mechanics/tests/test_kane3.py
index a9a1011560..3da81a0f86 100644
--- a/sympy/physics/mechanics/tests/test_kane3.py
+++ b/sympy/physics/mechanics/tests/test_kane3.py
@@ -1,18 +1,44 @@
-from sympy.core.evalf import evalf
+from sympy import cacheit, zeros, ImmutableMatrix
 from sympy.core.numbers import pi
 from sympy.core.symbol import symbols
 from sympy.functions.elementary.miscellaneous import sqrt
 from sympy.functions.elementary.trigonometric import acos, sin, cos
 from sympy.matrices.dense import Matrix
 from sympy.physics.mechanics import (ReferenceFrame, dynamicsymbols,
-    KanesMethod, inertia, msubs, Point, RigidBody, dot)
-from sympy.testing.pytest import slow, ON_CI, skip
+                                     KanesMethod, inertia, Point, RigidBody,
+                                     dot)
+from sympy.testing.pytest import slow
+
+
+# TODO : This is a temporary implementation and can be removed once we support
+# a new matrix solver A.CramerSolve(b).
+@cacheit
+def det_laplace(matrix):
+    n = matrix.shape[0]
+    if n == 1:
+        return matrix[0]
+    elif n == 2:
+        return matrix[0, 0] * matrix[1, 1] - matrix[0, 1] * matrix[1, 0]
+    else:
+        return sum((-1) ** i * matrix[0, i] *
+                   det_laplace(matrix.minor_submatrix(0, i)) for i in range(n))
+
+
+def cramer_solve(A, b, det_method=det_laplace):
+    def entry(i, j):
+        return b[i, sol] if j == col else A[i, j]
+
+    A_imu = ImmutableMatrix(A)  # Convert to immutable for cache purposes
+    det_A = det_method(A_imu)
+    x = zeros(*b.shape)
+    for sol in range(b.shape[1]):
+        for col in range(b.shape[0]):
+            x[col, sol] = det_method(ImmutableMatrix(*A.shape, entry)) / det_A
+    return x
 
 
 @slow
 def test_bicycle():
-    if ON_CI:
-        skip("Too slow for CI.")
     # Code to get equations of motion for a bicycle modeled as in:
     # J.P Meijaard, Jim M Papadopoulos, Andy Ruina and A.L Schwab. Linearized
     # dynamics equations for the balance and steer of a bicycle: a benchmark
@@ -175,7 +201,8 @@ def test_bicycle():
             q_dependent=[q4], configuration_constraints=conlist_coord,
             u_ind=[u2, u3, u5],
             u_dependent=[u1, u4, u6], velocity_constraints=conlist_speed,
-            kd_eqs=kd)
+            kd_eqs=kd,
+            constraint_solver=cramer_solve)
     (fr, frstar) = KM.kanes_equations(BL, FL)
 
     # This is the start of entering in the numerical values from the benchmark
@@ -186,23 +213,23 @@ def test_bicycle():
     # model.
     PaperRadRear                    =  0.3
     PaperRadFront                   =  0.35
-    HTA                             =  evalf.N(pi / 2 - pi / 10)
+    HTA                             =  (pi / 2 - pi / 10).evalf()
     TrailPaper                      =  0.08
-    rake                            =  evalf.N(-(TrailPaper*sin(HTA)-(PaperRadFront*cos(HTA))))
+    rake                            =  (-(TrailPaper*sin(HTA)-(PaperRadFront*cos(HTA)))).evalf()
     PaperWb                         =  1.02
     PaperFrameCgX                   =  0.3
     PaperFrameCgZ                   =  0.9
     PaperForkCgX                    =  0.9
     PaperForkCgZ                    =  0.7
-    FrameLength                     =  evalf.N(PaperWb*sin(HTA)-(rake-(PaperRadFront-PaperRadRear)*cos(HTA)))
-    FrameCGNorm                     =  evalf.N((PaperFrameCgZ - PaperRadRear-(PaperFrameCgX/sin(HTA))*cos(HTA))*sin(HTA))
-    FrameCGPar                      =  evalf.N(PaperFrameCgX / sin(HTA) + (PaperFrameCgZ - PaperRadRear - PaperFrameCgX / sin(HTA) * cos(HTA)) * cos(HTA))
-    tempa                           =  evalf.N(PaperForkCgZ - PaperRadFront)
-    tempb                           =  evalf.N(PaperWb-PaperForkCgX)
-    tempc                           =  evalf.N(sqrt(tempa**2+tempb**2))
-    PaperForkL                      =  evalf.N(PaperWb*cos(HTA)-(PaperRadFront-PaperRadRear)*sin(HTA))
-    ForkCGNorm                      =  evalf.N(rake+(tempc * sin(pi/2-HTA-acos(tempa/tempc))))
-    ForkCGPar                       =  evalf.N(tempc * cos((pi/2-HTA)-acos(tempa/tempc))-PaperForkL)
+    FrameLength                     =  (PaperWb*sin(HTA)-(rake-(PaperRadFront-PaperRadRear)*cos(HTA))).evalf()
+    FrameCGNorm                     =  ((PaperFrameCgZ - PaperRadRear-(PaperFrameCgX/sin(HTA))*cos(HTA))*sin(HTA)).evalf()
+    FrameCGPar                      =  (PaperFrameCgX / sin(HTA) + (PaperFrameCgZ - PaperRadRear - PaperFrameCgX / sin(HTA) * cos(HTA)) * cos(HTA)).evalf()
+    tempa                           =  (PaperForkCgZ - PaperRadFront)
+    tempb                           =  (PaperWb-PaperForkCgX)
+    tempc                           =  (sqrt(tempa**2+tempb**2)).evalf()
+    PaperForkL                      =  (PaperWb*cos(HTA)-(PaperRadFront-PaperRadRear)*sin(HTA)).evalf()
+    ForkCGNorm                      =  (rake+(tempc * sin(pi/2-HTA-acos(tempa/tempc)))).evalf()
+    ForkCGPar                       =  (tempc * cos((pi/2-HTA)-acos(tempa/tempc))-PaperForkL).evalf()
 
     # Here is the final assembly of the numerical values. The symbol 'v' is the
     # forward speed of the bicycle (a concept which only makes sense in the
@@ -257,27 +284,47 @@ def test_bicycle():
     # many rows as *total* coordinates and speeds, but only as many columns as
     # independent coordinates and speeds.
 
-    forcing_lin = KM.linearize()[0]
-
+    A, B, _ = KM.linearize(
+        A_and_B=True,
+        op_point={
+            # Operating points for the accelerations are required for the
+            # linearizer to eliminate u' terms showing up in the coefficient
+            # matrices.
+            u1.diff(): 0,
+            u2.diff(): 0,
+            u3.diff(): 0,
+            u4.diff(): 0,
+            u5.diff(): 0,
+            u6.diff(): 0,
+            u1: 0,
+            u2: 0,
+            u3: v / PaperRadRear,
+            u4: 0,
+            u5: 0,
+            u6: v / PaperRadFront,
+            q1: 0,
+            q2: 0,
+            q4: 0,
+            q5: 0,
+        },
+        linear_solver=cramer_solve,
+    )
     # As mentioned above, the size of the linearized forcing terms is expanded
     # to include both q's and u's, so the mass matrix must have this done as
     # well.  This will likely be changed to be part of the linearized process,
     # for future reference.
-    MM_full = KM.mass_matrix_full
-
-    MM_full_s = msubs(MM_full, val_dict)
-    forcing_lin_s = msubs(forcing_lin, KM.kindiffdict(), val_dict)
+    A_s = A.xreplace(val_dict)
+    B_s = B.xreplace(val_dict)
 
-    MM_full_s = MM_full_s.evalf()
-    forcing_lin_s = forcing_lin_s.evalf()
+    A_s = A_s.evalf()
+    B_s = B_s.evalf()
 
     # Finally, we construct an "A" matrix for the form xdot = A x (x being the
     # state vector, although in this case, the sizes are a little off). The
     # following line extracts only the minimum entries required for eigenvalue
     # analysis, which correspond to rows and columns for lean, steer, lean
     # rate, and steer rate.
-    Amat = MM_full_s.inv() * forcing_lin_s
-    A = Amat.extract([1, 2, 4, 6], [1, 2, 3, 5])
+    A = A_s.extract([1, 2, 3, 5], [1, 2, 3, 5])
 
     # Precomputed for comparison
     Res = Matrix([[               0,                                           0,                  1.0,                    0],
@@ -285,7 +332,6 @@ def test_bicycle():
                   [9.48977444677355, -0.891197738059089*v**2 - 0.571523173729245, -0.105522449805691*v, -0.330515398992311*v],
                   [11.7194768719633,   -1.97171508499972*v**2 + 30.9087533932407,   3.67680523332152*v,  -3.08486552743311*v]])
 
-
     # Actual eigenvalue comparison
     eps = 1.e-12
     for i in range(6):
diff --git a/sympy/physics/mechanics/tests/test_linearize.py b/sympy/physics/mechanics/tests/test_linearize.py
index 1c9c2aeed6..ff4ae07831 100644
--- a/sympy/physics/mechanics/tests/test_linearize.py
+++ b/sympy/physics/mechanics/tests/test_linearize.py
@@ -1,5 +1,5 @@
 from sympy.core.backend import (symbols, Matrix, cos, sin, atan, sqrt,
-    Rational, _simplify_matrix)
+    Rational, _simplify_matrix, USE_SYMENGINE)
 from sympy.core.sympify import sympify
 from sympy.simplify.simplify import simplify
 from sympy.solvers.solvers import solve
@@ -131,6 +131,15 @@ def test_linearize_rolling_disc_kane():
     # Check eigenvalues at critical speed are all zero:
     assert sympify(A.subs(upright_nominal).subs(q3d, 1/sqrt(3))).eigenvals() == {0: 8}
 
+    # Check whether alternative solvers work
+    if not USE_SYMENGINE:
+        # symengine doesn't support method='GJ'
+        linearizer = KM.to_linearizer(linear_solver='GJ')
+        A, B = linearizer.linearize(op_point=[q_op, u_op, qd_op, ud_op],
+                                    A_and_B=True, simplify=True)
+        assert A.subs(upright_nominal) == A_sol
+        assert B.subs(upright_nominal) == B_sol
+
 def test_linearize_pendulum_kane_minimal():
     q1 = dynamicsymbols('q1')                     # angle of pendulum
     u1 = dynamicsymbols('u1')                     # Angular velocity
@@ -232,6 +241,24 @@ def test_linearize_pendulum_kane_nonminimal():
     assert A.expand() == Matrix([[0, 1], [-9.8/L, 0]])
     assert B == Matrix([])
 
+
+    if not USE_SYMENGINE:
+        # symengine doesn't support method='GJ'
+        A, B, inp_vec = KM.linearize(op_point=[q_op, u_op, ud_op], A_and_B=True,
+                                    simplify=True, linear_solver='GJ')
+
+        assert A.expand() == Matrix([[0, 1], [-9.8/L, 0]])
+        assert B == Matrix([])
+
+    A, B, inp_vec = KM.linearize(op_point=[q_op, u_op, ud_op],
+                                 A_and_B=True,
+                                 simplify=True,
+                                 linear_solver=lambda A, b: A.LUsolve(b))
+
+    assert A.expand() == Matrix([[0, 1], [-9.8/L, 0]])
+    assert B == Matrix([])
+
+
 def test_linearize_pendulum_lagrange_minimal():
     q1 = dynamicsymbols('q1')                     # angle of pendulum
     q1d = dynamicsymbols('q1', 1)                 # Angular velocity
@@ -263,6 +290,16 @@ def test_linearize_pendulum_lagrange_minimal():
     assert _simplify_matrix(A) == Matrix([[0, 1], [-9.8*cos(q1)/L, 0]])
     assert B == Matrix([])
 
+    # Check an alternative solver
+    if not USE_SYMENGINE:
+        # symengine doesn't support method='GJ'
+        A, B, inp_vec = LM.linearize([q1], [q1d], A_and_B=True,
+                                     linear_solver='GJ')
+
+        assert _simplify_matrix(A) == Matrix([[0, 1], [-9.8*cos(q1)/L, 0]])
+        assert B == Matrix([])
+
+
 def test_linearize_pendulum_lagrange_nonminimal():
     q1, q2 = dynamicsymbols('q1:3')
     q1d, q2d = dynamicsymbols('q1:3', level=1)
@@ -296,6 +333,13 @@ def test_linearize_pendulum_lagrange_nonminimal():
     assert _simplify_matrix(A) == Matrix([[0, 1], [-9.8/L, 0]])
     assert B == Matrix([])
 
+    # Check if passing a function to linear_solver works
+    A, B, inp_vec = LM.linearize([q2], [q2d], [q1], [q1d], op_point=op_point,
+                                 A_and_B=True, linear_solver=lambda A, b:
+                                 A.LUsolve(b))
+    assert _simplify_matrix(A) == Matrix([[0, 1], [-9.8/L, 0]])
+    assert B == Matrix([])
+
 def test_linearize_rolling_disc_lagrange():
     q1, q2, q3 = q = dynamicsymbols('q1 q2 q3')
     q1d, q2d, q3d = qd = dynamicsymbols('q1 q2 q3', 1)