Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

Pull 140-v2: New matrices support in sympy.physics.matrices #515

Closed
wants to merge 3 commits into from

6 participants

@plaes

This is effectively pull request 140 which is rebased and original commits are squashed.

Original is here: #140

@smichr
Collaborator

The old request can be closed now. (It will still be accessible through the URL you cited above.)

@lazovich

Test results html report: http://pastehtml.com/view/b1vl2g758.html

Automatic review by sympy-bot.

sympy/matrices/tests/test_matrices.py
@@ -1168,11 +1164,9 @@ def test_nonvectorJacobian():
x, y, z = symbols('x y z')
X = Matrix([ [exp(x + y + z), exp(x + y + z)],
[exp(x + y + z), exp(x + y + z)] ])
- Y = Matrix([x, y, z])
- raises(TypeError, 'X.jacobian(Y)')
+ raises(TypeError, 'X.jacobian(Matrix([x, y, z])')

Looks like there is a missing paren inside the quotes here. I think it should be

raises(TypeError, 'X.jacobian(Matrix[x y, z]))')

@plaes
plaes added a note

Uh-oh.. I pushed this patch from wrong machine (I had this issue fixed on my laptop, but I ran tests on my desktop machine, and this is where I fixed it) :S

Fixed now..

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@lazovich

Test results html report: http://pastehtml.com/view/b1wjwwir0.html

Automatic review by sympy-bot.

@lazovich

@plaes: All tests pass now, thanks!

This looks good to me, and I think it was pretty thoroughly reviewed in the last PR, so I am +1 on getting it merged.

@asmeurer
Owner

Test results html report: http://pastehtml.com/view/b3ik6bni5.html

Summary: There do not appear to be any problems.

Please double check this against the test report given above before merging with master.

Automatic review by sympy-bot.

@plaes

I rebased and fixed new doctest issues 3**(1/2) vs sqrt(3).

@certik
Owner

I have rebased, fixed a few test failures and pushed it in. Thanks!

@certik certik closed this
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
This page is out of date. Refresh to see the latest.
View
2  sympy/matrices/__init__.py
@@ -5,6 +5,6 @@
from matrices import (Matrix, SparseMatrix, zeros, ones, eye, diag,
hessian, randMatrix, GramSchmidt, wronskian, casoratian,
list2numpy, matrix2numpy, DeferredVector, block_diag, symarray, ShapeError,
- NonSquareMatrixError)
+ NonSquareMatrixError, rot_axis1, rot_axis2, rot_axis3)
from expressions import *
View
106 sympy/matrices/matrices.py
@@ -1,14 +1,12 @@
from sympy import Basic, Symbol, Integer, C, S, Dummy, Rational, Add, Pow
-from sympy.core.numbers import Zero
from sympy.core.sympify import sympify, converter, SympifyError
from sympy.core.compatibility import is_sequence
from sympy.polys import Poly, roots, cancel
from sympy.simplify import simplify as sympy_simplify
-from sympy.utilities.iterables import flatten
from sympy.functions.elementary.miscellaneous import sqrt, Max, Min
-from sympy.functions.elementary.complexes import re, Abs
from sympy.printing import sstr
+from sympy.functions.elementary.trigonometric import cos, sin
from sympy.core.compatibility import callable, reduce
@@ -578,7 +576,7 @@ def LDLdecomposition(self):
"""
if not self.is_square:
- raise NonSquareMatrixException("Matrix must be square.")
+ raise NonSquareMatrixError("Matrix must be square.")
if not self.is_symmetric():
raise ValueError("Matrix must be symmetric.")
return self._LDLdecomposition()
@@ -606,7 +604,7 @@ def lower_triangular_solve(self, rhs):
"""
if not self.is_square:
- raise NonSquareMatrixException("Matrix must be square.")
+ raise NonSquareMatrixError("Matrix must be square.")
if rhs.rows != self.rows:
raise ShapeError("Matrices size mismatch.")
if not self.is_lower():
@@ -633,7 +631,7 @@ def upper_triangular_solve(self, rhs):
"""
if not self.is_square:
- raise NonSquareMatrixException("Matrix must be square.")
+ raise NonSquareMatrixError("Matrix must be square.")
if rhs.rows != self.rows:
raise TypeError("Matrix size mismatch.")
if not self.is_upper():
@@ -3277,9 +3275,93 @@ def symarray(prefix, shape):
arr[index] = Symbol('%s_%s' % (prefix, '_'.join(map(str, index))))
return arr
-def _separate_eig_results(res):
- eigvals = [item[0] for item in res]
- multiplicities = [item[1] for item in res]
- eigvals = flatten([[val]*mult for val, mult in zip(eigVals, multiplicities)])
- eigvects = flatten([item[2] for item in res])
- return eigvals, eigvects
+def rot_axis3(theta):
+ """Returns a rotation matrix for a rotation of theta (in radians) about
+ the 3-axis.
+
+ Examples
+ --------
+
+ >>> from sympy import pi
+ >>> from sympy.matrices import rot_axis3
+
+ A rotation of pi/3 (60 degrees):
+ >>> theta = pi/3
+ >>> rot_axis3(theta)
+ [ 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)
+ [ 0, 1, 0]
+ [-1, 0, 0]
+ [ 0, 0, 1]
+ """
+ ct = cos(theta)
+ st = sin(theta)
+ mat = ((ct,st,0),
+ (-st,ct,0),
+ (0,0,1))
+ return Matrix(mat)
+
+def rot_axis2(theta):
+ """Returns a rotation matrix for a rotation of theta (in radians) about
+ the 2-axis.
+
+ Examples
+ --------
+
+ >>> from sympy import pi
+ >>> from sympy.matrices import rot_axis2
+
+ A rotation of pi/3 (60 degrees):
+ >>> theta = pi/3
+ >>> rot_axis2(theta)
+ [ 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)
+ [0, 0, -1]
+ [0, 1, 0]
+ [1, 0, 0]
+ """
+ ct = cos(theta)
+ st = sin(theta)
+ mat = ((ct,0,-st),
+ (0,1,0),
+ (st,0,ct))
+ return Matrix(mat)
+
+def rot_axis1(theta):
+ """Returns a rotation matrix for a rotation of theta (in radians) about
+ the 1-axis.
+
+ Examples
+ --------
+
+ >>> from sympy import pi
+ >>> from sympy.matrices import rot_axis1
+
+ A rotation of pi/3 (60 degrees):
+ >>> theta = pi/3
+ >>> rot_axis1(theta)
+ [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)
+ [1, 0, 0]
+ [0, 0, 1]
+ [0, -1, 0]
+ """
+ ct = cos(theta)
+ st = sin(theta)
+ mat = ((1,0,0),
+ (0,ct,st),
+ (0,-st,ct))
+ return Matrix(mat)
+
View
52 sympy/matrices/tests/test_matrices.py
@@ -1,15 +1,11 @@
from sympy import (symbols, Matrix, SparseMatrix, eye, I, Symbol, Rational,
Float, wronskian, cos, sin, exp, hessian, sqrt, zeros, ones, randMatrix,
- Poly, S, pi, E, I, oo, trigsimp, Integer, block_diag, N, zeros, sympify,
+ Poly, S, pi, E, oo, trigsimp, Integer, N, sympify,
Pow, simplify, Min, Max, Abs)
from sympy.matrices.matrices import (ShapeError, MatrixError,
- matrix_multiply_elementwise, diag,
-
- SparseMatrix, SparseMatrix, NonSquareMatrixError, _dims_to_nm,
- matrix_multiply_elementwise)
+ matrix_multiply_elementwise, diag, NonSquareMatrixError, _dims_to_nm)
from sympy.utilities.pytest import raises
-#from sympy.functions.elementary.miscellaneous import Max, Min
-#from sympy.functions.elementary.miscellaneous import Max, Min
+from sympy.matrices import rot_axis1, rot_axis2, rot_axis3
def test_division():
x, y, z = symbols('x y z')
@@ -1168,11 +1164,9 @@ def test_nonvectorJacobian():
x, y, z = symbols('x y z')
X = Matrix([ [exp(x + y + z), exp(x + y + z)],
[exp(x + y + z), exp(x + y + z)] ])
- Y = Matrix([x, y, z])
- raises(TypeError, 'X.jacobian(Y)')
+ raises(TypeError, 'X.jacobian(Matrix([x, y, z]))')
X = X[0,:]
- Y = Matrix([ [x, y], [x,z] ])
- raises(TypeError, 'X.jacobian(Y)')
+ raises(TypeError, 'X.jacobian(Matrix([ [x, y], [x, z] ]))')
def test_vec():
m = Matrix([ [1,3], [2,4] ])
@@ -1199,10 +1193,8 @@ def test_vech():
assert m_vech[0] == y*x
def test_vech_errors():
- m = Matrix([ [1,3] ])
- raises(ShapeError, 'm.vech()')
- m = Matrix([ [1,3], [2,4] ])
- raises(ValueError, 'm.vech()')
+ raises(ShapeError, 'Matrix([ [1,3] ]).vech()')
+ raises(ValueError, 'Matrix([ [1,3], [2,4] ]).vech()')
def test_diag():
x, y, z = symbols("x y z")
@@ -1754,14 +1746,6 @@ def test_condition_number():
assert all(Float(1.).epsilon_eq(Mc.subs(x, val).evalf()) for val in \
[Rational(1,5), Rational(1, 2), Rational(1, 10), pi/2, pi, 7*pi/4 ])
-def test_len():
- assert len(Matrix()) == 0
- assert len(Matrix([[1, 2]])) == len(Matrix([[1], [2]])) == 2
- assert len(Matrix(0, 2, lambda i, j: 0)) == len(Matrix(2, 0, lambda i, j: 0)) == 0
- assert len(Matrix([[0, 1, 2], [3, 4, 5]])) == 6
- assert Matrix([1])
- assert not Matrix()
-
def test_equality():
A = Matrix(((1,2,3),(4,5,6),(7,8,9)))
B = Matrix(((9,8,7),(6,5,4),(3,2,1)))
@@ -1778,3 +1762,25 @@ def test_equality():
assert C == D
assert not C != D
+def test_rotation_matrices():
+ # This tests the rotation matrices by rotating about an axis and back.
+ theta = pi/3
+ r3_plus = rot_axis3(theta)
+ r3_minus = rot_axis3(-theta)
+ r2_plus = rot_axis2(theta)
+ r2_minus = rot_axis2(-theta)
+ r1_plus = rot_axis1(theta)
+ r1_minus = rot_axis1(-theta)
+ assert r3_minus*r3_plus*eye(3) == eye(3)
+ assert r2_minus*r2_plus*eye(3) == eye(3)
+ assert r1_minus*r1_plus*eye(3) == eye(3)
+
+ # Check the correctness of the trace of the rotation matrix
+ assert r1_plus.trace() == 1 + 2*cos(theta)
+ assert r2_plus.trace() == 1 + 2*cos(theta)
+ assert r3_plus.trace() == 1 + 2*cos(theta)
+
+ # Check that a rotation with zero angle doesn't change anything.
+ assert rot_axis1(0) == eye(3)
+ assert rot_axis2(0) == eye(3)
+ assert rot_axis3(0) == eye(3)
View
27 sympy/physics/matrices.py
@@ -26,6 +26,33 @@ def msigma(i):
raise IndexError("Invalid Pauli index")
return Matrix(mat)
+def pat_matrix(m, dx, dy, dz):
+ """Returns the Parallel Axis Theorem matrix to translate the inertia
+ matrix a distance of (dx, dy, dz) for a body of mass m.
+
+ Examples
+ --------
+ If the point we want the inertia about is a distance of 2 units of
+ length and 1 unit along the x-axis we get:
+ >>> from sympy.physics.matrices import pat_matrix
+ >>> pat_matrix(2,1,0,0)
+ [0, 0, 0]
+ [0, 2, 0]
+ [0, 0, 2]
+
+ In case we want to find the inertia along a vector of (1,1,1):
+ >>> pat_matrix(2,1,1,1)
+ [ 4, -2, -2]
+ [-2, 4, -2]
+ [-2, -2, 4]
+ """
+ dxdy = -dx*dy ; dydz = -dy*dz ; dzdx = -dz*dx
+ dxdx = dx**2 ; dydy = dy**2 ; dzdz = dz**2
+ mat = ((dydy + dzdz, dxdy, dzdx),
+ (dxdy, dxdx + dzdz, dydz),
+ (dzdx, dydz, dydy + dxdx))
+ return m*Matrix(mat)
+
def mgamma(mu,lower=False):
"""Returns a Dirac gamma matrix gamma^mu in the standard
(Dirac) representation.
View
26 sympy/physics/tests/test_physics_matrices.py
@@ -1,7 +1,29 @@
-from sympy.physics.matrices import msigma, mgamma, minkowski_tensor
-from sympy import zeros, eye, I
+from sympy.physics.matrices import msigma, mgamma, minkowski_tensor, pat_matrix
+from sympy import zeros, eye, I, Matrix
+def test_parallel_axis_theorem():
+ # This tests the parallel axis theorem matrix by comparing to test
+ # matrices.
+ # First case, 1 in all directions.
+ mat1 = Matrix(((2,-1,-1),(-1,2,-1),(-1,-1,2)))
+ assert pat_matrix(1,1,1,1) == mat1
+ assert pat_matrix(2,1,1,1) == 2*mat1
+
+ # Second case, 1 in x, 0 in all others
+ mat2 = Matrix(((0,0,0),(0,1,0),(0,0,1)))
+ assert pat_matrix(1,1,0,0) == mat2
+ assert pat_matrix(2,1,0,0) == 2*mat2
+
+ # Third case, 1 in y, 0 in all others
+ mat3 = Matrix(((1,0,0),(0,0,0),(0,0,1)))
+ assert pat_matrix(1,0,1,0) == mat3
+ assert pat_matrix(2,0,1,0) == 2*mat3
+
+ # Fourth case, 1 in z, 0 in all others
+ mat4 = Matrix(((1,0,0),(0,1,0),(0,0,0)))
+ assert pat_matrix(1,0,0,1) == mat4
+ assert pat_matrix(2,0,0,1) == 2*mat4
def test_Pauli():
#this and the following test are testing both Pauli and Dirac matrices
Something went wrong with that request. Please try again.