Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
##############################################################################
# ------------------------------nonlinsolve ---------------------------------#
##############################################################################
def substitution(system, symbols=None, result=[{}], known_symbols=[], exclude=[],
all_symbols=None):
r""" Solves the `system` using substitution method.
A helper function for `nonlinsolve`. This will be called from
`nonlinsolve` when any equation(s) is non polynomial equation.
Parameters
==========
system : list of equations
The target system of equations
symbols : list of symbols to be solved.
The variable(s) for which the system is solved
known_symbols : list of solved symbols
Values are known for these variable(s)
result : An empty list or list of dict
If No symbol values is known then empty list otherwise
symbol as keys and corresponding value in dict.
exclude : Set of expression.
Mostly denominator expression(s) of the equations of the system.
Final solution should not satisfy these expressions.
all_symbols : known_symbols + symbols(unsolved).
Returns
=======
A FiniteSet of ordered tuple of values of `all_symbols` for which the
`system` has solution. Order of value is same as symbols order in the
`all_symbols` list.
Please note that general FiniteSet is unordered, the solution returned
here is not simply a FiniteSet of solutions, rather it is a FiniteSet of
ordered tuple, i.e. the first & only argument to FiniteSet is a tuple of
solutions, which is ordered, & hence the returned solution is ordered.
Also note that solution could also have been returned as an ordered tuple,
FiniteSet is just a wrapper `{}` around the tuple. It has no other
significance except for the fact it is just used to maintain a consistent
output format throughout the solveset.
Examples
========
>>> from sympy.core.symbol import symbols
>>> x, y = symbols('x, y', real = True)
>>> from sympy.solvers.solveset import substitution
>>> substitution([x +y], [x], [{y: 1}], [y], set([]), [x, y])
{(-1, 1)}
* when you want soln should not satisfy eq `x + 1 = 0`
>>> substitution([x +y], [x], [{y: 1}], [y], set([x + 1]), [y, x])
EmptySet()
>>> substitution([x +y], [x], [{y: 1}], [y], set([x - 1]), [y, x])
{(1, -1)}
>>> substitution([x + y - 1, y - x**2 + 5], [x, y])
{(-3, 4), (2, -1)}
* Returns both real and complex solution
>>> x, y, z = symbols('x, y, z')
>>> from sympy import exp, sin
>>> substitution([exp(x) - sin(y), y**2 - 4], [x, y])
{(log(sin(2)), 2), (ImageSet(Lambda(_n, I*(2*_n*pi + pi) +
log(sin(2))), Integers()), -2), (ImageSet(Lambda(_n, 2*_n*I*pi +
Mod(log(sin(2)), 2*I*pi)), Integers()), 2)}
>>> eqs = [z**2 + exp(2*x) - sin(y), -3 + exp(-y)]
>>> substitution(eqs, [y, z])
{(-log(3), -sqrt(-exp(2*x) - sin(log(3)))),
(-log(3), sqrt(-exp(2*x) - sin(log(3)))),
(ImageSet(Lambda(_n, 2*_n*I*pi + Mod(-log(3), 2*I*pi)), Integers()),
ImageSet(Lambda(_n, -sqrt(-exp(2*x) + sin(2*_n*I*pi +
Mod(-log(3), 2*I*pi)))), Integers())),
(ImageSet(Lambda(_n, 2*_n*I*pi + Mod(-log(3), 2*I*pi)), Integers()),
ImageSet(Lambda(_n, sqrt(-exp(2*x) + sin(2*_n*I*pi +
Mod(-log(3), 2*I*pi)))), Integers()))}
"""
from sympy.core.compatibility import ordered, default_sort_key
from sympy import Complement
if not system:
return S.EmptySet
if not symbols:
raise ValueError('Symbols must be given, for which solution of the '
'system is to be found.')
if hasattr(symbols[0], '__iter__'):
symbols = symbols[0]
try:
sym = symbols[0].is_Symbol
except AttributeError:
sym = False
if not sym:
raise ValueError('Symbols or iterable of symbols must be given as '
'second argument, not type \
%s: %s' % (type(symbols[0]), symbols[0]))
# By default `all_symbols` will same as `symbols`
if all_symbols is None:
all_symbols = symbols
def _unsolved_syms(eq, sort=False):
"""Returns the unsolved symbol present
in the equation `eq`.
"""
free = eq.free_symbols
unsolved = (free - set(known_symbols)) & set(all_symbols)
if sort:
unsolved = list(unsolved)
unsolved.sort(key=default_sort_key)
return unsolved
old_result = result
# storing complements and intersection for particular symbol
complements = {}
intersections = {}
def _solve_using_known_values(result, solver):
"""Solves the system using already known solution
(result contains the dict <symbol: value>).
solver is `solveset_complex` or `solveset_real`.
"""
# sort such that equation with the fewest potential symbols is first.
# means eq with less variable first.
soln_imageset = None # stores imageset if it finds.
for eq in ordered(system, lambda _: len(_unsolved_syms(_))):
newresult = []
got_symbol = set() # symbols solved in one iteration
hit = False
original_imageset = {}
for res in result:
if soln_imageset:
# find the imageset and use it's expr.
for key_res, value_res in res.items():
if isinstance(value_res, ImageSet):
res[key_res] = value_res.lamda.expr
original_imageset[key_res] = value_res
soln_imageset = None
# update eq with everything that is known so far
eq2 = eq.subs(res)
unsolved_syms = _unsolved_syms(eq2, sort=True)
if not unsolved_syms:
if res:
newresult.append(res)
break # skip as it's independent of desired symbols
for sym in unsolved_syms:
not_solvable = False
try:
soln = solver(eq2, sym)
except NotImplementedError:
continue
if isinstance(soln, Complement):
# extract solution and complement
complements[sym] = soln.args[1]
soln = soln.args[0]
# complement should be added at the end
if isinstance(soln, Intersection):
# sometimes solveset returns Intersection with S.Real
# when #11174 is fixed,then this line should be removed
# Interval will be at 0th index
if soln.args[0] != Interval(-oo, oo):
intersections[sym] = soln.args[0]
soln = soln.args[1]
if isinstance(soln, ImageSet):
soln_imageset = soln
soln = FiniteSet(soln.lamda.expr)
elif isinstance(soln, ConditionSet):
soln = S.EmptySet
# dont do `continue` we may get soln
# in terms of other symbol(s)
not_solvable = True
for sol in soln:
free = sol.free_symbols
if got_symbol and any([
ss in free for ss in got_symbol
]):
# sol depends on previously solved symbols
# then continue
continue
rnew = res.copy()
# put each solution in res and append the new result
# in the new result list (solution for symbol `s`)
# along with old results.
for k, v in res.items():
if isinstance(v, Expr):
# if any unsolved symbol is present
# Then subs known value
rnew[k] = v.subs(sym, sol)
# and add this new solution
rnew[sym] = sol
# for check solve use imageset expr
if not any(checksol(d, rnew) for d in exclude):
# if sol was imageset then add imageset
if soln_imageset:
rnew[sym] = soln_imageset
# restore original imageset
restore_sym = set(rnew.keys()) & \
set(original_imageset.keys())
dummy_n = None
base_img = None
for key_sym in restore_sym:
img = original_imageset[key_sym]
rnew[key_sym] = img
dummy_n = (img.lamda.expr).atoms(Dummy).pop()
base_img = img.base_set
if dummy_n and base_img:
# also it means soln for `sym` contains
# dummy `n`, put this expr in ImageSet
lam = Lambda(dummy_n, rnew[sym])
rnew[sym] = ImageSet(lam, base_img)
# now append this solution
newresult.append(rnew)
hit = True
# solution got for s
if not not_solvable:
got_symbol.add(sym)
if not hit:
raise NotImplementedError('could not solve %s' % eq2)
else:
result = newresult
return result
# end def _solve_using_know_values
new_result_real = _solve_using_known_values(old_result, solveset_real)
new_result_complex = _solve_using_known_values(
old_result, solveset_complex)
# overall result
result = new_result_real + new_result_complex
result_finiteset = []
for res in result:
if not res:
# if {None : None} is present.
# No solution then first will be {None : None}
continue
# If length < len(all_symbols) means infinite soln.
# Some or all the soln is dependent on 1 symbol.
# eg. {x: y+2} then final soln is {x: y+2, y: y}
if len(res) < len(all_symbols):
solved_symbols = res.keys()
unsolved = list(filter(
lambda x: x not in solved_symbols, all_symbols))
rcopy = res.copy()
for unsolved_sym in unsolved:
rcopy[unsolved_sym] = unsolved_sym
# if symbol is not in `v` (unsolved symbol),
# means it can take any value.
# eg. if k, v => y, exp(x) then above lines will add {x: x}
# if we have another symbol `z` then add {z: z} using above line.
res = rcopy
result_finiteset.append(res)
if complements:
# If solveset have returned some complements for any symbol
# Mostly complements are Set.
result = []
for res in result_finiteset:
res_copy = res
for key_res, value_res in res.items():
for key_complement, value_complement in complements.items():
if key_complement == key_res:
res_copy[key_res] = \
FiniteSet(value_res) - value_complement
result.append(res_copy)
result_finiteset = result
if intersections:
# If solveset have returned some intersection for any symbol
# mostly intersection is Interval or Set
result = []
for res in result_finiteset:
res_copy = res
for key_res, value_res in res.items():
for key_intersec, value_intersec in intersections.items():
if key_intersec == key_res:
res_copy[key_res] = \
Intersection(FiniteSet(value_res), value_intersec)
result.append(res_copy)
result_finiteset = result
# convert to ordered tuple
result = S.EmptySet
for r in result_finiteset:
temp = [r[sym] for sym in all_symbols]
result += FiniteSet(tuple(temp))
return result
def nonlinsolve(system, symbols):
r"""
Solve system of N non linear equations with M variables, which means both
under - and overdetermined systems are supported. Positive dimensional
system is also supported (A system with infinitely many solutions is said
to be positive-dimensional). In Positive dimensional system solution will
be dependent on at least one symbol. If system doesn't have real solution
then complex solution will be returned in general form (in terms of `_n`,
where `_n` is any real number). The possible number of solutions is zero,
one or infinite.
Parameters
==========
system : list of equations
The target system of equations
symbols : list of Symbols
symbols should be given as a sequence eg. list
Returns
=======
A FiniteSet of ordered tuple of values of `symbols` for which the `system`
has solution. Order of value is same as symbols order in the
`all_symbols`list.
Please note that general FiniteSet is unordered, the solution returned
here is not simply a FiniteSet of solutions, rather it is a FiniteSet of
ordered tuple, i.e. the first & only argument to FiniteSet is a tuple of
solutions, which is ordered, & hence the returned solution is ordered.
Also note that solution could also have been returned as an ordered tuple,
FiniteSet is just a wrapper `{}` around the tuple. It has no other
significance except for the fact it is just used to maintain a consistent
output format throughout the solveset.
For the given set of Equations, the respective input types
are given below:
.. math:: x*y - 1 = 0
.. math:: 4*x**2 + y**2 - 5 = 0
`system = [x*y - 1, 4*x**2 + y**2 - 5]`
`symbols = [x, y]`
Examples
========
>>> from sympy.core.symbol import symbols
>>> from sympy.solvers.solveset import nonlinsolve
>>> x, y, z = symbols('x, y, z', real=True)
>>> nonlinsolve([x*y - 1, 4*x**2 + y**2 - 5], [x, y])
{(-1, -1), (-1/2, -2), (1/2, 2), (1, 1)}
* Positive dimensional system and complements:
>>> from sympy import pprint
>>> from sympy.polys.polytools import is_zero_dimensional
>>> a, b, c, d = symbols('a, b, c, d', real=True)
>>> foo = a + b + c + d
>>> bar = a*b + b*c + c*d + d*a
>>> foo_bar = a*b*c + b*c*d + c*d*a + d*a*b
>>> bar_foo = a*b*c*d -1
>>> system = [foo, bar, foo_bar, bar_foo]
>>> is_zero_dimensional(system)
False
>>> pprint(nonlinsolve(system, [a, b, c, d]), use_unicode=False)
-1 1
{(---, -d, -, {d} \ {0})}
d d
>>> nonlinsolve([(x+y)**2 - 4, x + y - 2], [x, y])
{(-y + 2, y)}
* If Non linear system having non polynomial equation have complex
solution for particular symbol and no real solution, then `nonlinsolve`
returns it's general solution.
>>> from sympy import exp, sin
>>> _n = symbols('_n')
>>> nonlinsolve([exp(x) - sin(y), y**2 - 4], [x, y])
{(log(sin(2)), 2), (ImageSet(Lambda(_n, I*(2*_n*pi + pi) +
log(sin(2))), Integers()), -2), (ImageSet(Lambda(_n, 2*_n*I*pi +
Mod(log(sin(2)), 2*I*pi)), Integers()), 2)}
* If Non linear system have all the equations polynomial, then it
returns both real and complex solutions.
>>> from sympy import sqrt
>>> nonlinsolve([x**2 - 2*y**2 -2, x*y - 2], [x, y])
{(-2, -1), (2, 1), (-sqrt(2)*I, sqrt(2)*I), (sqrt(2)*I, -sqrt(2)*I)}
* If system is linear and positive dimensional system, `nonlinsolve` can
solve this system (because it is using `groebner` function to get the
groebner basis and then `substitution` function basis as system). But it is
not recommended to solve linear system using `nonlinsolve`, because
`linsolve` is better for all kind of linear system.
>>> nonlinsolve([x + 2*y -z - 3, x - y - 4*z + 9 , y + z - 4], [x, y, z])
{(3*z - 5, -z + 4, z)}
* system having polynomial equations:
>>> e1 = sqrt(x**2 + y**2) - 10
>>> e2 = sqrt(y**2 + (-x + 10)**2) - 3
>>> nonlinsolve((e1, e2), (x, y))
{(191/20, -3*sqrt(391)/20), (191/20, 3*sqrt(391)/20)}
>>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [x, y])
{(1, 2), (1 + sqrt(5), -sqrt(5) + 2), (-sqrt(5) + 1, 2 + sqrt(5))}
>>> nonlinsolve([x**2 + 2/y - 2, x + y - 3], [y, x])
{(2, 1), (2 + sqrt(5), -sqrt(5) + 1), (-sqrt(5) + 2, 1 + sqrt(5))}
How nonlinsolve is better than old solver `_solve_system` :
===========================================================
* A positive dimensional system solver : nonlinsolve can return
solution for positive dimensional system. It finds the
Groebner Basis of the positive dimensional system(calling it as
basis) then we can start solving equation(having least number of
variable first in the basis) using solveset and substituting that
solved solutions into other equation(of basis) to get solution in
terms of minimum variables. Here the important thing is how we
are substituting the known values and in which equations.
* Real and Complex both solutions : nonlinsolve returns solution both real
and complex solution. If all the equations in the system are polynomial
then using `solve_poly_system` both real and complex solution is returned.
If all the equations in the system are not polynomial equation then solves
Polynomial equations(If any) using `solve_poly_system` and goes to
`substitution` method with this solution and non polynomial equation(s),
to solve for unsolved variables. Here to solve for particular variable
solveset_real and solveset_complex is used. For both real and complex
solution function `_solve_using_know_values` is used inside `substitution`
function.(`substitution` function will be called when there is non
polynomial equation(s) is present). When solution is valid then add its
general solution in the final result. There may the case when
`solve_poly_system` will not give complex solution (only real solution)
then that soln will be missing.If we want to get complex solution in
general form then use `substitution` function.(it always returns all the
solutions).
* Complements and Intersection will be added if any : nonlinsolve maintains
dict for complements and Intersections. If solveset find complements or
Intersection with any Interval or set in during the execution of
`substitution` function ,then complement and Intersection for that
variable is added before returning final solution.
"""
from sympy.solvers.solvers import _invert as _invert_solver
from sympy.utilities.iterables import subsets
from sympy.polys.polytools import is_zero_dimensional, groebner
if not system:
return S.EmptySet
if not symbols:
raise ValueError('Symbols must be given, for which solution of the '
'system is to be found.')
if hasattr(symbols[0], '__iter__'):
symbols = symbols[0]
try:
sym = symbols[0].is_Symbol
except AttributeError:
sym = False
if not sym:
raise ValueError('Symbols or iterable of symbols must be given as '
'second argument, not type \
%s: %s' % (type(symbols[0]), symbols[0]))
polys = []
nonpolys = []
denominators = set()
for eq in system:
# Store denom expression if it contains symbol
denominators.update(_simple_dens(eq, symbols))
# TODO : solveset `_invert is not for multiple symbol.
# using old solver's _invert for now
independent, dependent = _invert_solver(eq, *symbols)
eq = dependent - independent
eq = eq.as_numer_denom()[0]
# this will remove sqrt if symbols are under it
poly = eq.as_poly(*symbols, extension=True)
if poly is not None:
polys.append(poly)
else:
nonpolys.append(eq)
# If none of the equation is Poly
if not polys:
solved_syms = []
result = None
if len(symbols) == len(polys):
# If all the equations in the system is poly
if is_zero_dimensional(system):
try:
# try to solve 0 dimensional system using `solve_poly_system`
result = solve_poly_system(polys, *symbols)
# May be some extra soln is added because
# we did `poly = g.as_poly(*symbols, extension=True)`
# need to check that and remove if not a soln
result_update = S.EmptySet
for res in result:
dict_sym_value = dict(list(zip(symbols, res)))
if all(checksol(eq, dict_sym_value) for eq in system):
result_update += FiniteSet(res)
return result_update
except:
# Right now it doesn't fail for any polynomial system of
# equation if `solve_poly_system` fails then substitution
# method will handle it. pass it
pass
# positive dimensional system
# substitution method with groebner basis of the system
basis = groebner(polys, symbols, polys=True)
new_system = []
for poly_eq in basis:
new_system.append(poly_eq.as_expr())
# solved_symbols = []
result = [{}]
result = substitution(
new_system, symbols, result, [],
denominators, symbols)
else:
# If alll the equations are not polynomial.
result = []
depend_soln = {}
solved_syms = []
# first solve the polynomial equation(s) if present
if polys:
combinations = list(subsets(symbols, len(polys)))
# solving for combinations of symbols. So that we get
# solution in terms of other symbol(s).
for new_symbols in combinations:
try:
# solution for new_symbols in terms of other symbols
result_poly = solve_poly_system(polys, new_symbols)
# TODO : improve the result code so that we
# dont need to add r[unsolved symbols] = unsolved
# at the end ie. if we got x value then do
# FiniteSet(x_val, y) and update(dynamically update)
# the `y` when got the
# value otherwise return as it is(for +ve dimensional).
for res in result_poly:
skip = False
# check tuples of res
for res1 in res:
if depend_soln and \
any(
[k in res1.free_symbols
for k, v in depend_soln.items()]):
# sol depends on previously
# solved symbols: discard it
# eg : depend_soln=> {x : -y} and
# r1 = -x for y
# so we dont want both {x: -y}, {y: -x}
skip = True
if not skip:
for new_sym in new_symbols:
for res1 in res:
depend_soln[new_sym] = res1
# need to maintain dict <symbol : value>
dict_sym_value = dict(list(zip(new_symbols, res)))
result.append(dict_sym_value)
result_update = []
for res in result:
# If length < len(symbols) means infinite soln
# (one symbol soln is in terms of other symbol).
# Some or all the soln is dependent on 1 symbol
if len(res) < len(symbols):
rcopy = res.copy()
unsolved = FiniteSet(*[sym2 for sym2 in
list(filter(lambda x: x not in res.keys(), symbols))
])
for unsolv in unsolved:
rcopy[unsolv] = unsolv
# remove soln which is making denominator zero
if not any(
checksol(d, rcopy) for d in denominators
):
result_update.append(rcopy)
else:
# add if res isn't making denominator expr zero
if not any(
checksol(d, res) for d in denominators
):
result_update.append(res)
result = result_update
except:
pass
if depend_soln:
solved_syms = list(depend_soln)
# If we have solved all the symbols or can't solve further(no more
# eq), return that soln in ordered tuple
if not nonpolys:
# convert it to ordered tuple and return
result_tuple = S.EmptySet
for res in result:
rcopy = res.copy()
temp = [rcopy[s] for s in symbols]
result_tuple += FiniteSet(tuple(temp))
return result_tuple
# Polynomial eqs done. Now use substitution method to get the solution
# for unsolved symbols if non polynomial eq is there..
if nonpolys:
# if non polynomail equation is present
if not result:
# because we are passing list of dict
result = [{}]
unsolved_syms = list(
filter(lambda x: x not in solved_syms, symbols))
# Use substitution method for nonpoly eq
result = substitution(
nonpolys, unsolved_syms, list(result),
solved_syms, denominators, symbols
)
if result:
return result
else:
return S.EmptySet
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment