Comment #6 on issue 1850 by [email protected]: solve_linear_system contains
duplicate rref algorithm
http://code.google.com/p/sympy/issues/detail?id=1850
solve_linear_system() is not user-friendly, and has un-conventional
behavior. It and the related functions for solving linear
equations that should be replaced.
In many elementary algebraic settings, one can obtain a system of
linear equations that needs solving. It would be great if we could
just call solve_linear_system(equations, variables) and get our
solution. However, solve_linear_system() assumes our system is in an
implicit matrix form. We can obtain solutions using
solve(equations,variables) and solve_poly_system(equations,variables).
However, this relies on both functions being smart enough to recognize
the special linear properties of the system, which we are often
already aware of. We could also write code to transform our system to
matrix form for solve_linear_system(), but seems like much
un-necessary work.
It seems natural to me that solve_linear_system() or some equivalent
function should be able to recognize and pre-process all natural forms
of linear systems.
A related issue is that the behaviors of solve_linear_system and
solve_poly_system are inconsistent.
In [25]: from sympy import *
In [26]: from sympy.abc import x,y
In [27]: m = Matrix([[1,2,0],[3,-4,1]])
In [28]: solve_linear_system(m,x,y)
Out[28]: {x: 1/5, y: -1/10}
In [29]: solve_linear_system(m,[x,y])
Out[29]: []
In [32]: e = [x+2*y,3*x-4*y-1]
In [33]: solve_poly_system(e,x,y)
Out[33]: [(1/5, -1/10)]
In [34]: solve_poly_system(e,[x,y])
Out[34]: [(1/5, -1/10)]
Both approaches have merit, and return the correct answer.
solve_linear_system's returning of a dictionary makes potential
assignments explicit. However, it fails when the variables are given
inside an iterable. The output of solve_poly_system() is consistent.
And solve_poly_system also tolerates different input formats.
We further note that solve_linear_system() assumes that our matrix is
in the form of an augmented matrix, where the last column is on the
right of the equals sign when writen in matrix-vector form. However,
there is no way to tell that this is an augmented matrix rather than a
regular matrix, just from inspection of the input. Users might
naturally assume that solve_linear_system calculates the nullspace of
the given matrix, which would return the opposite sign of the correct
answer. In fact, solve_linear_system() can be replaced with a 1-liner
transforming the matrix to reduced row-eschelon form.
In [35]: dict(zip((x,y),list(m.rref()[0][:,-1])))
Out[35]: {x: 1/5, y: -1/10}
There is a 3rd function, solve_linear_system_LU,
In [37]: solve_linear_system_LU(m,[x,y])
Out[37]: {x: 1/5, y: -1/10}
solve_linear_system_LU() has it's own idiosyncratic behavior, as it only
accepts iterables of variables.
In [38]: solve_linear_system_LU(m,x,y)
---------------------------------------------------------------------------
TypeError Traceback (most recent call
last)
/home/junk/src/sympy/bin/<ipython-input-38-4683da0e9657> in <module>()
----> 1 solve_linear_system_LU(m,x,y)
TypeError: solve_linear_system_LU() takes exactly 2 arguments (3 given)
The 4th function, solve_linear(), seems to be intended for use in the
solution of systems of polynomial equations. However, it also has
non-generic behavior, sometimes returning numerator/denominator of
the argument, and sometimes returning a solution of a linear equation.
The exact reason for this behavior remains unclear to me, and I'd
appreciate anybody who can help me out in understanding this. To me,
the name atleast appears misleading.
Currently, all 4 functions are used sparesly within sympy's code base.
This is a surprise, in itself, given the frequency with which linear
systems are employed in mathematics, and suggests that they are not
fulfilling the niches they are named for. A consequence is that a
project-wide change would be a managable undertaking at this point in
time. The potential impact on existing code bases outside the project
is unclear.
There are several potential options.
1) Standardize on one particular form of argument style for solvers,
and reimplement. Deprecate obsolete versions.
2) Add code so that the existing solvers allow variety of argument styles.
3) Leave as-is. Pros: no changes needed. Cons: New users may be
confused by different styles, and may be less likely to adopt for
regular usage. Users who do adopt will be channeled toward "solve()"
which may be in-efficient in many simple cases. As usage spreads,
design changes breaking backward compatibility will get more
difficult.
The code below is my proposal for a replacement, and I am looking for
feedback. The function is just a wrapper around rref(), which does
all the hard work.
from sympy import *
def solve_lin_sys(eqs,xs):
"""
eqs is an iterable of linear equations
maybe Eq(*,*) or expr
if expr, we assume it is implicitly
set to zero.
xs is an iterable of variables we are solving for
"""
if isinstance(xs,Symbol):
xs = (xs,)
else:
xs = tuple(xs,)
if isinstance(eqs,Matrix):
assert eqs.shape[1] <= len(xs)+1, "too many matrix columns for
variables"
assert eqs.shape[1] >= len(xs)+1, "too few matrix columns for
variables"
else:
# move all equation arguments to one side
if not hasattr(eqs, "__iter__"):
eqs = ( eqs, )
new_eqs = []
for e in eqs:
if isinstance(e,Equality):
e = e.lhs-e.rhs
for x in xs:
assert diff(e,x,2) == 0, "not a linear system"
new_eqs.append(e)
eqs = tuple(new_eqs)
# transform from equations to matrix form
S = [ (x,0) for x in xs ]
M = zeros(len(eqs), len(xs)+1)
for j,e_j in enumerate(eqs):
for i,x_i in enumerate(xs):
M[j,i] = diff(e_j,x_i)
M[j,-1] = -1*e_j.subs(S)
eqs = M
# solve by row-reduction
eschelon, pivots = eqs.rref()
# construct the returnable form of the solutions
p = len(pivots)
if p > len(xs):
return None
sols = eschelon[0:p,p:]*Matrix([ [x] for x in xs[p:] ]+[ [1] ])
return dict(zip(xs,sols))
def example():
from sympy.abc import x,y,z
eqs = [x + y - 3, x - y - 4]
vs = [x,y]
# works for a list of expressions and a list of variables
print solve_lin_sys(eqs,vs)
# works for a list of equations
eqs = ( Eq(x + y, 3), Eq(x, y + 4) )
print solve_lin_sys(eqs,vs)
# works for a matrix and a list of variables
eqs = Matrix([[1,1,3],[1,-1,4]])
print solve_lin_sys(eqs,vs)
# works for simple equations
print solve_lin_sys(x+y**2 + 3*x,x)
# for over-determined systems
eqs = Matrix([[1,1,3],[1,-1,4],[2,-1,-5],[4,4,-2]])
print solve_lin_sys(eqs,vs)
# for under-determined systems, leave the last variables free
eqs = Matrix([[1,1,0,3],[1,-1,1,4]])
print solve_lin_sys(eqs,[x,y,z])
example()
--
You received this message because you are subscribed to the Google Groups
"sympy-issues" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/sympy-issues?hl=en.
