1091 lines
30 KiB
Python
1091 lines
30 KiB
Python
import random
|
|
|
|
from sympy.core.basic import Basic
|
|
from sympy.core.singleton import S
|
|
from sympy.core.symbol import Symbol
|
|
from sympy.core.sympify import sympify
|
|
from sympy.functions.elementary.trigonometric import cos, sin
|
|
from sympy.utilities.decorator import doctest_depends_on
|
|
from sympy.utilities.exceptions import sympy_deprecation_warning
|
|
from sympy.utilities.iterables import is_sequence
|
|
|
|
from .common import ShapeError
|
|
from .decompositions import _cholesky, _LDLdecomposition
|
|
from .matrices import MatrixBase
|
|
from .repmatrix import MutableRepMatrix, RepMatrix
|
|
from .solvers import _lower_triangular_solve, _upper_triangular_solve
|
|
|
|
|
|
def _iszero(x):
|
|
"""Returns True if x is zero."""
|
|
return x.is_zero
|
|
|
|
|
|
class DenseMatrix(RepMatrix):
|
|
"""Matrix implementation based on DomainMatrix as the internal representation"""
|
|
|
|
#
|
|
# DenseMatrix is a superclass for both MutableDenseMatrix and
|
|
# ImmutableDenseMatrix. Methods shared by both classes but not for the
|
|
# Sparse classes should be implemented here.
|
|
#
|
|
|
|
is_MatrixExpr = False # type: bool
|
|
|
|
_op_priority = 10.01
|
|
_class_priority = 4
|
|
|
|
@property
|
|
def _mat(self):
|
|
sympy_deprecation_warning(
|
|
"""
|
|
The private _mat attribute of Matrix is deprecated. Use the
|
|
.flat() method instead.
|
|
""",
|
|
deprecated_since_version="1.9",
|
|
active_deprecations_target="deprecated-private-matrix-attributes"
|
|
)
|
|
|
|
return self.flat()
|
|
|
|
def _eval_inverse(self, **kwargs):
|
|
return self.inv(method=kwargs.get('method', 'GE'),
|
|
iszerofunc=kwargs.get('iszerofunc', _iszero),
|
|
try_block_diag=kwargs.get('try_block_diag', False))
|
|
|
|
def as_immutable(self):
|
|
"""Returns an Immutable version of this Matrix
|
|
"""
|
|
from .immutable import ImmutableDenseMatrix as cls
|
|
return cls._fromrep(self._rep.copy())
|
|
|
|
def as_mutable(self):
|
|
"""Returns a mutable version of this matrix
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import ImmutableMatrix
|
|
>>> X = ImmutableMatrix([[1, 2], [3, 4]])
|
|
>>> Y = X.as_mutable()
|
|
>>> Y[1, 1] = 5 # Can set values in Y
|
|
>>> Y
|
|
Matrix([
|
|
[1, 2],
|
|
[3, 5]])
|
|
"""
|
|
return Matrix(self)
|
|
|
|
def cholesky(self, hermitian=True):
|
|
return _cholesky(self, hermitian=hermitian)
|
|
|
|
def LDLdecomposition(self, hermitian=True):
|
|
return _LDLdecomposition(self, hermitian=hermitian)
|
|
|
|
def lower_triangular_solve(self, rhs):
|
|
return _lower_triangular_solve(self, rhs)
|
|
|
|
def upper_triangular_solve(self, rhs):
|
|
return _upper_triangular_solve(self, rhs)
|
|
|
|
cholesky.__doc__ = _cholesky.__doc__
|
|
LDLdecomposition.__doc__ = _LDLdecomposition.__doc__
|
|
lower_triangular_solve.__doc__ = _lower_triangular_solve.__doc__
|
|
upper_triangular_solve.__doc__ = _upper_triangular_solve.__doc__
|
|
|
|
|
|
def _force_mutable(x):
|
|
"""Return a matrix as a Matrix, otherwise return x."""
|
|
if getattr(x, 'is_Matrix', False):
|
|
return x.as_mutable()
|
|
elif isinstance(x, Basic):
|
|
return x
|
|
elif hasattr(x, '__array__'):
|
|
a = x.__array__()
|
|
if len(a.shape) == 0:
|
|
return sympify(a)
|
|
return Matrix(x)
|
|
return x
|
|
|
|
|
|
class MutableDenseMatrix(DenseMatrix, MutableRepMatrix):
|
|
|
|
def simplify(self, **kwargs):
|
|
"""Applies simplify to the elements of a matrix in place.
|
|
|
|
This is a shortcut for M.applyfunc(lambda x: simplify(x, ratio, measure))
|
|
|
|
See Also
|
|
========
|
|
|
|
sympy.simplify.simplify.simplify
|
|
"""
|
|
from sympy.simplify.simplify import simplify as _simplify
|
|
for (i, j), element in self.todok().items():
|
|
self[i, j] = _simplify(element, **kwargs)
|
|
|
|
|
|
MutableMatrix = Matrix = MutableDenseMatrix
|
|
|
|
###########
|
|
# Numpy Utility Functions:
|
|
# list2numpy, matrix2numpy, symmarray
|
|
###########
|
|
|
|
|
|
def list2numpy(l, dtype=object): # pragma: no cover
|
|
"""Converts Python list of SymPy expressions to a NumPy array.
|
|
|
|
See Also
|
|
========
|
|
|
|
matrix2numpy
|
|
"""
|
|
from numpy import empty
|
|
a = empty(len(l), dtype)
|
|
for i, s in enumerate(l):
|
|
a[i] = s
|
|
return a
|
|
|
|
|
|
def matrix2numpy(m, dtype=object): # pragma: no cover
|
|
"""Converts SymPy's matrix to a NumPy array.
|
|
|
|
See Also
|
|
========
|
|
|
|
list2numpy
|
|
"""
|
|
from numpy import empty
|
|
a = empty(m.shape, dtype)
|
|
for i in range(m.rows):
|
|
for j in range(m.cols):
|
|
a[i, j] = m[i, j]
|
|
return a
|
|
|
|
|
|
###########
|
|
# Rotation matrices:
|
|
# rot_givens, rot_axis[123], rot_ccw_axis[123]
|
|
###########
|
|
|
|
|
|
def rot_givens(i, j, theta, dim=3):
|
|
r"""Returns a a Givens rotation matrix, a a rotation in the
|
|
plane spanned by two coordinates axes.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
The Givens rotation corresponds to a generalization of rotation
|
|
matrices to any number of dimensions, given by:
|
|
|
|
.. math::
|
|
G(i, j, \theta) =
|
|
\begin{bmatrix}
|
|
1 & \cdots & 0 & \cdots & 0 & \cdots & 0 \\
|
|
\vdots & \ddots & \vdots & & \vdots & & \vdots \\
|
|
0 & \cdots & c & \cdots & -s & \cdots & 0 \\
|
|
\vdots & & \vdots & \ddots & \vdots & & \vdots \\
|
|
0 & \cdots & s & \cdots & c & \cdots & 0 \\
|
|
\vdots & & \vdots & & \vdots & \ddots & \vdots \\
|
|
0 & \cdots & 0 & \cdots & 0 & \cdots & 1
|
|
\end{bmatrix}
|
|
|
|
Where $c = \cos(\theta)$ and $s = \sin(\theta)$ appear at the intersections
|
|
``i``\th and ``j``\th rows and columns.
|
|
|
|
For fixed ``i > j``\, the non-zero elements of a Givens matrix are
|
|
given by:
|
|
|
|
- $g_{kk} = 1$ for $k \ne i,\,j$
|
|
- $g_{kk} = c$ for $k = i,\,j$
|
|
- $g_{ji} = -g_{ij} = -s$
|
|
|
|
Parameters
|
|
==========
|
|
|
|
i : int between ``0`` and ``dim - 1``
|
|
Represents first axis
|
|
j : int between ``0`` and ``dim - 1``
|
|
Represents second axis
|
|
dim : int bigger than 1
|
|
Number of dimentions. Defaults to 3.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import pi, rot_givens
|
|
|
|
A counterclockwise rotation of pi/3 (60 degrees) around
|
|
the third axis (z-axis):
|
|
|
|
>>> rot_givens(1, 0, pi/3)
|
|
Matrix([
|
|
[ 1/2, -sqrt(3)/2, 0],
|
|
[sqrt(3)/2, 1/2, 0],
|
|
[ 0, 0, 1]])
|
|
|
|
If we rotate by pi/2 (90 degrees):
|
|
|
|
>>> rot_givens(1, 0, pi/2)
|
|
Matrix([
|
|
[0, -1, 0],
|
|
[1, 0, 0],
|
|
[0, 0, 1]])
|
|
|
|
This can be generalized to any number
|
|
of dimensions:
|
|
|
|
>>> rot_givens(1, 0, pi/2, dim=4)
|
|
Matrix([
|
|
[0, -1, 0, 0],
|
|
[1, 0, 0, 0],
|
|
[0, 0, 1, 0],
|
|
[0, 0, 0, 1]])
|
|
|
|
References
|
|
==========
|
|
|
|
.. [1] https://en.wikipedia.org/wiki/Givens_rotation
|
|
|
|
See Also
|
|
========
|
|
|
|
rot_axis1: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis (clockwise around the x axis)
|
|
rot_axis2: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis (clockwise around the y axis)
|
|
rot_axis3: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis (clockwise around the z axis)
|
|
rot_ccw_axis1: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis (counterclockwise around the x axis)
|
|
rot_ccw_axis2: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis (counterclockwise around the y axis)
|
|
rot_ccw_axis3: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis (counterclockwise around the z axis)
|
|
"""
|
|
if not isinstance(dim, int) or dim < 2:
|
|
raise ValueError('dim must be an integer biggen than one, '
|
|
'got {}.'.format(dim))
|
|
|
|
if i == j:
|
|
raise ValueError('i and j must be different, '
|
|
'got ({}, {})'.format(i, j))
|
|
|
|
for ij in [i, j]:
|
|
if not isinstance(ij, int) or ij < 0 or ij > dim - 1:
|
|
raise ValueError('i and j must be integers between 0 and '
|
|
'{}, got i={} and j={}.'.format(dim-1, i, j))
|
|
|
|
theta = sympify(theta)
|
|
c = cos(theta)
|
|
s = sin(theta)
|
|
M = eye(dim)
|
|
M[i, i] = c
|
|
M[j, j] = c
|
|
M[i, j] = s
|
|
M[j, i] = -s
|
|
return M
|
|
|
|
|
|
def rot_axis3(theta):
|
|
r"""Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
For a right-handed coordinate system, this corresponds to a
|
|
clockwise rotation around the `z`-axis, given by:
|
|
|
|
.. math::
|
|
|
|
R = \begin{bmatrix}
|
|
\cos(\theta) & \sin(\theta) & 0 \\
|
|
-\sin(\theta) & \cos(\theta) & 0 \\
|
|
0 & 0 & 1
|
|
\end{bmatrix}
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import pi, rot_axis3
|
|
|
|
A rotation of pi/3 (60 degrees):
|
|
|
|
>>> theta = pi/3
|
|
>>> rot_axis3(theta)
|
|
Matrix([
|
|
[ 1/2, sqrt(3)/2, 0],
|
|
[-sqrt(3)/2, 1/2, 0],
|
|
[ 0, 0, 1]])
|
|
|
|
If we rotate by pi/2 (90 degrees):
|
|
|
|
>>> rot_axis3(pi/2)
|
|
Matrix([
|
|
[ 0, 1, 0],
|
|
[-1, 0, 0],
|
|
[ 0, 0, 1]])
|
|
|
|
See Also
|
|
========
|
|
|
|
rot_givens: Returns a Givens rotation matrix (generalized rotation for
|
|
any number of dimensions)
|
|
rot_ccw_axis3: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis (counterclockwise around the z axis)
|
|
rot_axis1: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis (clockwise around the x axis)
|
|
rot_axis2: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis (clockwise around the y axis)
|
|
"""
|
|
return rot_givens(0, 1, theta, dim=3)
|
|
|
|
|
|
def rot_axis2(theta):
|
|
r"""Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
For a right-handed coordinate system, this corresponds to a
|
|
clockwise rotation around the `y`-axis, given by:
|
|
|
|
.. math::
|
|
|
|
R = \begin{bmatrix}
|
|
\cos(\theta) & 0 & -\sin(\theta) \\
|
|
0 & 1 & 0 \\
|
|
\sin(\theta) & 0 & \cos(\theta)
|
|
\end{bmatrix}
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import pi, rot_axis2
|
|
|
|
A rotation of pi/3 (60 degrees):
|
|
|
|
>>> theta = pi/3
|
|
>>> rot_axis2(theta)
|
|
Matrix([
|
|
[ 1/2, 0, -sqrt(3)/2],
|
|
[ 0, 1, 0],
|
|
[sqrt(3)/2, 0, 1/2]])
|
|
|
|
If we rotate by pi/2 (90 degrees):
|
|
|
|
>>> rot_axis2(pi/2)
|
|
Matrix([
|
|
[0, 0, -1],
|
|
[0, 1, 0],
|
|
[1, 0, 0]])
|
|
|
|
See Also
|
|
========
|
|
|
|
rot_givens: Returns a Givens rotation matrix (generalized rotation for
|
|
any number of dimensions)
|
|
rot_ccw_axis2: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis (clockwise around the y axis)
|
|
rot_axis1: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis (counterclockwise around the x axis)
|
|
rot_axis3: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis (counterclockwise around the z axis)
|
|
"""
|
|
return rot_givens(2, 0, theta, dim=3)
|
|
|
|
|
|
def rot_axis1(theta):
|
|
r"""Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
For a right-handed coordinate system, this corresponds to a
|
|
clockwise rotation around the `x`-axis, given by:
|
|
|
|
.. math::
|
|
|
|
R = \begin{bmatrix}
|
|
1 & 0 & 0 \\
|
|
0 & \cos(\theta) & \sin(\theta) \\
|
|
0 & -\sin(\theta) & \cos(\theta)
|
|
\end{bmatrix}
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import pi, rot_axis1
|
|
|
|
A rotation of pi/3 (60 degrees):
|
|
|
|
>>> theta = pi/3
|
|
>>> rot_axis1(theta)
|
|
Matrix([
|
|
[1, 0, 0],
|
|
[0, 1/2, sqrt(3)/2],
|
|
[0, -sqrt(3)/2, 1/2]])
|
|
|
|
If we rotate by pi/2 (90 degrees):
|
|
|
|
>>> rot_axis1(pi/2)
|
|
Matrix([
|
|
[1, 0, 0],
|
|
[0, 0, 1],
|
|
[0, -1, 0]])
|
|
|
|
See Also
|
|
========
|
|
|
|
rot_givens: Returns a Givens rotation matrix (generalized rotation for
|
|
any number of dimensions)
|
|
rot_ccw_axis1: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis (counterclockwise around the x axis)
|
|
rot_axis2: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis (clockwise around the y axis)
|
|
rot_axis3: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis (clockwise around the z axis)
|
|
"""
|
|
return rot_givens(1, 2, theta, dim=3)
|
|
|
|
|
|
def rot_ccw_axis3(theta):
|
|
r"""Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
For a right-handed coordinate system, this corresponds to a
|
|
counterclockwise rotation around the `z`-axis, given by:
|
|
|
|
.. math::
|
|
|
|
R = \begin{bmatrix}
|
|
\cos(\theta) & -\sin(\theta) & 0 \\
|
|
\sin(\theta) & \cos(\theta) & 0 \\
|
|
0 & 0 & 1
|
|
\end{bmatrix}
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import pi, rot_ccw_axis3
|
|
|
|
A rotation of pi/3 (60 degrees):
|
|
|
|
>>> theta = pi/3
|
|
>>> rot_ccw_axis3(theta)
|
|
Matrix([
|
|
[ 1/2, -sqrt(3)/2, 0],
|
|
[sqrt(3)/2, 1/2, 0],
|
|
[ 0, 0, 1]])
|
|
|
|
If we rotate by pi/2 (90 degrees):
|
|
|
|
>>> rot_ccw_axis3(pi/2)
|
|
Matrix([
|
|
[0, -1, 0],
|
|
[1, 0, 0],
|
|
[0, 0, 1]])
|
|
|
|
See Also
|
|
========
|
|
|
|
rot_givens: Returns a Givens rotation matrix (generalized rotation for
|
|
any number of dimensions)
|
|
rot_axis3: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis (clockwise around the z axis)
|
|
rot_ccw_axis1: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis (counterclockwise around the x axis)
|
|
rot_ccw_axis2: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis (counterclockwise around the y axis)
|
|
"""
|
|
return rot_givens(1, 0, theta, dim=3)
|
|
|
|
|
|
def rot_ccw_axis2(theta):
|
|
r"""Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
For a right-handed coordinate system, this corresponds to a
|
|
counterclockwise rotation around the `y`-axis, given by:
|
|
|
|
.. math::
|
|
|
|
R = \begin{bmatrix}
|
|
\cos(\theta) & 0 & \sin(\theta) \\
|
|
0 & 1 & 0 \\
|
|
-\sin(\theta) & 0 & \cos(\theta)
|
|
\end{bmatrix}
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import pi, rot_ccw_axis2
|
|
|
|
A rotation of pi/3 (60 degrees):
|
|
|
|
>>> theta = pi/3
|
|
>>> rot_ccw_axis2(theta)
|
|
Matrix([
|
|
[ 1/2, 0, sqrt(3)/2],
|
|
[ 0, 1, 0],
|
|
[-sqrt(3)/2, 0, 1/2]])
|
|
|
|
If we rotate by pi/2 (90 degrees):
|
|
|
|
>>> rot_ccw_axis2(pi/2)
|
|
Matrix([
|
|
[ 0, 0, 1],
|
|
[ 0, 1, 0],
|
|
[-1, 0, 0]])
|
|
|
|
See Also
|
|
========
|
|
|
|
rot_givens: Returns a Givens rotation matrix (generalized rotation for
|
|
any number of dimensions)
|
|
rot_axis2: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis (clockwise around the y axis)
|
|
rot_ccw_axis1: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis (counterclockwise around the x axis)
|
|
rot_ccw_axis3: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis (counterclockwise around the z axis)
|
|
"""
|
|
return rot_givens(0, 2, theta, dim=3)
|
|
|
|
|
|
def rot_ccw_axis1(theta):
|
|
r"""Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis.
|
|
|
|
Explanation
|
|
===========
|
|
|
|
For a right-handed coordinate system, this corresponds to a
|
|
counterclockwise rotation around the `x`-axis, given by:
|
|
|
|
.. math::
|
|
|
|
R = \begin{bmatrix}
|
|
1 & 0 & 0 \\
|
|
0 & \cos(\theta) & -\sin(\theta) \\
|
|
0 & \sin(\theta) & \cos(\theta)
|
|
\end{bmatrix}
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import pi, rot_ccw_axis1
|
|
|
|
A rotation of pi/3 (60 degrees):
|
|
|
|
>>> theta = pi/3
|
|
>>> rot_ccw_axis1(theta)
|
|
Matrix([
|
|
[1, 0, 0],
|
|
[0, 1/2, -sqrt(3)/2],
|
|
[0, sqrt(3)/2, 1/2]])
|
|
|
|
If we rotate by pi/2 (90 degrees):
|
|
|
|
>>> rot_ccw_axis1(pi/2)
|
|
Matrix([
|
|
[1, 0, 0],
|
|
[0, 0, -1],
|
|
[0, 1, 0]])
|
|
|
|
See Also
|
|
========
|
|
|
|
rot_givens: Returns a Givens rotation matrix (generalized rotation for
|
|
any number of dimensions)
|
|
rot_axis1: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 1-axis (clockwise around the x axis)
|
|
rot_ccw_axis2: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 2-axis (counterclockwise around the y axis)
|
|
rot_ccw_axis3: Returns a rotation matrix for a rotation of theta (in radians)
|
|
about the 3-axis (counterclockwise around the z axis)
|
|
"""
|
|
return rot_givens(2, 1, theta, dim=3)
|
|
|
|
|
|
@doctest_depends_on(modules=('numpy',))
|
|
def symarray(prefix, shape, **kwargs): # pragma: no cover
|
|
r"""Create a numpy ndarray of symbols (as an object array).
|
|
|
|
The created symbols are named ``prefix_i1_i2_``... You should thus provide a
|
|
non-empty prefix if you want your symbols to be unique for different output
|
|
arrays, as SymPy symbols with identical names are the same object.
|
|
|
|
Parameters
|
|
----------
|
|
|
|
prefix : string
|
|
A prefix prepended to the name of every symbol.
|
|
|
|
shape : int or tuple
|
|
Shape of the created array. If an int, the array is one-dimensional; for
|
|
more than one dimension the shape must be a tuple.
|
|
|
|
\*\*kwargs : dict
|
|
keyword arguments passed on to Symbol
|
|
|
|
Examples
|
|
========
|
|
These doctests require numpy.
|
|
|
|
>>> from sympy import symarray
|
|
>>> symarray('', 3)
|
|
[_0 _1 _2]
|
|
|
|
If you want multiple symarrays to contain distinct symbols, you *must*
|
|
provide unique prefixes:
|
|
|
|
>>> a = symarray('', 3)
|
|
>>> b = symarray('', 3)
|
|
>>> a[0] == b[0]
|
|
True
|
|
>>> a = symarray('a', 3)
|
|
>>> b = symarray('b', 3)
|
|
>>> a[0] == b[0]
|
|
False
|
|
|
|
Creating symarrays with a prefix:
|
|
|
|
>>> symarray('a', 3)
|
|
[a_0 a_1 a_2]
|
|
|
|
For more than one dimension, the shape must be given as a tuple:
|
|
|
|
>>> symarray('a', (2, 3))
|
|
[[a_0_0 a_0_1 a_0_2]
|
|
[a_1_0 a_1_1 a_1_2]]
|
|
>>> symarray('a', (2, 3, 2))
|
|
[[[a_0_0_0 a_0_0_1]
|
|
[a_0_1_0 a_0_1_1]
|
|
[a_0_2_0 a_0_2_1]]
|
|
<BLANKLINE>
|
|
[[a_1_0_0 a_1_0_1]
|
|
[a_1_1_0 a_1_1_1]
|
|
[a_1_2_0 a_1_2_1]]]
|
|
|
|
For setting assumptions of the underlying Symbols:
|
|
|
|
>>> [s.is_real for s in symarray('a', 2, real=True)]
|
|
[True, True]
|
|
"""
|
|
from numpy import empty, ndindex
|
|
arr = empty(shape, dtype=object)
|
|
for index in ndindex(shape):
|
|
arr[index] = Symbol('%s_%s' % (prefix, '_'.join(map(str, index))),
|
|
**kwargs)
|
|
return arr
|
|
|
|
|
|
###############
|
|
# Functions
|
|
###############
|
|
|
|
def casoratian(seqs, n, zero=True):
|
|
"""Given linear difference operator L of order 'k' and homogeneous
|
|
equation Ly = 0 we want to compute kernel of L, which is a set
|
|
of 'k' sequences: a(n), b(n), ... z(n).
|
|
|
|
Solutions of L are linearly independent iff their Casoratian,
|
|
denoted as C(a, b, ..., z), do not vanish for n = 0.
|
|
|
|
Casoratian is defined by k x k determinant::
|
|
|
|
+ a(n) b(n) . . . z(n) +
|
|
| a(n+1) b(n+1) . . . z(n+1) |
|
|
| . . . . |
|
|
| . . . . |
|
|
| . . . . |
|
|
+ a(n+k-1) b(n+k-1) . . . z(n+k-1) +
|
|
|
|
It proves very useful in rsolve_hyper() where it is applied
|
|
to a generating set of a recurrence to factor out linearly
|
|
dependent solutions and return a basis:
|
|
|
|
>>> from sympy import Symbol, casoratian, factorial
|
|
>>> n = Symbol('n', integer=True)
|
|
|
|
Exponential and factorial are linearly independent:
|
|
|
|
>>> casoratian([2**n, factorial(n)], n) != 0
|
|
True
|
|
|
|
"""
|
|
|
|
seqs = list(map(sympify, seqs))
|
|
|
|
if not zero:
|
|
f = lambda i, j: seqs[j].subs(n, n + i)
|
|
else:
|
|
f = lambda i, j: seqs[j].subs(n, i)
|
|
|
|
k = len(seqs)
|
|
|
|
return Matrix(k, k, f).det()
|
|
|
|
|
|
def eye(*args, **kwargs):
|
|
"""Create square identity matrix n x n
|
|
|
|
See Also
|
|
========
|
|
|
|
diag
|
|
zeros
|
|
ones
|
|
"""
|
|
|
|
return Matrix.eye(*args, **kwargs)
|
|
|
|
|
|
def diag(*values, strict=True, unpack=False, **kwargs):
|
|
"""Returns a matrix with the provided values placed on the
|
|
diagonal. If non-square matrices are included, they will
|
|
produce a block-diagonal matrix.
|
|
|
|
Examples
|
|
========
|
|
|
|
This version of diag is a thin wrapper to Matrix.diag that differs
|
|
in that it treats all lists like matrices -- even when a single list
|
|
is given. If this is not desired, either put a `*` before the list or
|
|
set `unpack=True`.
|
|
|
|
>>> from sympy import diag
|
|
|
|
>>> diag([1, 2, 3], unpack=True) # = diag(1,2,3) or diag(*[1,2,3])
|
|
Matrix([
|
|
[1, 0, 0],
|
|
[0, 2, 0],
|
|
[0, 0, 3]])
|
|
|
|
>>> diag([1, 2, 3]) # a column vector
|
|
Matrix([
|
|
[1],
|
|
[2],
|
|
[3]])
|
|
|
|
See Also
|
|
========
|
|
.common.MatrixCommon.eye
|
|
.common.MatrixCommon.diagonal
|
|
.common.MatrixCommon.diag
|
|
.expressions.blockmatrix.BlockMatrix
|
|
"""
|
|
return Matrix.diag(*values, strict=strict, unpack=unpack, **kwargs)
|
|
|
|
|
|
def GramSchmidt(vlist, orthonormal=False):
|
|
"""Apply the Gram-Schmidt process to a set of vectors.
|
|
|
|
Parameters
|
|
==========
|
|
|
|
vlist : List of Matrix
|
|
Vectors to be orthogonalized for.
|
|
|
|
orthonormal : Bool, optional
|
|
If true, return an orthonormal basis.
|
|
|
|
Returns
|
|
=======
|
|
|
|
vlist : List of Matrix
|
|
Orthogonalized vectors
|
|
|
|
Notes
|
|
=====
|
|
|
|
This routine is mostly duplicate from ``Matrix.orthogonalize``,
|
|
except for some difference that this always raises error when
|
|
linearly dependent vectors are found, and the keyword ``normalize``
|
|
has been named as ``orthonormal`` in this function.
|
|
|
|
See Also
|
|
========
|
|
|
|
.matrices.MatrixSubspaces.orthogonalize
|
|
|
|
References
|
|
==========
|
|
|
|
.. [1] https://en.wikipedia.org/wiki/Gram%E2%80%93Schmidt_process
|
|
"""
|
|
return MutableDenseMatrix.orthogonalize(
|
|
*vlist, normalize=orthonormal, rankcheck=True
|
|
)
|
|
|
|
|
|
def hessian(f, varlist, constraints=()):
|
|
"""Compute Hessian matrix for a function f wrt parameters in varlist
|
|
which may be given as a sequence or a row/column vector. A list of
|
|
constraints may optionally be given.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import Function, hessian, pprint
|
|
>>> from sympy.abc import x, y
|
|
>>> f = Function('f')(x, y)
|
|
>>> g1 = Function('g')(x, y)
|
|
>>> g2 = x**2 + 3*y
|
|
>>> pprint(hessian(f, (x, y), [g1, g2]))
|
|
[ d d ]
|
|
[ 0 0 --(g(x, y)) --(g(x, y)) ]
|
|
[ dx dy ]
|
|
[ ]
|
|
[ 0 0 2*x 3 ]
|
|
[ ]
|
|
[ 2 2 ]
|
|
[d d d ]
|
|
[--(g(x, y)) 2*x ---(f(x, y)) -----(f(x, y))]
|
|
[dx 2 dy dx ]
|
|
[ dx ]
|
|
[ ]
|
|
[ 2 2 ]
|
|
[d d d ]
|
|
[--(g(x, y)) 3 -----(f(x, y)) ---(f(x, y)) ]
|
|
[dy dy dx 2 ]
|
|
[ dy ]
|
|
|
|
References
|
|
==========
|
|
|
|
.. [1] https://en.wikipedia.org/wiki/Hessian_matrix
|
|
|
|
See Also
|
|
========
|
|
|
|
sympy.matrices.matrices.MatrixCalculus.jacobian
|
|
wronskian
|
|
"""
|
|
# f is the expression representing a function f, return regular matrix
|
|
if isinstance(varlist, MatrixBase):
|
|
if 1 not in varlist.shape:
|
|
raise ShapeError("`varlist` must be a column or row vector.")
|
|
if varlist.cols == 1:
|
|
varlist = varlist.T
|
|
varlist = varlist.tolist()[0]
|
|
if is_sequence(varlist):
|
|
n = len(varlist)
|
|
if not n:
|
|
raise ShapeError("`len(varlist)` must not be zero.")
|
|
else:
|
|
raise ValueError("Improper variable list in hessian function")
|
|
if not getattr(f, 'diff'):
|
|
# check differentiability
|
|
raise ValueError("Function `f` (%s) is not differentiable" % f)
|
|
m = len(constraints)
|
|
N = m + n
|
|
out = zeros(N)
|
|
for k, g in enumerate(constraints):
|
|
if not getattr(g, 'diff'):
|
|
# check differentiability
|
|
raise ValueError("Function `f` (%s) is not differentiable" % f)
|
|
for i in range(n):
|
|
out[k, i + m] = g.diff(varlist[i])
|
|
for i in range(n):
|
|
for j in range(i, n):
|
|
out[i + m, j + m] = f.diff(varlist[i]).diff(varlist[j])
|
|
for i in range(N):
|
|
for j in range(i + 1, N):
|
|
out[j, i] = out[i, j]
|
|
return out
|
|
|
|
|
|
def jordan_cell(eigenval, n):
|
|
"""
|
|
Create a Jordan block:
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import jordan_cell
|
|
>>> from sympy.abc import x
|
|
>>> jordan_cell(x, 4)
|
|
Matrix([
|
|
[x, 1, 0, 0],
|
|
[0, x, 1, 0],
|
|
[0, 0, x, 1],
|
|
[0, 0, 0, x]])
|
|
"""
|
|
|
|
return Matrix.jordan_block(size=n, eigenvalue=eigenval)
|
|
|
|
|
|
def matrix_multiply_elementwise(A, B):
|
|
"""Return the Hadamard product (elementwise product) of A and B
|
|
|
|
>>> from sympy import Matrix, matrix_multiply_elementwise
|
|
>>> A = Matrix([[0, 1, 2], [3, 4, 5]])
|
|
>>> B = Matrix([[1, 10, 100], [100, 10, 1]])
|
|
>>> matrix_multiply_elementwise(A, B)
|
|
Matrix([
|
|
[ 0, 10, 200],
|
|
[300, 40, 5]])
|
|
|
|
See Also
|
|
========
|
|
|
|
sympy.matrices.common.MatrixCommon.__mul__
|
|
"""
|
|
return A.multiply_elementwise(B)
|
|
|
|
|
|
def ones(*args, **kwargs):
|
|
"""Returns a matrix of ones with ``rows`` rows and ``cols`` columns;
|
|
if ``cols`` is omitted a square matrix will be returned.
|
|
|
|
See Also
|
|
========
|
|
|
|
zeros
|
|
eye
|
|
diag
|
|
"""
|
|
|
|
if 'c' in kwargs:
|
|
kwargs['cols'] = kwargs.pop('c')
|
|
|
|
return Matrix.ones(*args, **kwargs)
|
|
|
|
|
|
def randMatrix(r, c=None, min=0, max=99, seed=None, symmetric=False,
|
|
percent=100, prng=None):
|
|
"""Create random matrix with dimensions ``r`` x ``c``. If ``c`` is omitted
|
|
the matrix will be square. If ``symmetric`` is True the matrix must be
|
|
square. If ``percent`` is less than 100 then only approximately the given
|
|
percentage of elements will be non-zero.
|
|
|
|
The pseudo-random number generator used to generate matrix is chosen in the
|
|
following way.
|
|
|
|
* If ``prng`` is supplied, it will be used as random number generator.
|
|
It should be an instance of ``random.Random``, or at least have
|
|
``randint`` and ``shuffle`` methods with same signatures.
|
|
* if ``prng`` is not supplied but ``seed`` is supplied, then new
|
|
``random.Random`` with given ``seed`` will be created;
|
|
* otherwise, a new ``random.Random`` with default seed will be used.
|
|
|
|
Examples
|
|
========
|
|
|
|
>>> from sympy import randMatrix
|
|
>>> randMatrix(3) # doctest:+SKIP
|
|
[25, 45, 27]
|
|
[44, 54, 9]
|
|
[23, 96, 46]
|
|
>>> randMatrix(3, 2) # doctest:+SKIP
|
|
[87, 29]
|
|
[23, 37]
|
|
[90, 26]
|
|
>>> randMatrix(3, 3, 0, 2) # doctest:+SKIP
|
|
[0, 2, 0]
|
|
[2, 0, 1]
|
|
[0, 0, 1]
|
|
>>> randMatrix(3, symmetric=True) # doctest:+SKIP
|
|
[85, 26, 29]
|
|
[26, 71, 43]
|
|
[29, 43, 57]
|
|
>>> A = randMatrix(3, seed=1)
|
|
>>> B = randMatrix(3, seed=2)
|
|
>>> A == B
|
|
False
|
|
>>> A == randMatrix(3, seed=1)
|
|
True
|
|
>>> randMatrix(3, symmetric=True, percent=50) # doctest:+SKIP
|
|
[77, 70, 0],
|
|
[70, 0, 0],
|
|
[ 0, 0, 88]
|
|
"""
|
|
# Note that ``Random()`` is equivalent to ``Random(None)``
|
|
prng = prng or random.Random(seed)
|
|
|
|
if c is None:
|
|
c = r
|
|
|
|
if symmetric and r != c:
|
|
raise ValueError('For symmetric matrices, r must equal c, but %i != %i' % (r, c))
|
|
|
|
ij = range(r * c)
|
|
if percent != 100:
|
|
ij = prng.sample(ij, int(len(ij)*percent // 100))
|
|
|
|
m = zeros(r, c)
|
|
|
|
if not symmetric:
|
|
for ijk in ij:
|
|
i, j = divmod(ijk, c)
|
|
m[i, j] = prng.randint(min, max)
|
|
else:
|
|
for ijk in ij:
|
|
i, j = divmod(ijk, c)
|
|
if i <= j:
|
|
m[i, j] = m[j, i] = prng.randint(min, max)
|
|
|
|
return m
|
|
|
|
|
|
def wronskian(functions, var, method='bareiss'):
|
|
"""
|
|
Compute Wronskian for [] of functions
|
|
|
|
::
|
|
|
|
| f1 f2 ... fn |
|
|
| f1' f2' ... fn' |
|
|
| . . . . |
|
|
W(f1, ..., fn) = | . . . . |
|
|
| . . . . |
|
|
| (n) (n) (n) |
|
|
| D (f1) D (f2) ... D (fn) |
|
|
|
|
see: https://en.wikipedia.org/wiki/Wronskian
|
|
|
|
See Also
|
|
========
|
|
|
|
sympy.matrices.matrices.MatrixCalculus.jacobian
|
|
hessian
|
|
"""
|
|
|
|
functions = [sympify(f) for f in functions]
|
|
n = len(functions)
|
|
if n == 0:
|
|
return S.One
|
|
W = Matrix(n, n, lambda i, j: functions[i].diff(var, j))
|
|
return W.det(method)
|
|
|
|
|
|
def zeros(*args, **kwargs):
|
|
"""Returns a matrix of zeros with ``rows`` rows and ``cols`` columns;
|
|
if ``cols`` is omitted a square matrix will be returned.
|
|
|
|
See Also
|
|
========
|
|
|
|
ones
|
|
eye
|
|
diag
|
|
"""
|
|
|
|
if 'c' in kwargs:
|
|
kwargs['cols'] = kwargs.pop('c')
|
|
|
|
return Matrix.zeros(*args, **kwargs)
|