Coverage for /usr/lib/python3/dist-packages/sympy/matrices/dense.py: 32%
169 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1import random
3from sympy.core.basic import Basic
4from sympy.core.singleton import S
5from sympy.core.symbol import Symbol
6from sympy.core.sympify import sympify
7from sympy.functions.elementary.trigonometric import cos, sin
8from sympy.utilities.decorator import doctest_depends_on
9from sympy.utilities.exceptions import sympy_deprecation_warning
10from sympy.utilities.iterables import is_sequence
12from .common import ShapeError
13from .decompositions import _cholesky, _LDLdecomposition
14from .matrices import MatrixBase
15from .repmatrix import MutableRepMatrix, RepMatrix
16from .solvers import _lower_triangular_solve, _upper_triangular_solve
19def _iszero(x):
20 """Returns True if x is zero."""
21 return x.is_zero
24class DenseMatrix(RepMatrix):
25 """Matrix implementation based on DomainMatrix as the internal representation"""
27 #
28 # DenseMatrix is a superclass for both MutableDenseMatrix and
29 # ImmutableDenseMatrix. Methods shared by both classes but not for the
30 # Sparse classes should be implemented here.
31 #
33 is_MatrixExpr = False # type: bool
35 _op_priority = 10.01
36 _class_priority = 4
38 @property
39 def _mat(self):
40 sympy_deprecation_warning(
41 """
42 The private _mat attribute of Matrix is deprecated. Use the
43 .flat() method instead.
44 """,
45 deprecated_since_version="1.9",
46 active_deprecations_target="deprecated-private-matrix-attributes"
47 )
49 return self.flat()
51 def _eval_inverse(self, **kwargs):
52 return self.inv(method=kwargs.get('method', 'GE'),
53 iszerofunc=kwargs.get('iszerofunc', _iszero),
54 try_block_diag=kwargs.get('try_block_diag', False))
56 def as_immutable(self):
57 """Returns an Immutable version of this Matrix
58 """
59 from .immutable import ImmutableDenseMatrix as cls
60 return cls._fromrep(self._rep.copy())
62 def as_mutable(self):
63 """Returns a mutable version of this matrix
65 Examples
66 ========
68 >>> from sympy import ImmutableMatrix
69 >>> X = ImmutableMatrix([[1, 2], [3, 4]])
70 >>> Y = X.as_mutable()
71 >>> Y[1, 1] = 5 # Can set values in Y
72 >>> Y
73 Matrix([
74 [1, 2],
75 [3, 5]])
76 """
77 return Matrix(self)
79 def cholesky(self, hermitian=True):
80 return _cholesky(self, hermitian=hermitian)
82 def LDLdecomposition(self, hermitian=True):
83 return _LDLdecomposition(self, hermitian=hermitian)
85 def lower_triangular_solve(self, rhs):
86 return _lower_triangular_solve(self, rhs)
88 def upper_triangular_solve(self, rhs):
89 return _upper_triangular_solve(self, rhs)
91 cholesky.__doc__ = _cholesky.__doc__
92 LDLdecomposition.__doc__ = _LDLdecomposition.__doc__
93 lower_triangular_solve.__doc__ = _lower_triangular_solve.__doc__
94 upper_triangular_solve.__doc__ = _upper_triangular_solve.__doc__
97def _force_mutable(x):
98 """Return a matrix as a Matrix, otherwise return x."""
99 if getattr(x, 'is_Matrix', False):
100 return x.as_mutable()
101 elif isinstance(x, Basic):
102 return x
103 elif hasattr(x, '__array__'):
104 a = x.__array__()
105 if len(a.shape) == 0:
106 return sympify(a)
107 return Matrix(x)
108 return x
111class MutableDenseMatrix(DenseMatrix, MutableRepMatrix):
113 def simplify(self, **kwargs):
114 """Applies simplify to the elements of a matrix in place.
116 This is a shortcut for M.applyfunc(lambda x: simplify(x, ratio, measure))
118 See Also
119 ========
121 sympy.simplify.simplify.simplify
122 """
123 from sympy.simplify.simplify import simplify as _simplify
124 for (i, j), element in self.todok().items():
125 self[i, j] = _simplify(element, **kwargs)
128MutableMatrix = Matrix = MutableDenseMatrix
130###########
131# Numpy Utility Functions:
132# list2numpy, matrix2numpy, symmarray
133###########
136def list2numpy(l, dtype=object): # pragma: no cover
137 """Converts Python list of SymPy expressions to a NumPy array.
139 See Also
140 ========
142 matrix2numpy
143 """
144 from numpy import empty
145 a = empty(len(l), dtype)
146 for i, s in enumerate(l):
147 a[i] = s
148 return a
151def matrix2numpy(m, dtype=object): # pragma: no cover
152 """Converts SymPy's matrix to a NumPy array.
154 See Also
155 ========
157 list2numpy
158 """
159 from numpy import empty
160 a = empty(m.shape, dtype)
161 for i in range(m.rows):
162 for j in range(m.cols):
163 a[i, j] = m[i, j]
164 return a
167###########
168# Rotation matrices:
169# rot_givens, rot_axis[123], rot_ccw_axis[123]
170###########
173def rot_givens(i, j, theta, dim=3):
174 r"""Returns a a Givens rotation matrix, a a rotation in the
175 plane spanned by two coordinates axes.
177 Explanation
178 ===========
180 The Givens rotation corresponds to a generalization of rotation
181 matrices to any number of dimensions, given by:
183 .. math::
184 G(i, j, \theta) =
185 \begin{bmatrix}
186 1 & \cdots & 0 & \cdots & 0 & \cdots & 0 \\
187 \vdots & \ddots & \vdots & & \vdots & & \vdots \\
188 0 & \cdots & c & \cdots & -s & \cdots & 0 \\
189 \vdots & & \vdots & \ddots & \vdots & & \vdots \\
190 0 & \cdots & s & \cdots & c & \cdots & 0 \\
191 \vdots & & \vdots & & \vdots & \ddots & \vdots \\
192 0 & \cdots & 0 & \cdots & 0 & \cdots & 1
193 \end{bmatrix}
195 Where $c = \cos(\theta)$ and $s = \sin(\theta)$ appear at the intersections
196 ``i``\th and ``j``\th rows and columns.
198 For fixed ``i > j``\, the non-zero elements of a Givens matrix are
199 given by:
201 - $g_{kk} = 1$ for $k \ne i,\,j$
202 - $g_{kk} = c$ for $k = i,\,j$
203 - $g_{ji} = -g_{ij} = -s$
205 Parameters
206 ==========
208 i : int between ``0`` and ``dim - 1``
209 Represents first axis
210 j : int between ``0`` and ``dim - 1``
211 Represents second axis
212 dim : int bigger than 1
213 Number of dimentions. Defaults to 3.
215 Examples
216 ========
218 >>> from sympy import pi, rot_givens
220 A counterclockwise rotation of pi/3 (60 degrees) around
221 the third axis (z-axis):
223 >>> rot_givens(1, 0, pi/3)
224 Matrix([
225 [ 1/2, -sqrt(3)/2, 0],
226 [sqrt(3)/2, 1/2, 0],
227 [ 0, 0, 1]])
229 If we rotate by pi/2 (90 degrees):
231 >>> rot_givens(1, 0, pi/2)
232 Matrix([
233 [0, -1, 0],
234 [1, 0, 0],
235 [0, 0, 1]])
237 This can be generalized to any number
238 of dimensions:
240 >>> rot_givens(1, 0, pi/2, dim=4)
241 Matrix([
242 [0, -1, 0, 0],
243 [1, 0, 0, 0],
244 [0, 0, 1, 0],
245 [0, 0, 0, 1]])
247 References
248 ==========
250 .. [1] https://en.wikipedia.org/wiki/Givens_rotation
252 See Also
253 ========
255 rot_axis1: Returns a rotation matrix for a rotation of theta (in radians)
256 about the 1-axis (clockwise around the x axis)
257 rot_axis2: Returns a rotation matrix for a rotation of theta (in radians)
258 about the 2-axis (clockwise around the y axis)
259 rot_axis3: Returns a rotation matrix for a rotation of theta (in radians)
260 about the 3-axis (clockwise around the z axis)
261 rot_ccw_axis1: Returns a rotation matrix for a rotation of theta (in radians)
262 about the 1-axis (counterclockwise around the x axis)
263 rot_ccw_axis2: Returns a rotation matrix for a rotation of theta (in radians)
264 about the 2-axis (counterclockwise around the y axis)
265 rot_ccw_axis3: Returns a rotation matrix for a rotation of theta (in radians)
266 about the 3-axis (counterclockwise around the z axis)
267 """
268 if not isinstance(dim, int) or dim < 2:
269 raise ValueError('dim must be an integer biggen than one, '
270 'got {}.'.format(dim))
272 if i == j:
273 raise ValueError('i and j must be different, '
274 'got ({}, {})'.format(i, j))
276 for ij in [i, j]:
277 if not isinstance(ij, int) or ij < 0 or ij > dim - 1:
278 raise ValueError('i and j must be integers between 0 and '
279 '{}, got i={} and j={}.'.format(dim-1, i, j))
281 theta = sympify(theta)
282 c = cos(theta)
283 s = sin(theta)
284 M = eye(dim)
285 M[i, i] = c
286 M[j, j] = c
287 M[i, j] = s
288 M[j, i] = -s
289 return M
292def rot_axis3(theta):
293 r"""Returns a rotation matrix for a rotation of theta (in radians)
294 about the 3-axis.
296 Explanation
297 ===========
299 For a right-handed coordinate system, this corresponds to a
300 clockwise rotation around the `z`-axis, given by:
302 .. math::
304 R = \begin{bmatrix}
305 \cos(\theta) & \sin(\theta) & 0 \\
306 -\sin(\theta) & \cos(\theta) & 0 \\
307 0 & 0 & 1
308 \end{bmatrix}
310 Examples
311 ========
313 >>> from sympy import pi, rot_axis3
315 A rotation of pi/3 (60 degrees):
317 >>> theta = pi/3
318 >>> rot_axis3(theta)
319 Matrix([
320 [ 1/2, sqrt(3)/2, 0],
321 [-sqrt(3)/2, 1/2, 0],
322 [ 0, 0, 1]])
324 If we rotate by pi/2 (90 degrees):
326 >>> rot_axis3(pi/2)
327 Matrix([
328 [ 0, 1, 0],
329 [-1, 0, 0],
330 [ 0, 0, 1]])
332 See Also
333 ========
335 rot_givens: Returns a Givens rotation matrix (generalized rotation for
336 any number of dimensions)
337 rot_ccw_axis3: Returns a rotation matrix for a rotation of theta (in radians)
338 about the 3-axis (counterclockwise around the z axis)
339 rot_axis1: Returns a rotation matrix for a rotation of theta (in radians)
340 about the 1-axis (clockwise around the x axis)
341 rot_axis2: Returns a rotation matrix for a rotation of theta (in radians)
342 about the 2-axis (clockwise around the y axis)
343 """
344 return rot_givens(0, 1, theta, dim=3)
347def rot_axis2(theta):
348 r"""Returns a rotation matrix for a rotation of theta (in radians)
349 about the 2-axis.
351 Explanation
352 ===========
354 For a right-handed coordinate system, this corresponds to a
355 clockwise rotation around the `y`-axis, given by:
357 .. math::
359 R = \begin{bmatrix}
360 \cos(\theta) & 0 & -\sin(\theta) \\
361 0 & 1 & 0 \\
362 \sin(\theta) & 0 & \cos(\theta)
363 \end{bmatrix}
365 Examples
366 ========
368 >>> from sympy import pi, rot_axis2
370 A rotation of pi/3 (60 degrees):
372 >>> theta = pi/3
373 >>> rot_axis2(theta)
374 Matrix([
375 [ 1/2, 0, -sqrt(3)/2],
376 [ 0, 1, 0],
377 [sqrt(3)/2, 0, 1/2]])
379 If we rotate by pi/2 (90 degrees):
381 >>> rot_axis2(pi/2)
382 Matrix([
383 [0, 0, -1],
384 [0, 1, 0],
385 [1, 0, 0]])
387 See Also
388 ========
390 rot_givens: Returns a Givens rotation matrix (generalized rotation for
391 any number of dimensions)
392 rot_ccw_axis2: Returns a rotation matrix for a rotation of theta (in radians)
393 about the 2-axis (clockwise around the y axis)
394 rot_axis1: Returns a rotation matrix for a rotation of theta (in radians)
395 about the 1-axis (counterclockwise around the x axis)
396 rot_axis3: Returns a rotation matrix for a rotation of theta (in radians)
397 about the 3-axis (counterclockwise around the z axis)
398 """
399 return rot_givens(2, 0, theta, dim=3)
402def rot_axis1(theta):
403 r"""Returns a rotation matrix for a rotation of theta (in radians)
404 about the 1-axis.
406 Explanation
407 ===========
409 For a right-handed coordinate system, this corresponds to a
410 clockwise rotation around the `x`-axis, given by:
412 .. math::
414 R = \begin{bmatrix}
415 1 & 0 & 0 \\
416 0 & \cos(\theta) & \sin(\theta) \\
417 0 & -\sin(\theta) & \cos(\theta)
418 \end{bmatrix}
420 Examples
421 ========
423 >>> from sympy import pi, rot_axis1
425 A rotation of pi/3 (60 degrees):
427 >>> theta = pi/3
428 >>> rot_axis1(theta)
429 Matrix([
430 [1, 0, 0],
431 [0, 1/2, sqrt(3)/2],
432 [0, -sqrt(3)/2, 1/2]])
434 If we rotate by pi/2 (90 degrees):
436 >>> rot_axis1(pi/2)
437 Matrix([
438 [1, 0, 0],
439 [0, 0, 1],
440 [0, -1, 0]])
442 See Also
443 ========
445 rot_givens: Returns a Givens rotation matrix (generalized rotation for
446 any number of dimensions)
447 rot_ccw_axis1: Returns a rotation matrix for a rotation of theta (in radians)
448 about the 1-axis (counterclockwise around the x axis)
449 rot_axis2: Returns a rotation matrix for a rotation of theta (in radians)
450 about the 2-axis (clockwise around the y axis)
451 rot_axis3: Returns a rotation matrix for a rotation of theta (in radians)
452 about the 3-axis (clockwise around the z axis)
453 """
454 return rot_givens(1, 2, theta, dim=3)
457def rot_ccw_axis3(theta):
458 r"""Returns a rotation matrix for a rotation of theta (in radians)
459 about the 3-axis.
461 Explanation
462 ===========
464 For a right-handed coordinate system, this corresponds to a
465 counterclockwise rotation around the `z`-axis, given by:
467 .. math::
469 R = \begin{bmatrix}
470 \cos(\theta) & -\sin(\theta) & 0 \\
471 \sin(\theta) & \cos(\theta) & 0 \\
472 0 & 0 & 1
473 \end{bmatrix}
475 Examples
476 ========
478 >>> from sympy import pi, rot_ccw_axis3
480 A rotation of pi/3 (60 degrees):
482 >>> theta = pi/3
483 >>> rot_ccw_axis3(theta)
484 Matrix([
485 [ 1/2, -sqrt(3)/2, 0],
486 [sqrt(3)/2, 1/2, 0],
487 [ 0, 0, 1]])
489 If we rotate by pi/2 (90 degrees):
491 >>> rot_ccw_axis3(pi/2)
492 Matrix([
493 [0, -1, 0],
494 [1, 0, 0],
495 [0, 0, 1]])
497 See Also
498 ========
500 rot_givens: Returns a Givens rotation matrix (generalized rotation for
501 any number of dimensions)
502 rot_axis3: Returns a rotation matrix for a rotation of theta (in radians)
503 about the 3-axis (clockwise around the z axis)
504 rot_ccw_axis1: Returns a rotation matrix for a rotation of theta (in radians)
505 about the 1-axis (counterclockwise around the x axis)
506 rot_ccw_axis2: Returns a rotation matrix for a rotation of theta (in radians)
507 about the 2-axis (counterclockwise around the y axis)
508 """
509 return rot_givens(1, 0, theta, dim=3)
512def rot_ccw_axis2(theta):
513 r"""Returns a rotation matrix for a rotation of theta (in radians)
514 about the 2-axis.
516 Explanation
517 ===========
519 For a right-handed coordinate system, this corresponds to a
520 counterclockwise rotation around the `y`-axis, given by:
522 .. math::
524 R = \begin{bmatrix}
525 \cos(\theta) & 0 & \sin(\theta) \\
526 0 & 1 & 0 \\
527 -\sin(\theta) & 0 & \cos(\theta)
528 \end{bmatrix}
530 Examples
531 ========
533 >>> from sympy import pi, rot_ccw_axis2
535 A rotation of pi/3 (60 degrees):
537 >>> theta = pi/3
538 >>> rot_ccw_axis2(theta)
539 Matrix([
540 [ 1/2, 0, sqrt(3)/2],
541 [ 0, 1, 0],
542 [-sqrt(3)/2, 0, 1/2]])
544 If we rotate by pi/2 (90 degrees):
546 >>> rot_ccw_axis2(pi/2)
547 Matrix([
548 [ 0, 0, 1],
549 [ 0, 1, 0],
550 [-1, 0, 0]])
552 See Also
553 ========
555 rot_givens: Returns a Givens rotation matrix (generalized rotation for
556 any number of dimensions)
557 rot_axis2: Returns a rotation matrix for a rotation of theta (in radians)
558 about the 2-axis (clockwise around the y axis)
559 rot_ccw_axis1: Returns a rotation matrix for a rotation of theta (in radians)
560 about the 1-axis (counterclockwise around the x axis)
561 rot_ccw_axis3: Returns a rotation matrix for a rotation of theta (in radians)
562 about the 3-axis (counterclockwise around the z axis)
563 """
564 return rot_givens(0, 2, theta, dim=3)
567def rot_ccw_axis1(theta):
568 r"""Returns a rotation matrix for a rotation of theta (in radians)
569 about the 1-axis.
571 Explanation
572 ===========
574 For a right-handed coordinate system, this corresponds to a
575 counterclockwise rotation around the `x`-axis, given by:
577 .. math::
579 R = \begin{bmatrix}
580 1 & 0 & 0 \\
581 0 & \cos(\theta) & -\sin(\theta) \\
582 0 & \sin(\theta) & \cos(\theta)
583 \end{bmatrix}
585 Examples
586 ========
588 >>> from sympy import pi, rot_ccw_axis1
590 A rotation of pi/3 (60 degrees):
592 >>> theta = pi/3
593 >>> rot_ccw_axis1(theta)
594 Matrix([
595 [1, 0, 0],
596 [0, 1/2, -sqrt(3)/2],
597 [0, sqrt(3)/2, 1/2]])
599 If we rotate by pi/2 (90 degrees):
601 >>> rot_ccw_axis1(pi/2)
602 Matrix([
603 [1, 0, 0],
604 [0, 0, -1],
605 [0, 1, 0]])
607 See Also
608 ========
610 rot_givens: Returns a Givens rotation matrix (generalized rotation for
611 any number of dimensions)
612 rot_axis1: Returns a rotation matrix for a rotation of theta (in radians)
613 about the 1-axis (clockwise around the x axis)
614 rot_ccw_axis2: Returns a rotation matrix for a rotation of theta (in radians)
615 about the 2-axis (counterclockwise around the y axis)
616 rot_ccw_axis3: Returns a rotation matrix for a rotation of theta (in radians)
617 about the 3-axis (counterclockwise around the z axis)
618 """
619 return rot_givens(2, 1, theta, dim=3)
622@doctest_depends_on(modules=('numpy',))
623def symarray(prefix, shape, **kwargs): # pragma: no cover
624 r"""Create a numpy ndarray of symbols (as an object array).
626 The created symbols are named ``prefix_i1_i2_``... You should thus provide a
627 non-empty prefix if you want your symbols to be unique for different output
628 arrays, as SymPy symbols with identical names are the same object.
630 Parameters
631 ----------
633 prefix : string
634 A prefix prepended to the name of every symbol.
636 shape : int or tuple
637 Shape of the created array. If an int, the array is one-dimensional; for
638 more than one dimension the shape must be a tuple.
640 \*\*kwargs : dict
641 keyword arguments passed on to Symbol
643 Examples
644 ========
645 These doctests require numpy.
647 >>> from sympy import symarray
648 >>> symarray('', 3)
649 [_0 _1 _2]
651 If you want multiple symarrays to contain distinct symbols, you *must*
652 provide unique prefixes:
654 >>> a = symarray('', 3)
655 >>> b = symarray('', 3)
656 >>> a[0] == b[0]
657 True
658 >>> a = symarray('a', 3)
659 >>> b = symarray('b', 3)
660 >>> a[0] == b[0]
661 False
663 Creating symarrays with a prefix:
665 >>> symarray('a', 3)
666 [a_0 a_1 a_2]
668 For more than one dimension, the shape must be given as a tuple:
670 >>> symarray('a', (2, 3))
671 [[a_0_0 a_0_1 a_0_2]
672 [a_1_0 a_1_1 a_1_2]]
673 >>> symarray('a', (2, 3, 2))
674 [[[a_0_0_0 a_0_0_1]
675 [a_0_1_0 a_0_1_1]
676 [a_0_2_0 a_0_2_1]]
677 <BLANKLINE>
678 [[a_1_0_0 a_1_0_1]
679 [a_1_1_0 a_1_1_1]
680 [a_1_2_0 a_1_2_1]]]
682 For setting assumptions of the underlying Symbols:
684 >>> [s.is_real for s in symarray('a', 2, real=True)]
685 [True, True]
686 """
687 from numpy import empty, ndindex
688 arr = empty(shape, dtype=object)
689 for index in ndindex(shape):
690 arr[index] = Symbol('%s_%s' % (prefix, '_'.join(map(str, index))),
691 **kwargs)
692 return arr
695###############
696# Functions
697###############
699def casoratian(seqs, n, zero=True):
700 """Given linear difference operator L of order 'k' and homogeneous
701 equation Ly = 0 we want to compute kernel of L, which is a set
702 of 'k' sequences: a(n), b(n), ... z(n).
704 Solutions of L are linearly independent iff their Casoratian,
705 denoted as C(a, b, ..., z), do not vanish for n = 0.
707 Casoratian is defined by k x k determinant::
709 + a(n) b(n) . . . z(n) +
710 | a(n+1) b(n+1) . . . z(n+1) |
711 | . . . . |
712 | . . . . |
713 | . . . . |
714 + a(n+k-1) b(n+k-1) . . . z(n+k-1) +
716 It proves very useful in rsolve_hyper() where it is applied
717 to a generating set of a recurrence to factor out linearly
718 dependent solutions and return a basis:
720 >>> from sympy import Symbol, casoratian, factorial
721 >>> n = Symbol('n', integer=True)
723 Exponential and factorial are linearly independent:
725 >>> casoratian([2**n, factorial(n)], n) != 0
726 True
728 """
730 seqs = list(map(sympify, seqs))
732 if not zero:
733 f = lambda i, j: seqs[j].subs(n, n + i)
734 else:
735 f = lambda i, j: seqs[j].subs(n, i)
737 k = len(seqs)
739 return Matrix(k, k, f).det()
742def eye(*args, **kwargs):
743 """Create square identity matrix n x n
745 See Also
746 ========
748 diag
749 zeros
750 ones
751 """
753 return Matrix.eye(*args, **kwargs)
756def diag(*values, strict=True, unpack=False, **kwargs):
757 """Returns a matrix with the provided values placed on the
758 diagonal. If non-square matrices are included, they will
759 produce a block-diagonal matrix.
761 Examples
762 ========
764 This version of diag is a thin wrapper to Matrix.diag that differs
765 in that it treats all lists like matrices -- even when a single list
766 is given. If this is not desired, either put a `*` before the list or
767 set `unpack=True`.
769 >>> from sympy import diag
771 >>> diag([1, 2, 3], unpack=True) # = diag(1,2,3) or diag(*[1,2,3])
772 Matrix([
773 [1, 0, 0],
774 [0, 2, 0],
775 [0, 0, 3]])
777 >>> diag([1, 2, 3]) # a column vector
778 Matrix([
779 [1],
780 [2],
781 [3]])
783 See Also
784 ========
785 .common.MatrixCommon.eye
786 .common.MatrixCommon.diagonal
787 .common.MatrixCommon.diag
788 .expressions.blockmatrix.BlockMatrix
789 """
790 return Matrix.diag(*values, strict=strict, unpack=unpack, **kwargs)
793def GramSchmidt(vlist, orthonormal=False):
794 """Apply the Gram-Schmidt process to a set of vectors.
796 Parameters
797 ==========
799 vlist : List of Matrix
800 Vectors to be orthogonalized for.
802 orthonormal : Bool, optional
803 If true, return an orthonormal basis.
805 Returns
806 =======
808 vlist : List of Matrix
809 Orthogonalized vectors
811 Notes
812 =====
814 This routine is mostly duplicate from ``Matrix.orthogonalize``,
815 except for some difference that this always raises error when
816 linearly dependent vectors are found, and the keyword ``normalize``
817 has been named as ``orthonormal`` in this function.
819 See Also
820 ========
822 .matrices.MatrixSubspaces.orthogonalize
824 References
825 ==========
827 .. [1] https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process
828 """
829 return MutableDenseMatrix.orthogonalize(
830 *vlist, normalize=orthonormal, rankcheck=True
831 )
834def hessian(f, varlist, constraints=()):
835 """Compute Hessian matrix for a function f wrt parameters in varlist
836 which may be given as a sequence or a row/column vector. A list of
837 constraints may optionally be given.
839 Examples
840 ========
842 >>> from sympy import Function, hessian, pprint
843 >>> from sympy.abc import x, y
844 >>> f = Function('f')(x, y)
845 >>> g1 = Function('g')(x, y)
846 >>> g2 = x**2 + 3*y
847 >>> pprint(hessian(f, (x, y), [g1, g2]))
848 [ d d ]
849 [ 0 0 --(g(x, y)) --(g(x, y)) ]
850 [ dx dy ]
851 [ ]
852 [ 0 0 2*x 3 ]
853 [ ]
854 [ 2 2 ]
855 [d d d ]
856 [--(g(x, y)) 2*x ---(f(x, y)) -----(f(x, y))]
857 [dx 2 dy dx ]
858 [ dx ]
859 [ ]
860 [ 2 2 ]
861 [d d d ]
862 [--(g(x, y)) 3 -----(f(x, y)) ---(f(x, y)) ]
863 [dy dy dx 2 ]
864 [ dy ]
866 References
867 ==========
869 .. [1] https://en.wikipedia.org/wiki/Hessian_matrix
871 See Also
872 ========
874 sympy.matrices.matrices.MatrixCalculus.jacobian
875 wronskian
876 """
877 # f is the expression representing a function f, return regular matrix
878 if isinstance(varlist, MatrixBase):
879 if 1 not in varlist.shape:
880 raise ShapeError("`varlist` must be a column or row vector.")
881 if varlist.cols == 1:
882 varlist = varlist.T
883 varlist = varlist.tolist()[0]
884 if is_sequence(varlist):
885 n = len(varlist)
886 if not n:
887 raise ShapeError("`len(varlist)` must not be zero.")
888 else:
889 raise ValueError("Improper variable list in hessian function")
890 if not getattr(f, 'diff'):
891 # check differentiability
892 raise ValueError("Function `f` (%s) is not differentiable" % f)
893 m = len(constraints)
894 N = m + n
895 out = zeros(N)
896 for k, g in enumerate(constraints):
897 if not getattr(g, 'diff'):
898 # check differentiability
899 raise ValueError("Function `f` (%s) is not differentiable" % f)
900 for i in range(n):
901 out[k, i + m] = g.diff(varlist[i])
902 for i in range(n):
903 for j in range(i, n):
904 out[i + m, j + m] = f.diff(varlist[i]).diff(varlist[j])
905 for i in range(N):
906 for j in range(i + 1, N):
907 out[j, i] = out[i, j]
908 return out
911def jordan_cell(eigenval, n):
912 """
913 Create a Jordan block:
915 Examples
916 ========
918 >>> from sympy import jordan_cell
919 >>> from sympy.abc import x
920 >>> jordan_cell(x, 4)
921 Matrix([
922 [x, 1, 0, 0],
923 [0, x, 1, 0],
924 [0, 0, x, 1],
925 [0, 0, 0, x]])
926 """
928 return Matrix.jordan_block(size=n, eigenvalue=eigenval)
931def matrix_multiply_elementwise(A, B):
932 """Return the Hadamard product (elementwise product) of A and B
934 >>> from sympy import Matrix, matrix_multiply_elementwise
935 >>> A = Matrix([[0, 1, 2], [3, 4, 5]])
936 >>> B = Matrix([[1, 10, 100], [100, 10, 1]])
937 >>> matrix_multiply_elementwise(A, B)
938 Matrix([
939 [ 0, 10, 200],
940 [300, 40, 5]])
942 See Also
943 ========
945 sympy.matrices.common.MatrixCommon.__mul__
946 """
947 return A.multiply_elementwise(B)
950def ones(*args, **kwargs):
951 """Returns a matrix of ones with ``rows`` rows and ``cols`` columns;
952 if ``cols`` is omitted a square matrix will be returned.
954 See Also
955 ========
957 zeros
958 eye
959 diag
960 """
962 if 'c' in kwargs:
963 kwargs['cols'] = kwargs.pop('c')
965 return Matrix.ones(*args, **kwargs)
968def randMatrix(r, c=None, min=0, max=99, seed=None, symmetric=False,
969 percent=100, prng=None):
970 """Create random matrix with dimensions ``r`` x ``c``. If ``c`` is omitted
971 the matrix will be square. If ``symmetric`` is True the matrix must be
972 square. If ``percent`` is less than 100 then only approximately the given
973 percentage of elements will be non-zero.
975 The pseudo-random number generator used to generate matrix is chosen in the
976 following way.
978 * If ``prng`` is supplied, it will be used as random number generator.
979 It should be an instance of ``random.Random``, or at least have
980 ``randint`` and ``shuffle`` methods with same signatures.
981 * if ``prng`` is not supplied but ``seed`` is supplied, then new
982 ``random.Random`` with given ``seed`` will be created;
983 * otherwise, a new ``random.Random`` with default seed will be used.
985 Examples
986 ========
988 >>> from sympy import randMatrix
989 >>> randMatrix(3) # doctest:+SKIP
990 [25, 45, 27]
991 [44, 54, 9]
992 [23, 96, 46]
993 >>> randMatrix(3, 2) # doctest:+SKIP
994 [87, 29]
995 [23, 37]
996 [90, 26]
997 >>> randMatrix(3, 3, 0, 2) # doctest:+SKIP
998 [0, 2, 0]
999 [2, 0, 1]
1000 [0, 0, 1]
1001 >>> randMatrix(3, symmetric=True) # doctest:+SKIP
1002 [85, 26, 29]
1003 [26, 71, 43]
1004 [29, 43, 57]
1005 >>> A = randMatrix(3, seed=1)
1006 >>> B = randMatrix(3, seed=2)
1007 >>> A == B
1008 False
1009 >>> A == randMatrix(3, seed=1)
1010 True
1011 >>> randMatrix(3, symmetric=True, percent=50) # doctest:+SKIP
1012 [77, 70, 0],
1013 [70, 0, 0],
1014 [ 0, 0, 88]
1015 """
1016 # Note that ``Random()`` is equivalent to ``Random(None)``
1017 prng = prng or random.Random(seed)
1019 if c is None:
1020 c = r
1022 if symmetric and r != c:
1023 raise ValueError('For symmetric matrices, r must equal c, but %i != %i' % (r, c))
1025 ij = range(r * c)
1026 if percent != 100:
1027 ij = prng.sample(ij, int(len(ij)*percent // 100))
1029 m = zeros(r, c)
1031 if not symmetric:
1032 for ijk in ij:
1033 i, j = divmod(ijk, c)
1034 m[i, j] = prng.randint(min, max)
1035 else:
1036 for ijk in ij:
1037 i, j = divmod(ijk, c)
1038 if i <= j:
1039 m[i, j] = m[j, i] = prng.randint(min, max)
1041 return m
1044def wronskian(functions, var, method='bareiss'):
1045 """
1046 Compute Wronskian for [] of functions
1048 ::
1050 | f1 f2 ... fn |
1051 | f1' f2' ... fn' |
1052 | . . . . |
1053 W(f1, ..., fn) = | . . . . |
1054 | . . . . |
1055 | (n) (n) (n) |
1056 | D (f1) D (f2) ... D (fn) |
1058 see: https://en.wikipedia.org/wiki/Wronskian
1060 See Also
1061 ========
1063 sympy.matrices.matrices.MatrixCalculus.jacobian
1064 hessian
1065 """
1067 functions = [sympify(f) for f in functions]
1068 n = len(functions)
1069 if n == 0:
1070 return S.One
1071 W = Matrix(n, n, lambda i, j: functions[i].diff(var, j))
1072 return W.det(method)
1075def zeros(*args, **kwargs):
1076 """Returns a matrix of zeros with ``rows`` rows and ``cols`` columns;
1077 if ``cols`` is omitted a square matrix will be returned.
1079 See Also
1080 ========
1082 ones
1083 eye
1084 diag
1085 """
1087 if 'c' in kwargs:
1088 kwargs['cols'] = kwargs.pop('c')
1090 return Matrix.zeros(*args, **kwargs)