Coverage for sympy/solvers/solvers.py : 90%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
|
""" This module contain solvers for all kinds of equations:
- algebraic or transcendental, use solve()
- recurrence, use rsolve()
- differential, use dsolve()
- nonlinear (numerically), use nsolve() (you will need a good starting point)
"""
default_sort_key, range) Pow, Unequality) Derivative, AppliedUndef, UndefinedFunction, nfloat, Function, expand_power_exp, Lambda, _mexpand)
Abs, re, im, arg, sqrt, atan2) HyperbolicFunction) nsimplify, denom, logcombine)
"""Return (e, s, d) where e and s are versions of eqs and symbols in which any non-Symbol objects in symbols have been replaced with generic Dummy symbols and d is a dictionary that can be used to restore the original expressions.
Examples ========
>>> from sympy.solvers.solvers import recast_to_symbols >>> from sympy import symbols, Function >>> x, y = symbols('x y') >>> fx = Function('f')(x) >>> eqs, syms = [fx + 1, x, y], [fx, y] >>> e, s, d = recast_to_symbols(eqs, syms); (e, s, d) ([_X0 + 1, x, y], [_X0, y], {_X0: f(x)})
The original equations and symbols can be restored using d:
>>> assert [i.xreplace(d) for i in eqs] == eqs >>> assert [d.get(i, i) for i in s] == syms """ raise ValueError('Both eqs and symbols must be iterable')
"""Return True if e is a Pow or is exp."""
# when checking if a denominator is zero, we can just check the # base of powers with nonzero exponents since if the base is zero # the power will be zero, too. To keep it simple and fast, we # limit simplification to exponents that are Numbers
"""Return (recursively) set of all denominators that appear in eq that contain any symbol in ``symbols``; if ``symbols`` are not provided then all denominators will be returned.
Examples ========
>>> from sympy.solvers.solvers import denoms >>> from sympy.abc import x, y, z >>> from sympy import sqrt
>>> denoms(x/y) {y}
>>> denoms(x/(y*z)) {y, z}
>>> denoms(3/x + y/z) {x, z}
>>> denoms(x/2 + y/z) {2, z}
If `symbols` are provided then only denominators containing those symbols will be returned
>>> denoms(1/x + 1/y + 1/z, y, z) {y, z} """
"""Checks whether sol is a solution of equation f == 0.
Input can be either a single symbol and corresponding value or a dictionary of symbols and values. When given as a dictionary and flag ``simplify=True``, the values in the dictionary will be simplified. ``f`` can be a single equation or an iterable of equations. A solution must satisfy all equations in ``f`` to be considered valid; if a solution does not satisfy any equation, False is returned; if one or more checks are inconclusive (and none are False) then None is returned.
Examples ========
>>> from sympy import symbols >>> from sympy.solvers import checksol >>> x, y = symbols('x,y') >>> checksol(x**4 - 1, x, 1) True >>> checksol(x**4 - 1, x, 0) False >>> checksol(x**2 + y**2 - 5**2, {x: 3, y: 4}) True
To check if an expression is zero using checksol, pass it as ``f`` and send an empty dictionary for ``symbol``:
>>> checksol(x**2 + x - x*(x + 1), {}) True
None is returned if checksol() could not conclude.
flags: 'numerical=True (default)' do a fast numerical check if ``f`` has only one symbol. 'minimal=True (default is False)' a very fast, minimal testing. 'warn=True (default is False)' show a warning if checksol() could not conclude. 'simplify=True (default)' simplify solution before substituting into function and simplify the function before trying specific simplifications 'force=True (default is False)' make positive all symbols without assumptions regarding sign.
"""
else: msg = 'Expecting (sym, val) or ({sym: val}, None) but got (%s, %s)' raise ValueError(msg % (symbol, sol))
if not f: raise ValueError('no functions to check') rv = True for fi in f: check = checksol(fi, sol, **flags) if check: continue if check is False: return False rv = None # don't return, wait to see if there's a False return rv
f = f.as_expr() return else:
return True
# if f(y) == 0, x=3 does not set f(y) to zero...nor does it not
S.ComplexInfinity, S.Infinity, S.NegativeInfinity])
# there are free symbols -- simple expansion might work # start over without the failed expanded form, possibly # with a simplified solution # expansion may work now, so try again and check # we can decide now else: # if there are no radicals and no functions then this can't be # zero anymore -- can it? saw_pow_func = True # don't do a zero check with the positive assumptions in place # issue 5673: nz may be True even when False # so these are just hacks to keep a false positive # from being returned
# HACK 1: LambertW (issue 5673) # don't eval this to verify solution since if we got here, # numerical must be False return None
# add other HACKs here if necessary, otherwise we assume # the nz value is correct
return bool(val)
warnings.warn("\n\tWarning: could not verify solution %s." % sol) # returns None if it can't conclude # TODO: improve solution testing
"""Return a dictionary containing assumptions with values not matching those of the passed assumptions.
Examples ========
>>> from sympy import failing_assumptions, Symbol
>>> x = Symbol('x', real=True, positive=True) >>> y = Symbol('y') >>> failing_assumptions(6*x + y, real=True, positive=True) {'positive': None, 'real': None}
>>> failing_assumptions(x**2 - 1, positive=True) {'positive': None}
If all assumptions satisfy the `expr` an empty dictionary is returned.
>>> failing_assumptions(x**2, positive=True) {} """
"""Checks whether expression `expr` satisfies all assumptions.
`assumptions` is a dict of assumptions: {'assumption': True|False, ...}.
Examples ========
>>> from sympy import Symbol, pi, I, exp, check_assumptions
>>> check_assumptions(-5, integer=True) True >>> check_assumptions(pi, real=True, integer=False) True >>> check_assumptions(pi, real=True, negative=True) False >>> check_assumptions(exp(I*pi/7), real=False) True
>>> x = Symbol('x', real=True, positive=True) >>> check_assumptions(2*x + 1, real=True, positive=True) True >>> check_assumptions(-2*x - 5, real=True, positive=True) False
To check assumptions of ``expr`` against another variable or expression, pass the expression or variable as ``against``.
>>> check_assumptions(2*x + 1, x) True
`None` is returned if check_assumptions() could not conclude.
>>> check_assumptions(2*x - 1, real=True, positive=True) >>> z = Symbol('z') >>> check_assumptions(z, real=True)
See Also ======== failing_assumptions """
r""" Algebraically solves equations and systems of equations.
Currently supported are: - polynomial, - transcendental - piecewise combinations of the above - systems of linear and polynomial equations - systems containing relational expressions.
Input is formed as:
* f - a single Expr or Poly that must be zero, - an Equality - a Relational expression or boolean - iterable of one or more of the above
* symbols (object(s) to solve for) specified as - none given (other non-numeric objects will be used) - single symbol - denested list of symbols e.g. solve(f, x, y) - ordered iterable of symbols e.g. solve(f, [x, y])
* flags 'dict'=True (default is False) return list (perhaps empty) of solution mappings 'set'=True (default is False) return list of symbols and set of tuple(s) of solution(s) 'exclude=[] (default)' don't try to solve for any of the free symbols in exclude; if expressions are given, the free symbols in them will be extracted automatically. 'check=True (default)' If False, don't do any testing of solutions. This can be useful if one wants to include solutions that make any denominator zero. 'numerical=True (default)' do a fast numerical check if ``f`` has only one symbol. 'minimal=True (default is False)' a very fast, minimal testing. 'warn=True (default is False)' show a warning if checksol() could not conclude. 'simplify=True (default)' simplify all but polynomials of order 3 or greater before returning them and (if check is not False) use the general simplify function on the solutions and the expression obtained when they are substituted into the function which should be zero 'force=True (default is False)' make positive all symbols without assumptions regarding sign. 'rational=True (default)' recast Floats as Rational; if this option is not used, the system containing floats may fail to solve because of issues with polys. If rational=None, Floats will be recast as rationals but the answer will be recast as Floats. If the flag is False then nothing will be done to the Floats. 'manual=True (default is False)' do not use the polys/matrix method to solve a system of equations, solve them one at a time as you might "manually" 'implicit=True (default is False)' allows solve to return a solution for a pattern in terms of other functions that contain that pattern; this is only needed if the pattern is inside of some invertible function like cos, exp, .... 'particular=True (default is False)' instructs solve to try to find a particular solution to a linear system with as many zeros as possible; this is very expensive 'quick=True (default is False)' when using particular=True, use a fast heuristic instead to find a solution with many zeros (instead of using the very slow method guaranteed to find the largest number of zeros possible) 'cubics=True (default)' return explicit solutions when cubic expressions are encountered 'quartics=True (default)' return explicit solutions when quartic expressions are encountered 'quintics=True (default)' return explicit solutions (if possible) when quintic expressions are encountered
Examples ========
The output varies according to the input and can be seen by example::
>>> from sympy import solve, Poly, Eq, Function, exp >>> from sympy.abc import x, y, z, a, b >>> f = Function('f')
* boolean or univariate Relational
>>> solve(x < 3) (-oo < x) & (x < 3)
* to always get a list of solution mappings, use flag dict=True
>>> solve(x - 3, dict=True) [{x: 3}] >>> sol = solve([x - 3, y - 1], dict=True) >>> sol [{x: 3, y: 1}] >>> sol[0][x] 3 >>> sol[0][y] 1
* to get a list of symbols and set of solution(s) use flag set=True
>>> solve([x**2 - 3, y - 1], set=True) ([x, y], {(-sqrt(3), 1), (sqrt(3), 1)})
* single expression and single symbol that is in the expression
>>> solve(x - y, x) [y] >>> solve(x - 3, x) [3] >>> solve(Eq(x, 3), x) [3] >>> solve(Poly(x - 3), x) [3] >>> solve(x**2 - y**2, x, set=True) ([x], {(-y,), (y,)}) >>> solve(x**4 - 1, x, set=True) ([x], {(-1,), (1,), (-I,), (I,)})
* single expression with no symbol that is in the expression
>>> solve(3, x) [] >>> solve(x - 3, y) []
* single expression with no symbol given
In this case, all free symbols will be selected as potential symbols to solve for. If the equation is univariate then a list of solutions is returned; otherwise -- as is the case when symbols are given as an iterable of length > 1 -- a list of mappings will be returned.
>>> solve(x - 3) [3] >>> solve(x**2 - y**2) [{x: -y}, {x: y}] >>> solve(z**2*x**2 - z**2*y**2) [{x: -y}, {x: y}, {z: 0}] >>> solve(z**2*x - z**2*y**2) [{x: y**2}, {z: 0}]
* when an object other than a Symbol is given as a symbol, it is isolated algebraically and an implicit solution may be obtained. This is mostly provided as a convenience to save one from replacing the object with a Symbol and solving for that Symbol. It will only work if the specified object can be replaced with a Symbol using the subs method.
>>> solve(f(x) - x, f(x)) [x] >>> solve(f(x).diff(x) - f(x) - x, f(x).diff(x)) [x + f(x)] >>> solve(f(x).diff(x) - f(x) - x, f(x)) [-x + Derivative(f(x), x)] >>> solve(x + exp(x)**2, exp(x), set=True) ([exp(x)], {(-sqrt(-x),), (sqrt(-x),)})
>>> from sympy import Indexed, IndexedBase, Tuple, sqrt >>> A = IndexedBase('A') >>> eqs = Tuple(A[1] + A[2] - 3, A[1] - A[2] + 1) >>> solve(eqs, eqs.atoms(Indexed)) {A[1]: 1, A[2]: 2}
* To solve for a *symbol* implicitly, use 'implicit=True':
>>> solve(x + exp(x), x) [-LambertW(1)] >>> solve(x + exp(x), x, implicit=True) [-exp(x)]
* It is possible to solve for anything that can be targeted with subs:
>>> solve(x + 2 + sqrt(3), x + 2) [-sqrt(3)] >>> solve((x + 2 + sqrt(3), x + 4 + y), y, x + 2) {y: -2 + sqrt(3), x + 2: -sqrt(3)}
* Nothing heroic is done in this implicit solving so you may end up with a symbol still in the solution:
>>> eqs = (x*y + 3*y + sqrt(3), x + 4 + y) >>> solve(eqs, y, x + 2) {y: -sqrt(3)/(x + 3), x + 2: (-2*x - 6 + sqrt(3))/(x + 3)} >>> solve(eqs, y*x, x) {x: -y - 4, x*y: -3*y - sqrt(3)}
* if you attempt to solve for a number remember that the number you have obtained does not necessarily mean that the value is equivalent to the expression obtained:
>>> solve(sqrt(2) - 1, 1) [sqrt(2)] >>> solve(x - y + 1, 1) # /!\ -1 is targeted, too [x/(y - 1)] >>> [_.subs(z, -1) for _ in solve((x - y + 1).subs(-1, z), 1)] [-x + y]
* To solve for a function within a derivative, use dsolve.
* single expression and more than 1 symbol
* when there is a linear solution
>>> solve(x - y**2, x, y) [(y**2, y)] >>> solve(x**2 - y, x, y) [(x, x**2)] >>> solve(x**2 - y, x, y, dict=True) [{y: x**2}]
* when undetermined coefficients are identified
* that are linear
>>> solve((a + b)*x - b + 2, a, b) {a: -2, b: 2}
* that are nonlinear
>>> solve((a + b)*x - b**2 + 2, a, b, set=True) ([a, b], {(-sqrt(2), sqrt(2)), (sqrt(2), -sqrt(2))})
* if there is no linear solution then the first successful attempt for a nonlinear solution will be returned
>>> solve(x**2 - y**2, x, y, dict=True) [{x: -y}, {x: y}] >>> solve(x**2 - y**2/exp(x), x, y, dict=True) [{x: 2*LambertW(y/2)}] >>> solve(x**2 - y**2/exp(x), y, x) [(-x*sqrt(exp(x)), x), (x*sqrt(exp(x)), x)]
* iterable of one or more of the above
* involving relationals or bools
>>> solve([x < 3, x - 2]) Eq(x, 2) >>> solve([x > 3, x - 2]) False
* when the system is linear
* with a solution
>>> solve([x - 3], x) {x: 3} >>> solve((x + 5*y - 2, -3*x + 6*y - 15), x, y) {x: -3, y: 1} >>> solve((x + 5*y - 2, -3*x + 6*y - 15), x, y, z) {x: -3, y: 1} >>> solve((x + 5*y - 2, -3*x + 6*y - z), z, x, y) {x: -5*y + 2, z: 21*y - 6}
* without a solution
>>> solve([x + 3, x - 3]) []
* when the system is not linear
>>> solve([x**2 + y -2, y**2 - 4], x, y, set=True) ([x, y], {(-2, -2), (0, 2), (2, -2)})
* if no symbols are given, all free symbols will be selected and a list of mappings returned
>>> solve([x - 2, x**2 + y]) [{x: 2, y: -4}] >>> solve([x - 2, x**2 + f(x)], {f(x), x}) [{x: 2, f(x): -4}]
* if any equation doesn't depend on the symbol(s) given it will be eliminated from the equation set and an answer may be given implicitly in terms of variables that were not of interest
>>> solve([x - y, y - 3], x) {x: y}
Notes =====
solve() with check=True (default) will run through the symbol tags to elimate unwanted solutions. If no assumptions are included all possible solutions will be returned.
>>> from sympy import Symbol, solve >>> x = Symbol("x") >>> solve(x**2 - 1) [-1, 1]
By using the positive tag only one solution will be returned:
>>> pos = Symbol("pos", positive=True) >>> solve(pos**2 - 1) [1]
Assumptions aren't checked when `solve()` input involves relationals or bools.
When the solutions are checked, those that make any denominator zero are automatically excluded. If you do not want to exclude such solutions then use the check=False option:
>>> from sympy import sin, limit >>> solve(sin(x)/x) # 0 is excluded [pi]
If check=False then a solution to the numerator being zero is found: x = 0. In this case, this is a spurious solution since sin(x)/x has the well known limit (without dicontinuity) of 1 at x = 0:
>>> solve(sin(x)/x, check=False) [0, pi]
In the following case, however, the limit exists and is equal to the value of x = 0 that is excluded when check=True:
>>> eq = x**2*(1/x - z**2/x) >>> solve(eq, x) [] >>> solve(eq, x, check=False) [0] >>> limit(eq, x, 0, '-') 0 >>> limit(eq, x, 0, '+') 0
Disabling high-order, explicit solutions ----------------------------------------
When solving polynomial expressions, one might not want explicit solutions (which can be quite long). If the expression is univariate, CRootOf instances will be returned instead:
>>> solve(x**3 - x + 1) [-1/((-1/2 - sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)) - (-1/2 - sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3, -(-1/2 + sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3 - 1/((-1/2 + sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)), -(3*sqrt(69)/2 + 27/2)**(1/3)/3 - 1/(3*sqrt(69)/2 + 27/2)**(1/3)] >>> solve(x**3 - x + 1, cubics=False) [CRootOf(x**3 - x + 1, 0), CRootOf(x**3 - x + 1, 1), CRootOf(x**3 - x + 1, 2)]
If the expression is multivariate, no solution might be returned:
>>> solve(x**3 - x + a, x, cubics=False) []
Sometimes solutions will be obtained even when a flag is False because the expression could be factored. In the following example, the equation can be factored as the product of a linear and a quadratic factor so explicit solutions (which did not require solving a cubic expression) are obtained:
>>> eq = x**3 + 3*x**2 + x - 1 >>> solve(eq, cubics=False) [-1, -1 + sqrt(2), -sqrt(2) - 1]
Solving equations involving radicals ------------------------------------
Because of SymPy's use of the principle root (issue #8789), some solutions to radical equations will be missed unless check=False:
>>> from sympy import root >>> eq = root(x**3 - 3*x**2, 3) + 1 - x >>> solve(eq) [] >>> solve(eq, check=False) [1/3]
In the above example there is only a single solution to the equation. Other expressions will yield spurious roots which must be checked manually; roots which give a negative argument to odd-powered radicals will also need special checking:
>>> from sympy import real_root, S >>> eq = root(x, 3) - root(x, 5) + S(1)/7 >>> solve(eq) # this gives 2 solutions but misses a 3rd [CRootOf(7*_p**5 - 7*_p**3 + 1, 1)**15, CRootOf(7*_p**5 - 7*_p**3 + 1, 2)**15] >>> sol = solve(eq, check=False) >>> [abs(eq.subs(x,i).n(2)) for i in sol] [0.48, 0.e-110, 0.e-110, 0.052, 0.052]
The first solution is negative so real_root must be used to see that it satisfies the expression:
>>> abs(real_root(eq.subs(x, sol[0])).n(2)) 0.e-110
If the roots of the equation are not real then more care will be necessary to find the roots, especially for higher order equations. Consider the following expression:
>>> expr = root(x, 3) - root(x, 5)
We will construct a known value for this expression at x = 3 by selecting the 1-th root for each radical:
>>> expr1 = root(x, 3, 1) - root(x, 5, 1) >>> v = expr1.subs(x, -3)
The solve function is unable to find any exact roots to this equation:
>>> eq = Eq(expr, v); eq1 = Eq(expr1, v) >>> solve(eq, check=False), solve(eq1, check=False) ([], [])
The function unrad, however, can be used to get a form of the equation for which numerical roots can be found:
>>> from sympy.solvers.solvers import unrad >>> from sympy import nroots >>> e, (p, cov) = unrad(eq) >>> pvals = nroots(e) >>> inversion = solve(cov, x)[0] >>> xvals = [inversion.subs(p, i) for i in pvals]
Although eq or eq1 could have been used to find xvals, the solution can only be verified with expr1:
>>> z = expr - v >>> [xi.n(chop=1e-9) for xi in xvals if abs(z.subs(x, xi).n()) < 1e-9] [] >>> z1 = expr1 - v >>> [xi.n(chop=1e-9) for xi in xvals if abs(z1.subs(x, xi).n()) < 1e-9] [-3.0]
See Also ========
- rsolve() for solving recurrence relationships - dsolve() for solving differential equations
""" # keeping track of how f was passed since if it is a list # a dictionary of results will be returned. ###########################################################################
symbols[0] and (isinstance(symbols[0], Symbol) or is_sequence(symbols[0], include=GeneratorType) ) )
# preprocess symbol(s) ########################################################################### # get symbols from equations
# remove symbols the user is not interested in
# preprocess equation(s) ########################################################################### else: else: raise NotImplementedError(filldedent(''' Unanticipated argument of Eq when other arg is True or False. ''')) else:
# rewrite hyperbolics in terms of exp lambda w: w.rewrite(exp))
# if we have a Matrix, we need to iterate over its elements again
# if we can split it into real and imaginary parts then do so # accept as long as new re, im, arg or atan2 are not introduced i.atoms(re, im, arg, atan2) - had for i in (fr, fi)):
# real/imag handling ----------------------------- # Abs continue raise NotImplementedError('solving %s when the argument ' 'is not real or imaginary.' % a) piece(a.args[0]*S.ImaginaryUnit)))
# arg [atan(im(a.args[0])/re(a.args[0])) for a in _arg]))))
# save changes
# see if re(s) or im(s) appear # if re(s) or im(s) appear, the auxiliary equation must be present # end of real/imag handling -----------------------------
# we do this to make the results returned canonical in case f # contains a system of nonlinear equations; all other cases should # be unambiguous
# we can solve for non-symbol entities by replacing them with Dummy symbols
# this is needed in the next two events
# get rid of equations that have no symbols of interest; we don't # try to solve them because the user didn't ask and they might be # hard to solve; this means that solutions may be given in terms # of the eliminated equations e.g. solve((x-y, y-3), x) -> {x: y} # let the solver handle equations that.. # - have no symbols but are expressions # - have symbols of interest # - have no symbols of interest but are constant # but when an expression is not constant and has no symbols of # interest, it can't change what we obtain for a solution from # the remaining equations so we don't include it; and if it's # zero it can be removed and if it's not zero, there is no # solution for the equation set as a whole # # The reason for doing this filtering is to allow an answer # to be obtained to queries like solve((x - y, y), x); without # this mod the return value is [] else: else: ok = True
# mask off any Object that we aren't going to invert: Derivative, # Integral, etc... so that solving for anything that they contain will # give an implicit solution not p.args or p in symset or p.is_Add or p.is_Mul or p.is_Pow and not implicit or p.is_Function and not implicit) and p.func not in (re, im): else:
# Both xreplace and subs are needed below: xreplace to force substitution # inside Derivative, subs to handle non-straightforward substitutions
# rationalize Floats
# Any embedded piecewise functions need to be brought out to the # top level so that the appropriate strategy gets selected. # However, this is necessary only if one of the piecewise # functions depends on one of the symbols we are solving for.
# # try to get a solution ########################################################################### else:
# # postprocessing ########################################################################### # Restore masked-off objects
solution.items()]) in solution] else: else: raise NotImplementedError(filldedent(''' no handling of %s was implemented''' % solution))
# Restore original "symbols" if a dictionary is returned. # This is not necessary for # - the single univariate equation case # since the symbol will have been removed from the solution; # - the nonlinear poly_system since that only supports zero-dimensional # systems and those results come back as a list # # ** unless there were Derivatives with the symbols, but those were handled # above. for k, v in solution.items()]) for k, v in sol.items()])
# undo the dictionary solutions returned when the system was only partially # solved with poly-system if all symbols are present not flags.get('dict', False) and solution and ordered_symbols and not isinstance(solution, dict) and all(isinstance(sol, dict) for sol in solution) ): for r in solution]
# Get assumptions about symbols, to filter solutions. # Note that if assumptions about a solution can't be verified, it is still # returned.
# restore floats
# this has already been checked and is in as_set form got_None.append(sol) else: else: else: # list of expressions
a_None = True else: got_None.append(solution)
elif isinstance(solution, (Relational, And, Or)): if len(symbols) != 1: raise ValueError("Length should be 1") if warn and symbols[0].assumptions0: warnings.warn(filldedent(""" \tWarning: assumptions about variable '%s' are not handled currently.""" % symbols[0])) # TODO: check also variable assumptions for inequalities
else: raise TypeError('Unrecognized solution') # improve the checker
warnings.warn(filldedent(""" \tWarning: assumptions concerning following solution(s) can't be checked:""" + '\n\t' + ', '.join(str(s) for s in got_None)))
# # done ###########################################################################
# Make sure that a list of solutions is ordered in a canonical way.
# return a list of mappings or [] else: else: raise ValueError("Length should be 1") return [], set()
"""Return a checked solution for f in terms of one or more of the symbols. A list should be returned except for the case when a linear undetermined-coefficients equation is encountered (in which case a dictionary is returned).
If no method is implemented to solve the equation, a NotImplementedError will be raised. In the case that conversion of an expression to a Poly gives None a ValueError will be raised."""
# soln may come back as dict, list of dicts or tuples, or # tuple of symbol list and set of solution tuples except NotImplementedError: pass else: raise TypeError('unrecognized args in list') else: raise TypeError('unrecognized solution type') # find first successful solution # no need to check but we should simplify if desired # sol depends on previously solved symbols: discard it # sol depends on previously solved symbols: discard it else: raise NotImplementedError(not_impl_msg % f)
# /!\ capture this flag then set it to False so that no checking in # recursive calls will be done; only the final answer is checked
# build up solutions if f is a Mul # all solutions have been checked but now we must # check that the solutions do not set denominators # in any factor to zero all(not checksol(den, {symbol: s}, **flags) for den in dens)] # set flags for quick exit at end
raise NotImplementedError( 'solve cannot represent interval solutions') # the explicit condition for this expr is the current cond # and none of the previous conditions # an unconditional value was already there # unconditionally take the simplification of v ratio=2, measure=lambda x: 1) # incompatible type with condition(s) (candidate, v), (S.NaN, True))) # the solutions need not be checked again since they # were checked in finding the solutions for each piece else: # first see if it really depends on symbol and whether there # is only a linear solution # no need to check but simplify if desired
# Poly is generally robust enough to convert anything to # a polynomial and tell us the different generators that it # contains, so we will inspect the generators identified by # polys to figure out what to do.
# try to identify a single generator that will allow us to solve this # as a polynomial, followed (perhaps) by a change of variables if the # generator is not a symbol
raise ValueError('could not convert %s to Poly' % f_num) raise ValueError('expression appears to be a constant')
"""Return (b**e, q) for x = b**(p*e/q) where p/q is the leading Rational of the exponent of x, e.g. exp(-2*x/3) -> (exp(x), 3) """
# If there is more than one generator, it could be that the # generators have the same base but different powers, e.g. # >>> Poly(exp(x) + 1/exp(x)) # Poly(exp(-x) + exp(x), exp(-x), exp(x), domain='ZZ') # # If unrad was not disabled then there should be no rational # exponents appearing as in # >>> Poly(sqrt(x) + sqrt(sqrt(x))) # Poly(sqrt(x) + x**(1/4), sqrt(x), x**(1/4), domain='ZZ')
isinstance(_, TrigonometricFunction)])
# just a simple case - see if replacement of single function # clears all symbol-dependent functions, e.g. # log(x) - log(log(x) - 1) - 3 can be solved even though it has # two generators.
# perform the substitution
# if no Functions left, we can proceed with usual solve
else: # e.g. case where gens are exp(x), exp(-x) # this will be resolved by factor in _tsolve but we might # as well try a simple expansion here to get things in # order so something like the following will work now without # having to factor: # # >>> eq = (exp(I*(-x-2))+exp(I*(x+2))) # >>> eq.subs(exp(x),y) # fails # exp(I*(-x - 2)) + exp(I*(x + 2)) # >>> eq.expand().subs(exp(x),y) # works # y**I*exp(2*I) + y**(-I)*exp(-2*I) lambda w: w.is_Pow or isinstance(w, exp), _expand).subs(u, t)
# There is only one generator that we are interested in, but # there may have been more than one generator identified by # polys (e.g. for symbols other than the one we are interested # in) so recast the poly in terms of our generator of interest. # Also use composite=True with f_num since Poly won't update # poly as documented in issue 8810.
# if we aren't on the tsolve-pass, use roots ('cubics', 'quartics', 'quintics')]) # e.g. roots(32*x**5 + 400*x**4 + 2032*x**3 + # 5000*x**2 + 6250*x + 3189) -> {} # so all_roots is used and RootOf instances are # returned *unless* the system is multivariate # or high-order EX domain. raise NotImplementedError( filldedent(''' Neither high-order multivariate polynomials nor sorting of EX-domain polynomials is supported. If you want to see any results, pass keyword incomplete=True to solve; to see numerical values of roots for univariate expressions, use nroots. ''')) else: pass else:
# perhaps _tsolve can handle f_num else: # if the flag wasn't set then unset it since high-order # results are quite long. Perhaps one could base this # decision on a certain critical length of the # roots. In addition, wester test M2 has an expression # whose roots can be shown to be real with the # unsimplified form of the solution whereas only one of # the simplified forms appears to be real.
# fallback if above fails # ----------------------- # try unrad except (ValueError, NotImplementedError): u = False else: except NotImplementedError: rv = None # if the flag wasn't set then unset it since unrad results # can be quite long or of very high order else:
# try _tsolve # ----------- end of fallback ----------------------------
raise NotImplementedError('\n'.join([msg, not_impl_msg % f]))
# we just simplified the solution so we now set the flag to # False so the simplification doesn't happen again in checksol()
# reject any result that makes any denom. affirmatively 0; # if in doubt, keep it all(not checksol(d, {symbol: s}, **flags) for d in dens)] # keep only results if the check is not False checksol(f_num, {symbol: r}, **flags) is not False]
return []
else:
else:
# returns a dictionary ({symbols: values}) or None else: else: else:
else:
# returns [] or list of tuples of solutions for syms for ss in got_s]): # sol depends on previously # solved symbols: discard it except NotImplementedError: pass else: raise NotImplementedError('no valid subset found') else: # we don't know here if the symbols provided were given # or not, so let solve resolve that. A list of dictionaries # is going to always be returned from here. #
else:
# For each failed equation, see if we can solve for one of the # remaining symbols from that equation. If so, we update the # solution set and continue with the next failed equation, # repeating until we are done or we get an equation that can't # be solved.
# sort so equation with the fewest potential symbols is first # update eq with everything that is known so far # if check is True then we see if it satisfies this # equation, otherwise we just accept it # this solution is sufficient to know whether # it is valid or not so we either accept or # reject it, then continue else: # search for a symbol amongst those available that # can be solved for except NotImplementedError: continue # put each solution in r and append the now-expanded # result in the new result list; use copy since the # solution for s in being added in-place # sol depends on previously solved symbols: discard it # and add this new solution raise NotImplementedError('could not solve %s' % eq2) else: result.remove(b)
if not any(checksol(d, r, **flags) for d in dens)]
if not any(checksol(e, r, **flags) is False for e in exprs)]
r""" Return a tuple derived from f = lhs - rhs that is one of the following:
(0, 1) meaning that ``f`` is independent of the symbols in ``symbols`` that aren't in ``exclude``, e.g::
>>> from sympy.solvers.solvers import solve_linear >>> from sympy.abc import x, y, z >>> from sympy import cos, sin >>> eq = y*cos(x)**2 + y*sin(x)**2 - y # = y*(1 - 1) = 0 >>> solve_linear(eq) (0, 1) >>> eq = cos(x)**2 + sin(x)**2 # = 1 >>> solve_linear(eq) (0, 1) >>> solve_linear(x, exclude=[x]) (0, 1)
(0, 0) meaning that there is no solution to the equation amongst the symbols given.
(If the first element of the tuple is not zero then the function is guaranteed to be dependent on a symbol in ``symbols``.)
(symbol, solution) where symbol appears linearly in the numerator of ``f``, is in ``symbols`` (if given) and is not in ``exclude`` (if given). No simplification is done to ``f`` other than a ``mul=True`` expansion, so the solution will correspond strictly to a unique solution.
``(n, d)`` where ``n`` and ``d`` are the numerator and denominator of ``f`` when the numerator was not linear in any symbol of interest; ``n`` will never be a symbol unless a solution for that symbol was found (in which case the second element is the solution, not the denominator).
Examples ========
>>> from sympy.core.power import Pow >>> from sympy.polys.polytools import cancel
The variable ``x`` appears as a linear variable in each of the following:
>>> solve_linear(x + y**2) (x, -y**2) >>> solve_linear(1/x - y**2) (x, y**(-2))
When not linear in x or y then the numerator and denominator are returned.
>>> solve_linear(x**2/y**2 - 3) (x**2 - 3*y**2, y**2)
If the numerator of the expression is a symbol then (0, 0) is returned if the solution for that symbol would have set any denominator to 0:
>>> eq = 1/(1/x - 2) >>> eq.as_numer_denom() (x, -2*x + 1) >>> solve_linear(eq) (0, 0)
But automatic rewriting may cause a symbol in the denominator to appear in the numerator so a solution will be returned:
>>> (1/x)**-1 x >>> solve_linear((1/x)**-1) (x, 0)
Use an unevaluated expression to avoid this:
>>> solve_linear(Pow(1/x, -1, evaluate=False)) (0, 0)
If ``x`` is allowed to cancel in the following expression, then it appears to be linear in ``x``, but this sort of cancellation is not done by ``solve_linear`` so the solution will always satisfy the original expression without causing a division by zero error.
>>> eq = x**2*(1/x - z**2/x) >>> solve_linear(cancel(eq)) (x, 0) >>> solve_linear(eq) (x**2*(-z**2 + 1), x)
A list of symbols for which a solution is desired may be given:
>>> solve_linear(x + y + z, symbols=[y]) (y, -x - z)
A list of symbols to ignore may also be given:
>>> solve_linear(x + y + z, exclude=[x]) (y, -x - z)
(A solution for ``y`` is obtained because it is the first variable from the canonically sorted list of symbols that had a linear solution.)
""" If lhs is an Equality, rhs must be 0 but was %s''' % rhs))
else: else: eg = 'solve(%s, *%s)' % (eq, list(symbols)) solve_linear only handles symbols, not %s. To isolate non-symbols use solve, e.g. >>> %s <<<. ''' % (bad, eg)))
# derivatives are easy to do but tricky to analyze to see if they # are going to disallow a linear solution, so for simplicity we # just evaluate the ones that have the symbols of interest
# if there are derivatives in this var, calculate them now # dnewn_dxi can be nonzero if it survives differentation by any # of its free symbols for di in dens): # simplify any trivial integral i.function.is_number] # do a slight bit of simplification
r""" Find a particular solution to a linear system.
In particular, try to find a solution with the minimal possible number of non-zero variables using a naive algorithm with exponential complexity. If ``quick=True``, a heuristic is used. """ # Check if there are any non-zero solutions at all return s0 # We just solve the system and try to heuristically find a nice # solution. # NOTE sort by default_sort_key to get deterministic result key=lambda x: (len(x.free_symbols), default_sort_key(x))) else: else: determined[x] = val else: # We try to select n variables which we want to be non-zero. # All others will be assumed zero. We try to solve the modified system. # If there is a non-trivial solution, just set the free variables to # one. If we do this for increasing n, trying all combinations of # variables, we will find an optimal solution. # We speed up slightly by starting at one less than the number of # variables the quick method manages. else: break
r""" Solve system of N linear equations with M variables, which means both under- and overdetermined systems are supported. The possible number of solutions is zero, one or infinite. Respectively, this procedure will return None or a dictionary with solutions. In the case of underdetermined systems, all arbitrary parameters are skipped. This may cause a situation in which an empty dictionary is returned. In that case, all symbols can be assigned arbitrary values.
Input to this functions is a Nx(M+1) matrix, which means it has to be in augmented form. If you prefer to enter N equations and M unknowns then use `solve(Neqs, *Msymbols)` instead. Note: a local copy of the matrix is made by this routine so the matrix that is passed will not be modified.
The algorithm used here is fraction-free Gaussian elimination, which results, after elimination, in an upper-triangular matrix. Then solutions are found using back-substitution. This approach is more efficient and compact than the Gauss-Jordan method.
>>> from sympy import Matrix, solve_linear_system >>> from sympy.abc import x, y
Solve the following system::
x + 4 y == 2 -2 x + y == 14
>>> system = Matrix(( (1, 4, 2), (-2, 1, 14))) >>> solve_linear_system(system, x, y) {x: -6, y: 2}
A degenerate system returns an empty dictionary.
>>> system = Matrix(( (0,0,0), (0,0,0) )) >>> solve_linear_system(system, x, y) {}
"""
# well behaved n-equations and n-unknowns # non-trivial solution
# an overdetermined system else: # remove trailing rows
# there is no pivot in current column # so try to find one in other columns else: # We need to know if this is always zero or not. We # assume that if there are free symbols that it is not # identically zero (or that there is more than one way # to make this zero). Otherwise, if there are none, this # is a constant and we assume that it does not simplify # to zero XXX are there better (fast) ways to test this? # The .equals(0) method could be used but that can be # slow; numerical testing is prone to errors of scaling.
# A row of zeros with a non-zero rhs can only be accepted # if there is another equivalent row. Any such rows will # be deleted. # do we need to see if the rhs of j # is a constant multiple of i's rhs? if ip is None: _, ip = rowi[-1].as_content_primitive() _, jp = rowj[-1].as_content_primitive() if not (simplify(jp - ip) or simplify(jp + ip)): matrix.row_del(j)
# no solution return None # zero row or was a linear combination of # other rows or was a row with a symbolic # expression that matched other rows, e.g. [0, 0, x - y] # so now we can safely skip it # every choice of variable values is a solution # so we return an empty dict instead of None return dict() continue
# we want to change the order of columns so # the order of variables must also change
# divide all elements in the current row by the pivot
# subtract from the current row the row containing # pivot and multiplied by extracted coefficient
# if there weren't any problems, augmented matrix is now # in row-echelon form so we can check how many solutions # there are and extract them using back substitution
# this system is Cramer equivalent so there is # exactly one solution to this system of equations
# run back-substitution for variables
else: solutions[syms[k]] = content
# this system will have infinite number of solutions # dependent on exactly len(syms) - i parameters
# run back-substitution for variables
# run back-substitution for parameters
else: solutions[syms[k]] = content
else: return [] # no solutions
"""Solve equation of a type p(x; a_1, ..., a_k) == q(x) where both p, q are univariate polynomials and f depends on k parameters. The result of this functions is a dictionary with symbolic values of those parameters with respect to coefficients in q.
This functions accepts both Equations class instances and ordinary SymPy expressions. Specification of parameters and variable is obligatory for efficiency and simplicity reason.
>>> from sympy import Eq >>> from sympy.abc import a, b, c, x >>> from sympy.solvers import solve_undetermined_coeffs
>>> solve_undetermined_coeffs(Eq(2*a*x + a+b, x), [a, b], x) {a: 1/2, b: -1/2}
>>> solve_undetermined_coeffs(Eq(a*c*x + a+b, x), [a, b], x) {a: 1/c, b: -1/c}
""" # got equation, so move all the # terms to the left hand side equ = equ.lhs - equ.rhs
# consecutive powers in the input expressions have # been successfully collected, so solve remaining # system using Gaussian elimination algorithm else: return None # no solutions
""" Solves the augmented matrix system using LUsolve and returns a dictionary in which solutions are keyed to the symbols of syms *as ordered*.
The matrix must be invertible.
Examples ========
>>> from sympy import Matrix >>> from sympy.abc import x, y, z >>> from sympy.solvers.solvers import solve_linear_system_LU
>>> solve_linear_system_LU(Matrix([ ... [1, 2, 0, 1], ... [3, 2, 2, 1], ... [2, 0, 0, 1]]), [x, y, z]) {x: 1/2, y: 1/4, z: -1/2}
See Also ========
sympy.matrices.LUsolve
""" raise ValueError("Rows should be equal to columns - 1")
"""Return the det(``M``) by using permutations to select factors. For size larger than 8 the number of permutations becomes prohibitively large, or if there are no symbols in the matrix, it is better to use the standard determinant routines, e.g. `M.det()`.
See Also ======== det_minor det_quick """ except AttributeError: list = flatten(M.tolist())
"""Return the ``det(M)`` computed from minors without introducing new nesting in products.
See Also ======== det_perm det_quick """ else: Add.make_args(det_minor(M.minor_submatrix(0, i)))]) if M[0, i] else S.Zero for i in range(n)])
"""Return ``det(M)`` assuming that either there are lots of zeros or the size of the matrix is small. If this assumption is not met, then the normal Matrix.det function will be used with method = ``method``.
See Also ======== det_minor det_perm """ else:
"""Return the inverse of ``M``, assuming that either there are lots of zeros or the size of the matrix is small. """ else: else: raise ValueError("Matrix det == 0; not invertible.")
# these are functions that have multiple inverse values per period sin: lambda x: (asin(x), S.Pi - asin(x)), cos: lambda x: (acos(x), 2*S.Pi - acos(x)), }
""" Helper for _solve that solves a transcendental equation with respect to the given symbol. Various equations containing powers and logarithms, can be solved.
There is currently no guarantee that all solutions will be returned or that a real solution will be favored over a complex one.
Either a list of potential solutions will be returned or None will be returned (in the case that no method was known to get a solution for the equation). All other errors (like the inability to cast an expression as a Poly) are unhandled.
Examples ========
>>> from sympy import log >>> from sympy.solvers.solvers import _tsolve as tsolve >>> from sympy.abc import x
>>> tsolve(3**(2*x + 5) - 4, x) [-5/2 + log(2)/log(3), (-5*log(3)/2 + log(2) + I*pi)/log(3)]
>>> tsolve(log(x) + 2*x, x) [LambertW(2)/2]
""" else:
# it's time to try factoring; powdenest is used # to try get powers in standard form for better factoring
# f(x)**g(x) only has solutions where f(x) == 0 and g(x) != 0 at # the same place return sol_base # no solutions to remove so return now _solve(lhs.exp, sym, **flags)))) lhs.base.is_positive and lhs.exp.is_real):
# sin(x) = 1/3 -> x - asin(1/3) & x - (pi - asin(1/3))
# maybe it is a lambert pattern # lambert forms may need some help being recognized, e.g. changing # 2**(3*x) + x**3*log(2)**3 + 3*x**2*log(2)**2 + 3*x*log(2) + 1 # to 2**(3*x) + (x*log(2) + 1)**3 dict(list(zip(up_or_log, [0]*len(up_or_log))))) # maybe it's a convoluted function raise NotImplementedError for i in inversion for s in sol]))) else:
else: else:
# TODO: option for calculating J numerically
def nsolve(*args, **kwargs): r""" Solve a nonlinear equation system numerically::
nsolve(f, [args,] x0, modules=['mpmath'], **kwargs)
f is a vector function of symbolic expressions representing the system. args are the variables. If there is only one variable, this argument can be omitted. x0 is a starting vector close to a solution.
Use the modules keyword to specify which modules should be used to evaluate the function and the Jacobian matrix. Make sure to use a module that supports matrices. For more information on the syntax, please see the docstring of lambdify.
If the keyword arguments contain 'dict'=True (default is False) nsolve will return a list (perhaps empty) of solution mappings. This might be especially useful if you want to use nsolve as a fallback to solve since using the dict argument for both methods produces return values of consistent type structure. Please note: to keep this consistency with solve, the solution will be returned in a list even though nsolve (currently at least) only finds one solution at a time.
Overdetermined systems are supported.
>>> from sympy import Symbol, nsolve >>> import sympy >>> import mpmath >>> mpmath.mp.dps = 15 >>> x1 = Symbol('x1') >>> x2 = Symbol('x2') >>> f1 = 3 * x1**2 - 2 * x2**2 - 1 >>> f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8 >>> print(nsolve((f1, f2), (x1, x2), (-1, 1))) Matrix([[-1.19287309935246], [1.27844411169911]])
For one-dimensional functions the syntax is simplified:
>>> from sympy import sin, nsolve >>> from sympy.abc import x >>> nsolve(sin(x), x, 2) 3.14159265358979 >>> nsolve(sin(x), 2) 3.14159265358979
To solve with higher precision than the default, use the prec argument.
>>> from sympy import cos >>> nsolve(cos(x) - x, 1) 0.739085133215161 >>> nsolve(cos(x) - x, 1, prec=50) 0.73908513321516064165531208767387340401341175890076 >>> cos(_) 0.73908513321516064165531208767387340401341175890076
To solve for complex roots of real functions, a nonreal initial point must be specified:
>>> from sympy import I >>> nsolve(x**2 + 2, I) 1.4142135623731*I
mpmath.findroot is used and you can find there more extensive documentation, especially concerning keyword parameters and available solvers. Note, however, that functions which are very steep near the root the verification of the solution may fail. In this case you should use the flag `verify=False` and independently verify the solution.
>>> from sympy import cos, cosh >>> from sympy.abc import i >>> f = cos(x)*cosh(x) - 1 >>> nsolve(f, 3.14*100) Traceback (most recent call last): ... ValueError: Could not find root within given tolerance. (1.39267e+230 > 2.1684e-19) >>> ans = nsolve(f, 3.14*100, verify=False); ans 312.588469032184 >>> f.subs(x, ans).n(2) 2.1e+121 >>> (f/f.diff(x)).subs(x, ans).n(2) 7.4e-15
One might safely skip the verification if bounds of the root are known and a bisection method is used:
>>> bounds = lambda i: (3.14*i, 3.14*(i + 1)) >>> nsolve(f, bounds(100), solver='bisect', verify=False) 315.730061685774
Alternatively, a function may be better behaved when the denominator is ignored. Since this is not always the case, however, the decision of what function to use is left to the discretion of the user.
>>> eq = x**2/(1 - x)/(1 - 2*x)**2 - 100 >>> nsolve(eq, 0.46) Traceback (most recent call last): ... ValueError: Could not find root within given tolerance. (10000 > 2.1684e-19) Try another starting point or tweak arguments. >>> nsolve(eq.as_numer_denom()[0], 0.46) 0.46792545969349058 """ # there are several other SymPy functions that use method= so # guard against that here Keyword "method" should not be used in this context. When using some mpmath solvers directly, the keyword "method" is used, but when using nsolve (and findroot) the keyword to use is "solver".'''))
else:
# keyword argument to return result as a dictionary
# interpret arguments % (len(fargs), len(x0))) % len(args)) else: % len(args)) # assume it's a sympy expression raise ValueError(filldedent(''' expected a one-dimensional and numerical function'''))
# the function is much better behaved if there is no denominator # but sending the numerator is left to the user since sometimes # the function is better behaved when the denominator is present # e.g., issue 11768
raise NotImplementedError(filldedent(''' need at least as many equations as variables''')) print('f(x):') print(f) # derive Jacobian print('J(x):') print(J) # create functions # solve the system numerically
"""Return tuple (i, d) where ``i`` is independent of ``symbols`` and ``d`` contains symbols. ``i`` and ``d`` are obtained after recursively using algebraic inversion until an uninvertible ``d`` remains. If there are no free symbols then ``d`` will be zero. Some (but not necessarily all) solutions to the expression ``i - d`` will be related to the solutions of the original expression.
Examples ========
>>> from sympy.solvers.solvers import _invert as invert >>> from sympy import sqrt, cos >>> from sympy.abc import x, y >>> invert(x - 3) (3, x) >>> invert(3) (3, 0) >>> invert(2*cos(x) - 1) (1/2, cos(x)) >>> invert(sqrt(x) - 3) (3, sqrt(x)) >>> invert(sqrt(x) + y, x) (-y, sqrt(x)) >>> invert(sqrt(x) + y, y) (-sqrt(x), y) >>> invert(sqrt(x) + y, x, y) (0, sqrt(x) + y)
If there is more than one symbol in a power's base and the exponent is not an Integer, then the principal root will be used for the inversion:
>>> invert(sqrt(x + y) - 2) (4, x + y) >>> invert(sqrt(x + y) - 2) (4, x + y)
If the exponent is an integer, setting ``integer_power`` to True will force the principal root to be selected:
>>> invert(x**2 - 4, integer_power=True) (2, x)
"""
# dep + indep == rhs # this indicates we have done it all
# dep * indep == rhs else: # this indicates we have done it all
# collect like-terms in symbols else:
# if it's a two-term Add with rhs = 0 and two powers we can get the # dependent terms together, e.g. 3*f(x) + 2*g(x) -> f(x)/g(x) = -2/3 not lhs.is_polynomial(*symbols): # a = -b else: isinstance(i, Function) for i in (ad, bd)) and \ ad.func == bd.func and len(ad.args) == len(bd.args): else: # should be able to solve # f(x, y) == f(2, 3) -> x == 2 # f(x, x + y) == f(2, 3) -> x == 2 or x == 3 - y raise NotImplementedError('equal function with more than 1 argument')
# -1 # f(x) = g -> x = f (g) # # /!\ inverse should not be defined if there are multiple values # for the function -- these are handled in _tsolve #
# base**a = b -> base = b**(1/a) if # a is an Integer and dointpow=True (this gives real branch of root) # a is not an Integer and the equation is multivariate and the # base has more than 1 symbol in it # The rationale for this is that right now the multi-system solvers # doesn't try to resolve generators to see, for example, if the whole # system is written in terms of sqrt(x + y) so it will just fail, so we # do that step here. lhs.exp.is_Integer and dointpow or not lhs.exp.is_Integer and len(symbols) > 1 and len(lhs.base.free_symbols & set(symbols)) > 1):
""" Remove radicals with symbolic arguments and return (eq, cov), None or raise an error:
None is returned if there are no radicals to remove.
NotImplementedError is raised if there are radicals and they cannot be removed or if the relationship between the original symbols and the change of variable needed to rewrite the system as a polynomial cannot be solved.
Otherwise the tuple, ``(eq, cov)``, is returned where::
``eq``, ``cov`` ``eq`` is an equation without radicals (in the symbol(s) of interest) whose solutions are a superset of the solutions to the original expression. ``eq`` might be re-written in terms of a new variable; the relationship to the original variables is given by ``cov`` which is a list containing ``v`` and ``v**p - b`` where ``p`` is the power needed to clear the radical and ``b`` is the radical now expressed as a polynomial in the symbols of interest. For example, for sqrt(2 - x) the tuple would be ``(c, c**2 - 2 + x)``. The solutions of ``eq`` will contain solutions to the original equation (if there are any).
``syms`` an iterable of symbols which, if provided, will limit the focus of radical removal: only radicals with one or more of the symbols of interest will be cleared. All free symbols are used if ``syms`` is not set.
``flags`` are used internally for communication during recursive calls. Two options are also recognized::
``take``, when defined, is interpreted as a single-argument function that returns True if a given Pow should be handled.
Radicals can be removed from an expression if::
* all bases of the radicals are the same; a change of variables is done in this case. * if all radicals appear in one term of the expression * there are only 4 terms with sqrt() factors or there are less than four terms having sqrt() factors * there are only two terms with radicals
Examples ========
>>> from sympy.solvers.solvers import unrad >>> from sympy.abc import x >>> from sympy import sqrt, Rational, root, real_roots, solve
>>> unrad(sqrt(x)*x**Rational(1, 3) + 2) (x**5 - 64, []) >>> unrad(sqrt(x) + root(x + 1, 3)) (x**3 - x**2 - 2*x - 1, []) >>> eq = sqrt(x) + root(x, 3) - 2 >>> unrad(eq) (_p**3 + _p**2 - 2, [_p, _p**6 - x])
"""
# XXX - uncovered oldp, olde = cov if Poly(e, p).degree(p) in (1, 2): cov[:] = [p, olde.subs(oldp, _solve(e, p, **uflags)[0])] else: raise NotImplementedError else:
# change symbol to vanilla so no solutions are eliminated
# remove constants and powers of factors since these don't change # the location of the root; XXX should factor or factor_terms be used? else:
# make the sign canonical
# return leading Rational of denominator of Pow's exponent return S.One
# define the _take method that will determine whether a term is of interest # return True if coefficient of any factor's exponent's den is not 1
sorted(dict(cov=[], n=None, rpt=0).items())]
# preconditioning
# check for trivial case # - already a polynomial in integer powers # - an exponent has a symbol of interest (don't handle) return
# if all the bases are the same or all the radicals are in one # term, `lcm` will be the lcm of the denominators of the # exponents of the radicals
return
# only keep in syms symbols that actually appear in radicals; # and update gens
# get terms together that have common generators else:
# the output will depend on the order terms are processed, so # make it canonical quickly
else: rterms = list(rterms[0].args) free = b.free_symbols x = {g for g in gens if g.is_Symbol} & free if not x: x = free x = ordered(x) else: raise NotImplementedError except NotImplementedError: pass else: # no longer consider integer powers as generators
# the lcm-is-power-of-two case is handled below r0, r1 = rterms if flags.get('_reverse', False): r1, r0 = r0, r1 i0 = _rads0, _bases0, lcm0 = _rads_bases_lcm(r0.as_poly()) i1 = _rads1, _bases1, lcm1 = _rads_bases_lcm(r1.as_poly()) for reverse in range(2): if reverse: i0, i1 = i1, i0 r0, r1 = r1, r0 _rads1, _, lcm1 = i1 _rads1 = Mul(*_rads1) t1 = _rads1**lcm1 c = covsym**lcm1 - t1 for x in syms: try: sol = _solve(c, x, **uflags) if not sol: raise NotImplementedError neweq = r0.subs(x, sol[0]) + covsym*r1/_rads1 + \ others tmp = unrad(neweq, covsym) if tmp: eq, newcov = tmp if newcov: newp, newc = newcov _cov(newp, c.subs(covsym, _solve(newc, covsym, **uflags)[0])) else: _cov(covsym, c) else: eq = neweq _cov(covsym, c) ok = True break except NotImplementedError: if reverse: raise NotImplementedError( 'no successful change of variable found') else: pass if ok: break # two cube roots and another with order less than 5 # (so an analytical solution can be found) or a base # that matches one of the cube root bases elif info[1][LCM] != 3: info.append(info.pop(1)) rterms.append(rterms.pop(1)) if info[1][BASES] != info[2][BASES]: info[0], info[1] = info[1], info[0] rterms[0], rterms[1] = rterms[1], rterms[0] if info[1][BASES] == info[2][BASES]: eq = rterms[0]**3 + (rterms[1] + rterms[2] + others)**3 ok = True elif info[2][LCM] < 5: # a*root(A, 3) + b*root(B, 3) + others = c a, b, c, d, A, B = [Dummy(i) for i in 'abcdAB'] # zz represents the unraded expression into which the # specifics for this case are substituted zz = (c - d)*(A**3*a**9 + 3*A**2*B*a**6*b**3 - 3*A**2*a**6*c**3 + 9*A**2*a**6*c**2*d - 9*A**2*a**6*c*d**2 + 3*A**2*a**6*d**3 + 3*A*B**2*a**3*b**6 + 21*A*B*a**3*b**3*c**3 - 63*A*B*a**3*b**3*c**2*d + 63*A*B*a**3*b**3*c*d**2 - 21*A*B*a**3*b**3*d**3 + 3*A*a**3*c**6 - 18*A*a**3*c**5*d + 45*A*a**3*c**4*d**2 - 60*A*a**3*c**3*d**3 + 45*A*a**3*c**2*d**4 - 18*A*a**3*c*d**5 + 3*A*a**3*d**6 + B**3*b**9 - 3*B**2*b**6*c**3 + 9*B**2*b**6*c**2*d - 9*B**2*b**6*c*d**2 + 3*B**2*b**6*d**3 + 3*B*b**3*c**6 - 18*B*b**3*c**5*d + 45*B*b**3*c**4*d**2 - 60*B*b**3*c**3*d**3 + 45*B*b**3*c**2*d**4 - 18*B*b**3*c*d**5 + 3*B*b**3*d**6 - c**9 + 9*c**8*d - 36*c**7*d**2 + 84*c**6*d**3 - 126*c**5*d**4 + 126*c**4*d**5 - 84*c**3*d**6 + 36*c**2*d**7 - 9*c*d**8 + d**9) def _t(i): b = Mul(*info[i][RAD]) return cancel(rterms[i]/b), Mul(*info[i][BASES]) aa, AA = _t(0) bb, BB = _t(1) cc = -rterms[2] dd = others eq = zz.xreplace(dict(zip( (a, A, b, B, c, d), (aa, AA, bb, BB, cc, dd)))) ok = True # handle power-of-2 cases len(rterms) == 4 or len(rterms) < 4):
# (r0+r1)**2 - (r2+r3)**2 r0, r1, r2, r3 = rterms eq = _norm2(r0, r1) - _norm2(r2, r3) ok = True # (r1+r2)**2 - (r0+others)**2 # r0**2 - (r1+others)**2
nwas is not None and len(rterms) == nwas and new_depth is not None and new_depth == depth and rpt > 3): raise NotImplementedError('Cannot remove all radicals')
bivariate_type, _solve_lambert, _filtered_gens) |