Skip to content

Instantly share code, notes, and snippets.

# Shekharrajak/nonlinsolve_till_24june_2016.py Created Jun 24, 2016

 ############################################################################## # ------------------------------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, '__iter__'): symbols = symbols try: sym = symbols.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), symbols)) # 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 ). solver is `solveset_complex` or `solveset_real`. """ # stores imageset . soln_imageset = {} # sort such that equation with the fewest potential symbols is first. # means eq with less variable first for index, eq in enumerate( ordered(system, lambda _: len(_unsolved_syms(_)))): newresult = [] got_symbol = set() # symbols solved in one iteration hit = False original_imageset = {} # if imageset expr is used to solve other symbol imgset_yes = False 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 dummy_n = value_res.lamda.expr.atoms(Dummy).pop() base = value_res.base_set imgset_yes = (dummy_n, base) soln_imageset = {} # 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) soln_new = S.EmptySet if isinstance(soln, Complement): # extract solution and complement complements[sym] = soln.args soln = soln.args # complement will be added at the end if isinstance(soln, Intersection): # sometimes solveset returns Intersection with # S.Reals when #11174 is fixed,then # this if block should be removed # Interval will be at 0th index if soln.args != Interval(-oo, oo): intersections[sym] = soln.args soln_new += soln.args soln = soln_new if soln_new else soln if index > 0 and solver == solveset_real: # 1 symbol's real soln , another symbol may have # corresponding complex soln. if not isinstance(soln, ImageSet): soln += solveset_complex(eq2, sym) except NotImplementedError: continue if isinstance(soln, ConditionSet): soln = S.EmptySet # dont do `continue` we may get soln # in terms of other symbol(s) not_solvable = True if isinstance(soln, Complement): # extract solution and complement complements[sym] = soln.args soln = soln.args # complement will be added at the end if isinstance(soln, ImageSet): soln_imagest = soln expr2 = soln.lamda.expr soln = FiniteSet(expr2) soln_imageset[expr2] = soln_imagest # if there is union of Imageset in soln # no testcase for this if block if isinstance(soln, Union): soln_args = soln.args soln = [] # need to use list. # (finiteset union imageset).args is # (finiteset, Imageset). But we want only iteration # so append finteset elements and imageset. for soln_arg2 in soln_args: if isinstance(soln_arg2, FiniteSet): list_finiteset = list( [sol_1 for sol_1 in soln_arg2]) soln += list_finiteset else: # ImageSet, Intersection, complement # append directly soln.append(soln_arg2) for sol in soln: # if there is union, then need to check # complement, intersection, Imageset # order should not be changed. if isinstance(soln, Complement): # extract solution and complement complements[sym] = soln.args soln = soln.args # complement will be added at the end if isinstance(sol, Intersection): # sometimes solveset returns Intersection with # S.Reals when #11174 is fixed,then # this if block should be removed # Interval will be at 0th index if sol.args != Interval(-oo, oo): intersections[sym] = sol.args sol = sol.args # after intersection and complement Imageset should # be checked. if isinstance(sol, ImageSet): soln_imagest = sol sol = sol.lamda.expr soln_imageset[sol] = soln_imagest 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 local_n = None if imgset_yes: local_n = imgset_yes base = imgset_yes # use ImageSet, we have dummy in sol lam = Lambda(local_n, sol) rnew[sym] = ImageSet(lam, base) if soln_imageset: rnew[sym] = soln_imageset[sol] img_expr = soln_imageset[sol].lamda.expr local_n = img_expr.atoms(Dummy).pop() # 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 expr = img.lamda.expr dummy_n = expr.atoms(Dummy).pop() base_img = img.base_set if local_n and restore_sym: # also it means soln for `sym` contains # dummy `n`, put this expr in ImageSet # use only one dummy img_expr = rnew[sym].lamda.expr lam = Lambda( dummy_n, img_expr.subs(local_n, dummy_n) ) 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 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 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 # 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. Returns both real solution and complex solution. 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) >>> eq1 = a + b + c + d >>> eq2 = a*b + b*c + c*d + d*a >>> eq3 = a*b*c + b*c*d + c*d*a + d*a*b >>> eq4 = a*b*c*d -1 >>> system = [eq1, eq2, eq3, eq4] >>> 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 some of the equations are non polynomial eq then `nonlinsolve` will call `substitution` function and returns real and complex solutions . >>> from sympy import exp, sin >>> 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 polynomial and zero dimensional system, then it returns both real and complex solutions using `solve_poly_system`. >>> 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)} * `nonlinsolve` can solve some linear( zero or positive dimensional)system solve this system (because it is using `groebner` function to get the groebner basis and then `substitution` function basis as the 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.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, '__iter__'): symbols = symbols try: sym = symbols.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), symbols)) polys = [] polys_expr = [] 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() # this will remove sqrt if symbols are under it poly = eq.as_poly(*symbols, extension=True) if poly is not None: polys.append(poly) polys_expr.append(poly.as_expr()) else: nonpolys.append(eq) 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. # Use substitution method for nonpoly + polys eq result = substitution( polys_expr + nonpolys, symbols, exclude=denominators) if result: return result else: return S.EmptySet
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.