Repository: sympy/sympy

Issue:
KanesMethod numerical error due to LUsolve of the velocity constraints
`KanesMethod` uses and [`LUsolve` ](https://github.com/sympy/sympy/blob/4fb6b42dd0941f05178aa570edf9ef61162ce960/sympy/physics/mechanics/kane.py#L252) to solve the velocity constraints. While this matrix is invertible it does cause numerical instabilities from time to time, caused by zero divisions. This problem is illustrated in the [following notebook](https://github.com/sympy/sympy/files/10825869/whipple_bicycle.zip), which models a Whipple bicycle model as proposed in [this comment](https://github.com/sympy/sympy/issues/24679#issuecomment-1443682326).

While it can be seen that [this solution](https://github.com/sympy/sympy/issues/24679#issuecomment-1424672119) proposed by @oscarbenjamin does result in numerical stable solutions, it does increase the number of operations massively (see [this comment](https://github.com/sympy/sympy/issues/24679#issuecomment-1443413703)).

Another possible solution is to use [`MatrixSolve`](https://github.com/sympy/sympy/blob/4fb6b42dd0941f05178aa570edf9ef61162ce960/sympy/codegen/matrix_nodes.py#L27), which would move the solve towards the codegen. This would however result in non-differentiable equations of motion, which is undesirable in several optimization algorithms. A second slight disadvantage of using `MatrixSolve` would be the requirement to use `BlockMatrix` when creating for example the `KanesMethod.mass_matrix`.


Conclusion from below on what needs to be done:
- [x] Support the use of different solvers in `KanesMethod` for solving the constraints and the kinematic differential equations.
- [x] Support the use of different linear solvers in `Linearize`
- [ ] Implement `cramer_solve`.
- [ ] Decide on the default determinant of `cramer_solve`
- [ ] Add support for a numerical solver in `KanesMethod`

Discussion:

Patch:
diff --git a/sympy/matrices/determinant.py b/sympy/matrices/determinant.py
index 7986f28a04..3ebeee2e98 100644
--- a/sympy/matrices/determinant.py
+++ b/sympy/matrices/determinant.py
@@ -1,5 +1,6 @@
 from types import FunctionType
 
+from sympy.core.cache import cacheit
 from sympy.core.numbers import Float, Integer
 from sympy.core.singleton import S
 from sympy.core.symbol import uniquely_named_symbol
@@ -7,6 +8,7 @@
 from sympy.polys import PurePoly, cancel
 from sympy.functions.combinatorial.numbers import nC
 from sympy.polys.matrices.domainmatrix import DomainMatrix
+from sympy.polys.matrices.ddm import DDM
 
 from .common import NonSquareMatrixError
 from .utilities import (
@@ -303,8 +305,8 @@ def _adjugate(M, method="berkowitz"):
     ==========
 
     method : string, optional
-        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
-        "lu".
+        Method to use to find the cofactors, can be "bareiss", "berkowitz",
+        "bird", "laplace" or "lu".
 
     Examples
     ========
@@ -421,8 +423,8 @@ def _cofactor(M, i, j, method="berkowitz"):
     ==========
 
     method : string, optional
-        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
-        "lu".
+        Method to use to find the cofactors, can be "bareiss", "berkowitz",
+        "bird", "laplace" or "lu".
 
     Examples
     ========
@@ -453,8 +455,8 @@ def _cofactor_matrix(M, method="berkowitz"):
     ==========
 
     method : string, optional
-        Method to use to find the cofactors, can be "bareiss", "berkowitz" or
-        "lu".
+        Method to use to find the cofactors, can be "bareiss", "berkowitz",
+        "bird", "laplace" or "lu".
 
     Examples
     ========
@@ -569,6 +571,10 @@ def _det(M, method="bareiss", iszerofunc=None):
 
         If it is set to ``'berkowitz'``, Berkowitz' algorithm will be used.
 
+        If it is set to ``'bird'``, Bird's algorithm will be used [1]_.
+
+        If it is set to ``'laplace'``, Laplace's algorithm will be used.
+
         Otherwise, if it is set to ``'lu'``, LU decomposition will be used.
 
         .. note::
@@ -617,6 +623,13 @@ def _det(M, method="bareiss", iszerofunc=None):
     True
     >>> M.det(method="domain-ge")
     -2
+
+    References
+    ==========
+
+    .. [1] Bird, R. S. (2011). A simple division-free algorithm for computing
+           determinants. Inf. Process. Lett., 111(21), 1072-1074. doi:
+           10.1016/j.ipl.2011.08.006
     """
 
     # sanitize `method`
@@ -627,7 +640,8 @@ def _det(M, method="bareiss", iszerofunc=None):
     elif method == "det_lu":
         method = "lu"
 
-    if method not in ("bareiss", "berkowitz", "lu", "domain-ge"):
+    if method not in ("bareiss", "berkowitz", "lu", "domain-ge", "bird",
+                      "laplace"):
         raise ValueError("Determinant method '%s' unrecognized" % method)
 
     if iszerofunc is None:
@@ -668,6 +682,10 @@ def _det(M, method="bareiss", iszerofunc=None):
             det = M[b, b]._eval_det_berkowitz()
         elif method == "lu":
             det = M[b, b]._eval_det_lu(iszerofunc=iszerofunc)
+        elif method == "bird":
+            det = M[b, b]._eval_det_bird()
+        elif method == "laplace":
+            det = M[b, b]._eval_det_laplace()
         dets.append(det)
     return Mul(*dets)
 
@@ -819,6 +837,89 @@ def _det_LU(M, iszerofunc=_iszero, simpfunc=None):
     return det
 
 
+@cacheit
+def __det_laplace(M):
+    """Compute the determinant of a matrix using Laplace expansion.
+
+    This is a recursive function, and it should not be called directly.
+    Use _det_laplace() instead. The reason for splitting this function
+    into two is to allow caching of determinants of submatrices. While
+    one could also define this function inside _det_laplace(), that
+    would remove the advantage of using caching in Cramer Solve.
+    """
+    n = M.shape[0]
+    if n == 1:
+        return M[0]
+    elif n == 2:
+        return M[0, 0] * M[1, 1] - M[0, 1] * M[1, 0]
+    else:
+        return sum((-1) ** i * M[0, i] *
+                   __det_laplace(M.minor_submatrix(0, i)) for i in range(n))
+
+
+def _det_laplace(M):
+    """Compute the determinant of a matrix using Laplace expansion.
+
+    While Laplace expansion is not the most efficient method of computing
+    a determinant, it is a simple one, and it has the advantage of
+    being division free. To improve efficiency, this function uses
+    caching to avoid recomputing determinants of submatrices.
+    """
+    if not M.is_square:
+        raise NonSquareMatrixError()
+    if M.shape[0] == 0:
+        return M.one
+        # sympy/matrices/tests/test_matrices.py contains a test that
+        # suggests that the determinant of a 0 x 0 matrix is one, by
+        # convention.
+    return __det_laplace(M.as_immutable())
+
+
+def _det_bird(M):
+    r"""Compute the determinant of a matrix using Bird's algorithm.
+
+    Bird's algorithm is a simple division-free algorithm for computing, which
+    is of lower order than the Laplace's algorithm. It is described in [1]_.
+
+    References
+    ==========
+
+    .. [1] Bird, R. S. (2011). A simple division-free algorithm for computing
+           determinants. Inf. Process. Lett., 111(21), 1072-1074. doi:
+           10.1016/j.ipl.2011.08.006
+    """
+    def mu(X):
+        n = X.shape[0]
+        zero = X.domain.zero
+
+        total = zero
+        diag_sums = [zero]
+        for i in reversed(range(1, n)):
+            total -= X[i][i]
+            diag_sums.append(total)
+        diag_sums = diag_sums[::-1]
+
+        elems = [[zero] * i + [diag_sums[i]] + X_i[i + 1:] for i, X_i in
+                 enumerate(X)]
+        return DDM(elems, X.shape, X.domain)
+
+    Mddm = M._rep.to_dense().rep
+    n = M.shape[0]
+    if n == 0:
+        return M.one
+        # sympy/matrices/tests/test_matrices.py contains a test that
+        # suggests that the determinant of a 0 x 0 matrix is one, by
+        # convention.
+    Fn1 = Mddm
+    for _ in range(n - 1):
+        Fn1 = mu(Fn1).matmul(Mddm)
+    detA = Fn1[0][0]
+    if n % 2 == 0:
+        detA = -detA
+
+    return Mddm.domain.to_sympy(detA)
+
+
 def _minor(M, i, j, method="berkowitz"):
     """Return the (i,j) minor of ``M``.  That is,
     return the determinant of the matrix obtained by deleting
@@ -832,7 +933,7 @@ def _minor(M, i, j, method="berkowitz"):
 
     method : string, optional
         Method to use to find the determinant of the submatrix, can be
-        "bareiss", "berkowitz" or "lu".
+        "bareiss", "berkowitz", "bird", "laplace" or "lu".
 
     Examples
     ========
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index 9254c5773d..f68333306b 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -33,7 +33,8 @@
 from .determinant import (
     _find_reasonable_pivot, _find_reasonable_pivot_naive,
     _adjugate, _charpoly, _cofactor, _cofactor_matrix, _per,
-    _det, _det_bareiss, _det_berkowitz, _det_LU, _minor, _minor_submatrix)
+    _det, _det_bareiss, _det_berkowitz, _det_bird, _det_laplace, _det_LU,
+    _minor, _minor_submatrix)
 
 from .reductions import _is_echelon, _echelon_form, _rank, _rref
 from .subspaces import _columnspace, _nullspace, _rowspace, _orthogonalize
@@ -58,7 +59,7 @@
 from .solvers import (
     _diagonal_solve, _lower_triangular_solve, _upper_triangular_solve,
     _cholesky_solve, _LDLsolve, _LUsolve, _QRsolve, _gauss_jordan_solve,
-    _pinv_solve, _solve, _solve_least_squares)
+    _pinv_solve, _cramer_solve, _solve, _solve_least_squares)
 
 from .inverse import (
     _pinv, _inv_mod, _inv_ADJ, _inv_GE, _inv_LU, _inv_CH, _inv_LDL, _inv_QR,
@@ -109,6 +110,12 @@ def _eval_det_berkowitz(self):
     def _eval_det_lu(self, iszerofunc=_iszero, simpfunc=None):
         return _det_LU(self, iszerofunc=iszerofunc, simpfunc=simpfunc)
 
+    def _eval_det_bird(self):
+        return _det_bird(self)
+
+    def _eval_det_laplace(self):
+        return _det_laplace(self)
+
     def _eval_determinant(self): # for expressions.determinant.Determinant
         return _det(self)
 
@@ -140,6 +147,8 @@ def minor_submatrix(self, i, j):
     _find_reasonable_pivot_naive.__doc__ = _find_reasonable_pivot_naive.__doc__
     _eval_det_bareiss.__doc__            = _det_bareiss.__doc__
     _eval_det_berkowitz.__doc__          = _det_berkowitz.__doc__
+    _eval_det_bird.__doc__            = _det_bird.__doc__
+    _eval_det_laplace.__doc__            = _det_laplace.__doc__
     _eval_det_lu.__doc__                 = _det_LU.__doc__
     _eval_determinant.__doc__            = _det.__doc__
     adjugate.__doc__                     = _adjugate.__doc__
@@ -2274,6 +2283,9 @@ def gauss_jordan_solve(self, B, freevar=False):
     def pinv_solve(self, B, arbitrary_matrix=None):
         return _pinv_solve(self, B, arbitrary_matrix=arbitrary_matrix)
 
+    def cramer_solve(self, rhs, det_method="laplace"):
+        return _cramer_solve(self, rhs, det_method=det_method)
+
     def solve(self, rhs, method='GJ'):
         return _solve(self, rhs, method=method)
 
@@ -2344,6 +2356,7 @@ def strongly_connected_components_decomposition(self, lower=True):
     QRsolve.__doc__                = _QRsolve.__doc__
     gauss_jordan_solve.__doc__     = _gauss_jordan_solve.__doc__
     pinv_solve.__doc__             = _pinv_solve.__doc__
+    cramer_solve.__doc__           = _cramer_solve.__doc__
     solve.__doc__                  = _solve.__doc__
     solve_least_squares.__doc__    = _solve_least_squares.__doc__
 
diff --git a/sympy/matrices/solvers.py b/sympy/matrices/solvers.py
index ddfe4473ef..8cb6133c9a 100644
--- a/sympy/matrices/solvers.py
+++ b/sympy/matrices/solvers.py
@@ -31,6 +31,7 @@ def _diagonal_solve(M, rhs):
     LUsolve
     QRsolve
     pinv_solve
+    cramer_solve
     """
 
     if not M.is_diagonal():
@@ -56,6 +57,7 @@ def _lower_triangular_solve(M, rhs):
     LUsolve
     QRsolve
     pinv_solve
+    cramer_solve
     """
 
     from .dense import MutableDenseMatrix
@@ -94,6 +96,7 @@ def _lower_triangular_solve_sparse(M, rhs):
     LUsolve
     QRsolve
     pinv_solve
+    cramer_solve
     """
 
     if not M.is_square:
@@ -136,6 +139,7 @@ def _upper_triangular_solve(M, rhs):
     LUsolve
     QRsolve
     pinv_solve
+    cramer_solve
     """
 
     from .dense import MutableDenseMatrix
@@ -174,6 +178,7 @@ def _upper_triangular_solve_sparse(M, rhs):
     LUsolve
     QRsolve
     pinv_solve
+    cramer_solve
     """
 
     if not M.is_square:
@@ -219,6 +224,7 @@ def _cholesky_solve(M, rhs):
     LUsolve
     QRsolve
     pinv_solve
+    cramer_solve
     """
 
     if M.rows < M.cols:
@@ -276,6 +282,7 @@ def _LDLsolve(M, rhs):
     LUsolve
     QRsolve
     pinv_solve
+    cramer_solve
     """
 
     if M.rows < M.cols:
@@ -324,6 +331,7 @@ def _LUsolve(M, rhs, iszerofunc=_iszero):
     QRsolve
     pinv_solve
     LUdecomposition
+    cramer_solve
     """
 
     if rhs.rows != M.rows:
@@ -400,6 +408,7 @@ def _QRsolve(M, b):
     LUsolve
     pinv_solve
     QRdecomposition
+    cramer_solve
     """
 
     dps  = _get_intermediate_simp(expand_mul, expand_mul)
@@ -699,6 +708,73 @@ def _pinv_solve(M, B, arbitrary_matrix=None):
             A_pinv.multiply(A)).multiply(arbitrary_matrix)
 
 
+def _cramer_solve(M, rhs, det_method="laplace"):
+    """Solves system of linear equations using Cramer's rule.
+
+    This method is relatively inefficient compared to other methods.
+    However it only uses a single division, assuming a division-free determinant
+    method is provided. This is helpful to minimize the chance of divide-by-zero
+    cases in symbolic solutions to linear systems.
+
+    Parameters
+    ==========
+    M : Matrix
+        The matrix representing the left hand side of the equation.
+    rhs : Matrix
+        The matrix representing the right hand side of the equation.
+    det_method : str or callable
+        The method to use to calculate the determinant of the matrix.
+        The default is ``'laplace'``.  If a callable is passed, it should take a
+        single argument, the matrix, and return the determinant of the matrix.
+
+    Returns
+    =======
+    x : Matrix
+        The matrix that will satisfy ``Ax = B``.  Will have as many rows as
+        matrix A has columns, and as many columns as matrix B.
+
+    Examples
+    ========
+
+    >>> from sympy import Matrix
+    >>> A = Matrix([[0, -6, 1], [0, -6, -1], [-5, -2, 3]])
+    >>> B = Matrix([[-30, -9], [-18, -27], [-26, 46]])
+    >>> x = A.cramer_solve(B)
+    >>> x
+    Matrix([
+    [ 0, -5],
+    [ 4,  3],
+    [-6,  9]])
+
+    References
+    ==========
+
+    .. [1] https://en.wikipedia.org/wiki/Cramer%27s_rule#Explicit_formulas_for_small_systems
+
+    """
+    from .dense import zeros
+
+    def entry(i, j):
+        return rhs[i, sol] if j == col else M[i, j]
+
+    if det_method == "bird":
+        from .determinant import _det_bird
+        det = _det_bird
+    elif det_method == "laplace":
+        from .determinant import _det_laplace
+        det = _det_laplace
+    elif isinstance(det_method, str):
+        det = lambda matrix: matrix.det(method=det_method)
+    else:
+        det = det_method
+    det_M = det(M)
+    x = zeros(*rhs.shape)
+    for sol in range(rhs.shape[1]):
+        for col in range(rhs.shape[0]):
+            x[col, sol] = det(M.__class__(*M.shape, entry)) / det_M
+    return M.__class__(x)
+
+
 def _solve(M, rhs, method='GJ'):
     """Solves linear equation where the unique solution exists.
 
@@ -718,6 +794,8 @@ def _solve(M, rhs, method='GJ'):
 
         If set to ``'PINV'``, ``pinv_solve`` routine will be used.
 
+        If set to ``'CRAMER'``, ``cramer_solve`` routine will be used.
+
         It also supports the methods available for special linear systems
 
         For positive definite systems:
@@ -769,6 +847,8 @@ def _solve(M, rhs, method='GJ'):
         return M.LDLsolve(rhs)
     elif method == 'PINV':
         return M.pinv_solve(rhs)
+    elif method == 'CRAMER':
+        return M.cramer_solve(rhs)
     else:
         return M.inv(method=method).multiply(rhs)
 
diff --git a/sympy/physics/mechanics/kane.py b/sympy/physics/mechanics/kane.py
index e1e126a6ec..abeb7b0de2 100644
--- a/sympy/physics/mechanics/kane.py
+++ b/sympy/physics/mechanics/kane.py
@@ -90,6 +90,11 @@ class KanesMethod(_Methods):
     number of operations. The weakness of this method is that it can result in zero
     division errors.
 
+    If zero divisions are encountered, a possible solver which may solve the problem
+    is ``"CRAMER"``. This method uses Cramer's rule to solve the system. This method
+    is slower and results in more operations than the default solver. However it only
+    uses a single division by default per entry of the solution.
+
     While a valid list of solvers can be found at
     :meth:`sympy.matrices.matrices.MatrixBase.solve`, it is also possible to supply a
     `callable`. This way it is possible to use a different solver routine. If the


Test patch:
diff --git a/sympy/matrices/tests/test_determinant.py b/sympy/matrices/tests/test_determinant.py
index 7690e57bd6..4c7f56dd6b 100644
--- a/sympy/matrices/tests/test_determinant.py
+++ b/sympy/matrices/tests/test_determinant.py
@@ -1,4 +1,5 @@
 import random
+import pytest
 from sympy.core.numbers import I
 from sympy.core.numbers import Rational
 from sympy.core.symbol import (Symbol, symbols)
@@ -11,159 +12,80 @@
 from sympy.functions.combinatorial.factorials import factorial, subfactorial
 
 
-def test_determinant():
-    M = Matrix()
-
-    assert M.det() == 1
+@pytest.mark.parametrize("method", [
     # Evaluating these directly because they are never reached via M.det()
-    assert M._eval_det_bareiss() == 1
-    assert M._eval_det_berkowitz() == 1
-    assert M._eval_det_lu() == 1
-
-    M = Matrix([ [0] ])
-
-    assert M.det() == 0
-    assert M._eval_det_bareiss() == 0
-    assert M._eval_det_berkowitz() == 0
-    assert M._eval_det_lu() == 0
-
-    M = Matrix([ [5] ])
-
-    assert M.det() == 5
-    assert M._eval_det_bareiss() == 5
-    assert M._eval_det_berkowitz() == 5
-    assert M._eval_det_lu() == 5
-
-    M = Matrix(( (-3,  2),
-                 ( 8, -5) ))
-
-    assert M.det(method="domain-ge") == -1
-    assert M.det(method="bareiss") == -1
-    assert M.det(method="berkowitz") == -1
-    assert M.det(method="lu") == -1
-
-    M = Matrix(( (x,   1),
-                 (y, 2*y) ))
-
-    assert M.det(method="domain-ge") == 2*x*y - y
-    assert M.det(method="bareiss") == 2*x*y - y
-    assert M.det(method="berkowitz") == 2*x*y - y
-    assert M.det(method="lu") == 2*x*y - y
-
-    M = Matrix(( (1, 1, 1),
-                 (1, 2, 3),
-                 (1, 3, 6) ))
-
-    assert M.det(method="domain-ge") == 1
-    assert M.det(method="bareiss") == 1
-    assert M.det(method="berkowitz") == 1
-    assert M.det(method="lu") == 1
-
-    M = Matrix(( ( 3, -2,  0, 5),
-                 (-2,  1, -2, 2),
-                 ( 0, -2,  5, 0),
-                 ( 5,  0,  3, 4) ))
-
-    assert M.det(method="domain-ge") == -289
-    assert M.det(method="bareiss") == -289
-    assert M.det(method="berkowitz") == -289
-    assert M.det(method="lu") == -289
-
-    M = Matrix(( ( 1,  2,  3,  4),
-                 ( 5,  6,  7,  8),
-                 ( 9, 10, 11, 12),
-                 (13, 14, 15, 16) ))
-
-    assert M.det(method="domain-ge") == 0
-    assert M.det(method="bareiss") == 0
-    assert M.det(method="berkowitz") == 0
-    assert M.det(method="lu") == 0
-
-    M = Matrix(( (3, 2, 0, 0, 0),
-                 (0, 3, 2, 0, 0),
-                 (0, 0, 3, 2, 0),
-                 (0, 0, 0, 3, 2),
-                 (2, 0, 0, 0, 3) ))
-
-    assert M.det(method="domain-ge") == 275
-    assert M.det(method="bareiss") == 275
-    assert M.det(method="berkowitz") == 275
-    assert M.det(method="lu") == 275
-
-    M = Matrix(( ( 3,  0,  0, 0),
-                 (-2,  1,  0, 0),
-                 ( 0, -2,  5, 0),
-                 ( 5,  0,  3, 4) ))
-
-    assert M.det(method="domain-ge") == 60
-    assert M.det(method="bareiss") == 60
-    assert M.det(method="berkowitz") == 60
-    assert M.det(method="lu") == 60
-
-    M = Matrix(( ( 1,  0,  0,  0),
-                 ( 5,  0,  0,  0),
-                 ( 9, 10, 11, 0),
-                 (13, 14, 15, 16) ))
-
-    assert M.det(method="domain-ge") == 0
-    assert M.det(method="bareiss") == 0
-    assert M.det(method="berkowitz") == 0
-    assert M.det(method="lu") == 0
-
-    M = Matrix(( (3, 2, 0, 0, 0),
-                 (0, 3, 2, 0, 0),
-                 (0, 0, 3, 2, 0),
-                 (0, 0, 0, 3, 2),
-                 (0, 0, 0, 0, 3) ))
-
-    assert M.det(method="domain-ge") == 243
-    assert M.det(method="bareiss") == 243
-    assert M.det(method="berkowitz") == 243
-    assert M.det(method="lu") == 243
-
-    M = Matrix(( (1, 0,  1,  2, 12),
-                 (2, 0,  1,  1,  4),
-                 (2, 1,  1, -1,  3),
-                 (3, 2, -1,  1,  8),
-                 (1, 1,  1,  0,  6) ))
-
-    assert M.det(method="domain-ge") == -55
-    assert M.det(method="bareiss") == -55
-    assert M.det(method="berkowitz") == -55
-    assert M.det(method="lu") == -55
-
-    M = Matrix(( (-5,  2,  3,  4,  5),
-                 ( 1, -4,  3,  4,  5),
-                 ( 1,  2, -3,  4,  5),
-                 ( 1,  2,  3, -2,  5),
-                 ( 1,  2,  3,  4, -1) ))
-
-    assert M.det(method="domain-ge") == 11664
-    assert M.det(method="bareiss") == 11664
-    assert M.det(method="berkowitz") == 11664
-    assert M.det(method="lu") == 11664
-
-    M = Matrix(( ( 2,  7, -1, 3, 2),
-                 ( 0,  0,  1, 0, 1),
-                 (-2,  0,  7, 0, 2),
-                 (-3, -2,  4, 5, 3),
-                 ( 1,  0,  0, 0, 1) ))
-
-    assert M.det(method="domain-ge") == 123
-    assert M.det(method="bareiss") == 123
-    assert M.det(method="berkowitz") == 123
-    assert M.det(method="lu") == 123
-
-    M = Matrix(( (x, y, z),
-                 (1, 0, 0),
-                 (y, z, x) ))
-
-    assert M.det(method="domain-ge") == z**2 - x*y
-    assert M.det(method="bareiss") == z**2 - x*y
-    assert M.det(method="berkowitz") == z**2 - x*y
-    assert M.det(method="lu") == z**2 - x*y
-
-    # issue 13835
+    Matrix._eval_det_bareiss, Matrix._eval_det_berkowitz,
+    Matrix._eval_det_bird, Matrix._eval_det_laplace, Matrix._eval_det_lu
+])
+@pytest.mark.parametrize("M, sol", [
+    (Matrix(), 1),
+    (Matrix([[0]]), 0),
+    (Matrix([[5]]), 5),
+])
+def test_eval_determinant(method, M, sol):
+    assert method(M) == sol
+
+
+@pytest.mark.parametrize("method", [
+    "domain-ge", "bareiss", "berkowitz", "bird", "laplace", "lu"])
+@pytest.mark.parametrize("M, sol", [
+    (Matrix(( (-3,  2),
+              ( 8, -5) )), -1),
+    (Matrix(( (x,   1),
+              (y, 2*y) )), 2*x*y - y),
+    (Matrix(( (1, 1, 1),
+              (1, 2, 3),
+              (1, 3, 6) )), 1),
+    (Matrix(( ( 3, -2,  0, 5),
+              (-2,  1, -2, 2),
+              ( 0, -2,  5, 0),
+              ( 5,  0,  3, 4) )), -289),
+    (Matrix(( ( 1,  2,  3,  4),
+              ( 5,  6,  7,  8),
+              ( 9, 10, 11, 12),
+              (13, 14, 15, 16) )), 0),
+    (Matrix(( (3, 2, 0, 0, 0),
+              (0, 3, 2, 0, 0),
+              (0, 0, 3, 2, 0),
+              (0, 0, 0, 3, 2),
+              (2, 0, 0, 0, 3) )), 275),
+    (Matrix(( ( 3,  0,  0, 0),
+              (-2,  1,  0, 0),
+              ( 0, -2,  5, 0),
+              ( 5,  0,  3, 4) )), 60),
+    (Matrix(( ( 1,  0,  0,  0),
+              ( 5,  0,  0,  0),
+              ( 9, 10, 11, 0),
+              (13, 14, 15, 16) )), 0),
+    (Matrix(( (3, 2, 0, 0, 0),
+              (0, 3, 2, 0, 0),
+              (0, 0, 3, 2, 0),
+              (0, 0, 0, 3, 2),
+              (0, 0, 0, 0, 3) )), 243),
+    (Matrix(( (1, 0,  1,  2, 12),
+              (2, 0,  1,  1,  4),
+              (2, 1,  1, -1,  3),
+              (3, 2, -1,  1,  8),
+              (1, 1,  1,  0,  6) )), -55),
+    (Matrix(( (-5,  2,  3,  4,  5),
+              ( 1, -4,  3,  4,  5),
+              ( 1,  2, -3,  4,  5),
+              ( 1,  2,  3, -2,  5),
+              ( 1,  2,  3,  4, -1) )), 11664),
+    (Matrix(( ( 2,  7, -1, 3, 2),
+              ( 0,  0,  1, 0, 1),
+              (-2,  0,  7, 0, 2),
+              (-3, -2,  4, 5, 3),
+              ( 1,  0,  0, 0, 1) )), 123),
+    (Matrix(( (x, y, z),
+              (1, 0, 0),
+              (y, z, x) )), z**2 - x*y),
+])
+def test_determinant(method, M, sol):
+    assert M.det(method=method) == sol
+
+
+def test_issue_13835():
     a = symbols('a')
     M = lambda n: Matrix([[i + a*j for i in range(n)]
                           for j in range(n)])
@@ -171,6 +93,7 @@ def test_determinant():
     assert M(6).det() == 0
     assert M(7).det() == 0
 
+
 def test_issue_14517():
     M = Matrix([
         [   0, 10*I,    10*I,       0],
@@ -182,86 +105,25 @@ def test_issue_14517():
     test_ev = random.choice(list(ev.keys()))
     assert (M - test_ev*eye(4)).det() == 0
 
-def test_legacy_det():
+
+@pytest.mark.parametrize("method", [
+    "bareis", "det_lu", "det_LU", "Bareis", "BAREISS", "BERKOWITZ", "LU"])
+@pytest.mark.parametrize("M, sol", [
+    (Matrix(( ( 3, -2,  0, 5),
+              (-2,  1, -2, 2),
+              ( 0, -2,  5, 0),
+              ( 5,  0,  3, 4) )), -289),
+    (Matrix(( (-5,  2,  3,  4,  5),
+              ( 1, -4,  3,  4,  5),
+              ( 1,  2, -3,  4,  5),
+              ( 1,  2,  3, -2,  5),
+              ( 1,  2,  3,  4, -1) )), 11664),
+])
+def test_legacy_det(method, M, sol):
     # Minimal support for legacy keys for 'method' in det()
     # Partially copied from test_determinant()
+    assert M.det(method=method) == sol
 
-    M = Matrix(( ( 3, -2,  0, 5),
-                 (-2,  1, -2, 2),
-                 ( 0, -2,  5, 0),
-                 ( 5,  0,  3, 4) ))
-
-    assert M.det(method="bareis") == -289
-    assert M.det(method="det_lu") == -289
-    assert M.det(method="det_LU") == -289
-
-    M = Matrix(( (3, 2, 0, 0, 0),
-                 (0, 3, 2, 0, 0),
-                 (0, 0, 3, 2, 0),
-                 (0, 0, 0, 3, 2),
-                 (2, 0, 0, 0, 3) ))
-
-    assert M.det(method="bareis") == 275
-    assert M.det(method="det_lu") == 275
-    assert M.det(method="Bareis") == 275
-
-    M = Matrix(( (1, 0,  1,  2, 12),
-                 (2, 0,  1,  1,  4),
-                 (2, 1,  1, -1,  3),
-                 (3, 2, -1,  1,  8),
-                 (1, 1,  1,  0,  6) ))
-
-    assert M.det(method="bareis") == -55
-    assert M.det(method="det_lu") == -55
-    assert M.det(method="BAREISS") == -55
-
-    M = Matrix(( ( 3,  0,  0, 0),
-                 (-2,  1,  0, 0),
-                 ( 0, -2,  5, 0),
-                 ( 5,  0,  3, 4) ))
-
-    assert M.det(method="bareiss") == 60
-    assert M.det(method="berkowitz") == 60
-    assert M.det(method="lu") == 60
-
-    M = Matrix(( ( 1,  0,  0,  0),
-                 ( 5,  0,  0,  0),
-                 ( 9, 10, 11, 0),
-                 (13, 14, 15, 16) ))
-
-    assert M.det(method="bareiss") == 0
-    assert M.det(method="berkowitz") == 0
-    assert M.det(method="lu") == 0
-
-    M = Matrix(( (3, 2, 0, 0, 0),
-                 (0, 3, 2, 0, 0),
-                 (0, 0, 3, 2, 0),
-                 (0, 0, 0, 3, 2),
-                 (0, 0, 0, 0, 3) ))
-
-    assert M.det(method="bareiss") == 243
-    assert M.det(method="berkowitz") == 243
-    assert M.det(method="lu") == 243
-
-    M = Matrix(( (-5,  2,  3,  4,  5),
-                 ( 1, -4,  3,  4,  5),
-                 ( 1,  2, -3,  4,  5),
-                 ( 1,  2,  3, -2,  5),
-                 ( 1,  2,  3,  4, -1) ))
-
-    assert M.det(method="bareis") == 11664
-    assert M.det(method="det_lu") == 11664
-    assert M.det(method="BERKOWITZ") == 11664
-
-    M = Matrix(( ( 2,  7, -1, 3, 2),
-                 ( 0,  0,  1, 0, 1),
-                 (-2,  0,  7, 0, 2),
-                 (-3, -2,  4, 5, 3),
-                 ( 1,  0,  0, 0, 1) ))
-
-    assert M.det(method="bareis") == 123
-    assert M.det(method="det_lu") == 123
-    assert M.det(method="LU") == 123
 
 def eye_Determinant(n):
     return Matrix(n, n, lambda i, j: int(i == j))
@@ -329,6 +191,8 @@ def test_adjugate():
     assert e.adjugate() == adj
     assert e.adjugate(method='bareiss') == adj
     assert e.adjugate(method='berkowitz') == adj
+    assert e.adjugate(method='bird') == adj
+    assert e.adjugate(method='laplace') == adj
 
     a = Matrix(2, 3, [1, 2, 3, 4, 5, 6])
     raises(NonSquareMatrixError, lambda: a.adjugate())
@@ -378,6 +242,8 @@ def test_cofactor_and_minors():
     assert e.cofactor_matrix() == cm
     assert e.cofactor_matrix(method="bareiss") == cm
     assert e.cofactor_matrix(method="berkowitz") == cm
+    assert e.cofactor_matrix(method="bird") == cm
+    assert e.cofactor_matrix(method="laplace") == cm
 
     raises(ValueError, lambda: e.cofactor(4, 5))
     raises(ValueError, lambda: e.minor(4, 5))
diff --git a/sympy/matrices/tests/test_solvers.py b/sympy/matrices/tests/test_solvers.py
index c555978c99..2b7c6375df 100644
--- a/sympy/matrices/tests/test_solvers.py
+++ b/sympy/matrices/tests/test_solvers.py
@@ -1,3 +1,4 @@
+import pytest
 from sympy.core.function import expand_mul
 from sympy.core.numbers import (I, Rational)
 from sympy.core.singleton import S
@@ -7,8 +8,10 @@
 from sympy.matrices.matrices import (ShapeError, NonSquareMatrixError)
 from sympy.matrices import (
     ImmutableMatrix, Matrix, eye, ones, ImmutableDenseMatrix, dotprodsimp)
+from sympy.matrices.determinant import _det_laplace
 from sympy.testing.pytest import raises
 from sympy.matrices.common import NonInvertibleMatrixError
+from sympy.polys.matrices.exceptions import DMShapeError
 from sympy.solvers.solveset import linsolve
 from sympy.abc import x, y
 
@@ -569,6 +572,27 @@ def test_linsolve_underdetermined_AND_gauss_jordan_solve():
     assert sol_2 == Matrix([[0], [0], [0]])
 
 
+@pytest.mark.parametrize("det_method", ["bird", "laplace"])
+@pytest.mark.parametrize("M, rhs", [
+    (Matrix([[2, 3, 5], [3, 6, 2], [8, 3, 6]]), Matrix(3, 1, [3, 7, 5])),
+    (Matrix([[2, 3, 5], [3, 6, 2], [8, 3, 6]]),
+     Matrix([[1, 2], [3, 4], [5, 6]])),
+    (Matrix(2, 2, symbols("a:4")), Matrix(2, 1, symbols("b:2"))),
+])
+def test_cramer_solve(det_method, M, rhs):
+    assert simplify(M.cramer_solve(rhs, det_method=det_method) - M.LUsolve(rhs)
+                    ) == Matrix.zeros(M.rows, rhs.cols)
+
+
+@pytest.mark.parametrize("det_method, error", [
+    ("bird", DMShapeError), (_det_laplace, NonSquareMatrixError)])
+def test_cramer_solve_errors(det_method, error):
+    # Non-square matrix
+    A = Matrix([[0, -1, 2], [5, 10, 7]])
+    b = Matrix([-2, 15])
+    raises(error, lambda: A.cramer_solve(b, det_method=det_method))
+
+
 def test_solve():
     A = Matrix([[1,2], [2,4]])
     b = Matrix([[3], [4]])
diff --git a/sympy/physics/mechanics/tests/test_kane3.py b/sympy/physics/mechanics/tests/test_kane3.py
index 431bd45d14..438759451c 100644
--- a/sympy/physics/mechanics/tests/test_kane3.py
+++ b/sympy/physics/mechanics/tests/test_kane3.py
@@ -1,4 +1,3 @@
-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
@@ -10,33 +9,6 @@
 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():
     # Code to get equations of motion for a bicycle modeled as in:
@@ -206,7 +178,7 @@ def test_bicycle():
             u_ind=[u2, u3, u5],
             u_dependent=[u1, u4, u6], velocity_constraints=conlist_speed,
             kd_eqs=kd,
-            constraint_solver=cramer_solve)
+            constraint_solver="CRAMER")
     (fr, frstar) = KM.kanes_equations(BL, FL)
 
     # This is the start of entering in the numerical values from the benchmark
@@ -311,7 +283,7 @@ def test_bicycle():
             q4: 0,
             q5: 0,
         },
-        linear_solver=cramer_solve,
+        linear_solver="CRAMER",
     )
     # 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