gh-39451: Fixes incorrect behaviour when taking the derivative of a constant matrix
    
<!-- ^ Please provide a concise and informative title. -->
<!-- ^ Don't put issue numbers in the title, do this in the PR
description below. -->
<!-- ^ For example, instead of "Fixes #12345" use "Introduce new method
to calculate 1 + 2". -->
<!-- v Describe your changes below in detail. -->
<!-- v Why is this change required? What problem does it solve? -->
<!-- v If this PR resolves an open issue, please link to it here. For
example, "Fixes #12345". -->

Fixes #15067.

Fixed error in the _derivative methods of matrix_dense.pyx and
matrix_sparse.pyx. These functions attempted to call the derivative
method belonging to each element of the matrix. This caused incorrect
behaviour when the element did not have a derivative method. This was
fixed by altering the code to call the derivative function defined in
/src/sage/calculus/functional.py with the element as a parameter,
instead of the derivative method belonging to the element. Doctests were
added to verify this change.

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. -->

- [X] The title is concise and informative.
- [X] The description explains in detail what this PR is about.
- [X] I have linked a relevant issue or discussion.
- [X] I have created tests covering the changes.
- [X] I have updated the documentation and checked the documentation
preview.

### ⌛ Dependencies

<!-- List all open PRs that this PR logically depends on. For example,
-->
<!-- - #12345: short description why this is a dependency -->
<!-- - #34567: ... -->
    
URL: #39451
Reported by: Caleb Van't Land
Reviewer(s):

Author: Release Manager

src/sage/matrix/matrix_dense.pyx

@@ -11,6 +11,7 @@ TESTS::
 cimport sage.matrix.matrix as matrix
 
 from sage.structure.richcmp cimport richcmp_item, rich_to_bool
+from sage.calculus.functional import derivative
 import sage.matrix.matrix_space
 import sage.structure.sequence
 
@@ -275,11 +276,19 @@ cdef class Matrix_dense(matrix.Matrix):
             sage: m._derivative(x)                                                      # needs sage.symbolic
             [    0     1]
             [  2*x 3*x^2]
+
+        TESTS:
+
+        Verify that :issue:`15067` is fixed::
+
+            sage: u = matrix(1, 2, [-1, 1])
+            sage: derivative(u, x)
+            [0 0]
         """
         # We could just use apply_map
         if self._nrows==0 or self._ncols==0:
             return self.__copy__()
-        v = [z.derivative(var) for z in self.list()]
+        v = [sage.calculus.functional.derivative(z, var) for z in self.list()]
         if R is None:
             v = sage.structure.sequence.Sequence(v)
             R = v.universe()

src/sage/matrix/matrix_sparse.pyx

@@ -16,6 +16,7 @@ from cysignals.signals cimport sig_check
 cimport sage.matrix.matrix as matrix
 cimport sage.matrix.matrix0 as matrix0
 from sage.categories.rings import Rings
+from sage.calculus.functional import derivative
 from sage.structure.element cimport Element, Vector
 from sage.structure.richcmp cimport richcmp_item, rich_to_bool
 
@@ -810,13 +811,23 @@ cdef class Matrix_sparse(matrix.Matrix):
             sage: m._derivative(x)                                                      # needs sage.symbolic
             [    0     1]
             [  2*x 3*x^2]
+
+        TESTS:
+
+        Verify that :issue:`15067` is fixed::
+
+            sage: m = matrix(3, 3, {(1, 1): 2, (0,2): 5})
+            sage: derivative(m, x)
+            [0 0 0]
+            [0 0 0]
+            [0 0 0]
         """
         # We would just use apply_map, except that Cython does not
         # allow lambda functions
 
         if self._nrows==0 or self._ncols==0:
             return self.__copy__()
-        v = [(ij, z.derivative(var)) for ij, z in self.dict().iteritems()]
+        v = [(ij, sage.calculus.functional.derivative(z, var)) for ij, z in self.dict().iteritems()]
         if R is None:
             w = [x for _, x in v]
             w = sage.structure.sequence.Sequence(w)