details: http://hg.sympy.org/sympy/rev/5fc85c34d2ae
changeset: 1857:5fc85c34d2ae
user: Henrik Johansson <[EMAIL PROTECTED]>
date: Thu Nov 06 22:23:24 2008 +0100
description:
matrices: optional zero finder function argument in the inv method implemented.
plus a note about the pivoting problem in the docstring.
Signed-off-by: Henrik Johansson <[EMAIL PROTECTED]>
Signed-off-by: Ondrej Certik <[EMAIL PROTECTED]>
diffs (163 lines):
diff -r 4d0cd200f57b -r 5fc85c34d2ae sympy/matrices/matrices.py
--- a/sympy/matrices/matrices.py Wed Nov 05 13:18:07 2008 +0100
+++ b/sympy/matrices/matrices.py Thu Nov 06 22:23:24 2008 +0100
@@ -45,6 +45,9 @@
raise ValueError("Matrix dimensions should be a two-element tuple of
ints or a single int!")
return n, m
+
+def _iszero(x):
+ return x == 0
class Matrix(object):
@@ -469,7 +472,7 @@
def __repr__(self):
return StrPrinter.doprint(self)
- def inv(self, method="GE"):
+ def inv(self, method="GE", iszerofunc=_iszero):
"""
Calculates the matrix inverse.
@@ -479,12 +482,18 @@
LU .... inverse_LU()
ADJ ... inverse_ADJ()
+ Note, the GE and LU methods may require the matrix to be simplified
+ before it is inverted in order to properly detect zeros during
+ pivoting. In difficult cases a custom zero detection function can
+ be provided by setting the iszerosfunc argument to a function that
+ should return True if its argument is zero.
+
"""
assert self.cols==self.lines
if method == "GE":
- return self.inverse_GE()
+ return self.inverse_GE(iszerofunc=iszerofunc)
elif method == "LU":
- return self.inverse_LU()
+ return self.inverse_LU(iszerofunc=iszerofunc)
elif method == "ADJ":
return self.inverse_ADJ()
else:
@@ -756,13 +765,13 @@
s+="]\n"
print s
- def LUsolve(self, rhs):
+ def LUsolve(self, rhs, iszerofunc=_iszero):
"""
Solve the linear system Ax = b.
self is the coefficient matrix A and rhs is the right side b.
"""
assert rhs.lines == self.lines
- A, perm = self.LUdecomposition_Simple()
+ A, perm = self.LUdecomposition_Simple(iszerofunc=_iszero)
n = self.lines
b = rhs.permuteFwd(perm)
# forward substitution, all diag entries are scaled to 1
@@ -776,11 +785,11 @@
b.row(i, lambda x,k: x / A[i,i])
return b
- def LUdecomposition(self):
+ def LUdecomposition(self, iszerofunc=_iszero):
"""
Returns the decompositon LU and the row swaps p.
"""
- combined, p = self.LUdecomposition_Simple()
+ combined, p = self.LUdecomposition_Simple(iszerofunc=_iszero)
L = self.zeros(self.lines)
U = self.zeros(self.lines)
for i in range(self.lines):
@@ -793,7 +802,7 @@
U[i,j] = combined[i,j]
return L, U, p
- def LUdecomposition_Simple(self):
+ def LUdecomposition_Simple(self, iszerofunc=_iszero):
"""
Returns A compused of L,U (L's diag entries are 1) and
p which is the list of the row swaps (in order).
@@ -812,14 +821,14 @@
for k in range(j):
A[i,j] = A[i,j] - A[i,k]*A[k,j]
# find the first non-zero pivot, includes any expression
- if pivot == -1 and A[i,j] != 0:
+ if pivot == -1 and not iszerofunc(A[i,j]):
pivot = i
if pivot < 0:
raise "Error: non-invertible matrix passed to
LUdecomposition_Simple()"
if pivot != j: # row must be swapped
A.row_swap(pivot,j)
p.append([pivot,j])
- assert A[j,j] != 0
+ assert not iszerofunc(A[j,j])
scale = 1 / A[j,j]
for i in range(j+1,n):
A[i,j] = A[i,j] * scale
@@ -1145,20 +1154,20 @@
return self.cofactorMatrix(method).T
- def inverse_LU(self):
+ def inverse_LU(self, iszerofunc=_iszero):
"""
Calculates the inverse using LU decomposition.
"""
- return self.LUsolve(self.eye(self.lines))
+ return self.LUsolve(self.eye(self.lines), iszerofunc=_iszero)
- def inverse_GE(self):
+ def inverse_GE(self, iszerofunc=_iszero):
"""
Calculates the inverse using Gaussian elimination.
"""
assert self.lines == self.cols
assert self.det() != 0
big = self.row_join(self.eye(self.lines))
- red = big.rref()
+ red = big.rref(iszerofunc=iszerofunc)
return red[0][:,big.lines:]
def inverse_ADJ(self):
@@ -1170,7 +1179,7 @@
assert d != 0
return self.adjugate()/d
- def rref(self,simplified=False):
+ def rref(self,simplified=False, iszerofunc=_iszero):
"""
Take any matrix and return reduced row-echelon form and indices of
pivot vars
@@ -1184,13 +1193,13 @@
break
if simplified:
r[pivots,i] = simplify(r[pivots,i])
- if r[pivots,i] == 0:
+ if iszerofunc(r[pivots,i]):
for k in range(pivots, r.lines):
if simplified and k>pivots:
r[k,i] = simplify(r[k,i])
- if r[k,i] != 0:
+ if not iszerofunc(r[k,i]):
break
- if k == r.lines - 1 and r[k,i] == 0:
+ if k == r.lines - 1 and iszerofunc(r[k,i]):
continue
r.row_swap(pivots,k)
scale = r[pivots,i]
diff -r 4d0cd200f57b -r 5fc85c34d2ae sympy/matrices/tests/test_matrices.py
--- a/sympy/matrices/tests/test_matrices.py Wed Nov 05 13:18:07 2008 +0100
+++ b/sympy/matrices/tests/test_matrices.py Thu Nov 06 22:23:24 2008 +0100
@@ -938,3 +938,9 @@
assert Matrix([[(exp(x)-1)/x, 2*x + y*x, x**x ],
[1/x, abs(x) , abs(sin(x+1))]]).limit(x, 0) == Matrix([[1,
0, 1],[oo, 0, sin(1)]])
assert a.integrate(x) == Matrix([[Rational(1,3)*x**3,
y*x**2/2],[x**2*sin(y)/2, x**2*cos(y)/2]])
+
+def test_inv_iszerofunc():
+ A = eye(4)
+ A.col_swap(0,1)
+ for method in "GE", "LU":
+ assert A.inv(method, iszerofunc=lambda x: x==0) == A.inv("ADJ")
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy-commits" group.
To post to this group, send email to sympy-commits@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sympy-commits?hl=en
-~----------~----~----~----~------~----~------~--~---
