Coverage for /usr/lib/python3/dist-packages/sympy/solvers/solvers.py: 4%
1622 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1"""
2This module contain solvers for all kinds of equations:
4 - algebraic or transcendental, use solve()
6 - recurrence, use rsolve()
8 - differential, use dsolve()
10 - nonlinear (numerically), use nsolve()
11 (you will need a good starting point)
13"""
14from __future__ import annotations
16from sympy.core import (S, Add, Symbol, Dummy, Expr, Mul)
17from sympy.core.assumptions import check_assumptions
18from sympy.core.exprtools import factor_terms
19from sympy.core.function import (expand_mul, expand_log, Derivative,
20 AppliedUndef, UndefinedFunction, nfloat,
21 Function, expand_power_exp, _mexpand, expand,
22 expand_func)
23from sympy.core.logic import fuzzy_not
24from sympy.core.numbers import ilcm, Float, Rational, _illegal
25from sympy.core.power import integer_log, Pow
26from sympy.core.relational import Eq, Ne
27from sympy.core.sorting import ordered, default_sort_key
28from sympy.core.sympify import sympify, _sympify
29from sympy.core.traversal import preorder_traversal
30from sympy.logic.boolalg import And, BooleanAtom
32from sympy.functions import (log, exp, LambertW, cos, sin, tan, acos, asin, atan,
33 Abs, re, im, arg, sqrt, atan2)
34from sympy.functions.combinatorial.factorials import binomial
35from sympy.functions.elementary.hyperbolic import HyperbolicFunction
36from sympy.functions.elementary.piecewise import piecewise_fold, Piecewise
37from sympy.functions.elementary.trigonometric import TrigonometricFunction
38from sympy.integrals.integrals import Integral
39from sympy.ntheory.factor_ import divisors
40from sympy.simplify import (simplify, collect, powsimp, posify, # type: ignore
41 powdenest, nsimplify, denom, logcombine, sqrtdenest, fraction,
42 separatevars)
43from sympy.simplify.sqrtdenest import sqrt_depth
44from sympy.simplify.fu import TR1, TR2i
45from sympy.matrices.common import NonInvertibleMatrixError
46from sympy.matrices import Matrix, zeros
47from sympy.polys import roots, cancel, factor, Poly
48from sympy.polys.polyerrors import GeneratorsNeeded, PolynomialError
49from sympy.polys.solvers import sympy_eqs_to_ring, solve_lin_sys
50from sympy.utilities.lambdify import lambdify
51from sympy.utilities.misc import filldedent, debugf
52from sympy.utilities.iterables import (connected_components,
53 generate_bell, uniq, iterable, is_sequence, subsets, flatten)
54from sympy.utilities.decorator import conserve_mpmath_dps
56from mpmath import findroot
58from sympy.solvers.polysys import solve_poly_system
60from types import GeneratorType
61from collections import defaultdict
62from itertools import combinations, product
64import warnings
67def recast_to_symbols(eqs, symbols):
68 """
69 Return (e, s, d) where e and s are versions of *eqs* and
70 *symbols* in which any non-Symbol objects in *symbols* have
71 been replaced with generic Dummy symbols and d is a dictionary
72 that can be used to restore the original expressions.
74 Examples
75 ========
77 >>> from sympy.solvers.solvers import recast_to_symbols
78 >>> from sympy import symbols, Function
79 >>> x, y = symbols('x y')
80 >>> fx = Function('f')(x)
81 >>> eqs, syms = [fx + 1, x, y], [fx, y]
82 >>> e, s, d = recast_to_symbols(eqs, syms); (e, s, d)
83 ([_X0 + 1, x, y], [_X0, y], {_X0: f(x)})
85 The original equations and symbols can be restored using d:
87 >>> assert [i.xreplace(d) for i in eqs] == eqs
88 >>> assert [d.get(i, i) for i in s] == syms
90 """
91 if not iterable(eqs) and iterable(symbols):
92 raise ValueError('Both eqs and symbols must be iterable')
93 orig = list(symbols)
94 symbols = list(ordered(symbols))
95 swap_sym = {}
96 i = 0
97 for j, s in enumerate(symbols):
98 if not isinstance(s, Symbol) and s not in swap_sym:
99 swap_sym[s] = Dummy('X%d' % i)
100 i += 1
101 new_f = []
102 for i in eqs:
103 isubs = getattr(i, 'subs', None)
104 if isubs is not None:
105 new_f.append(isubs(swap_sym))
106 else:
107 new_f.append(i)
108 restore = {v: k for k, v in swap_sym.items()}
109 return new_f, [swap_sym.get(i, i) for i in orig], restore
112def _ispow(e):
113 """Return True if e is a Pow or is exp."""
114 return isinstance(e, Expr) and (e.is_Pow or isinstance(e, exp))
117def _simple_dens(f, symbols):
118 # when checking if a denominator is zero, we can just check the
119 # base of powers with nonzero exponents since if the base is zero
120 # the power will be zero, too. To keep it simple and fast, we
121 # limit simplification to exponents that are Numbers
122 dens = set()
123 for d in denoms(f, symbols):
124 if d.is_Pow and d.exp.is_Number:
125 if d.exp.is_zero:
126 continue # foo**0 is never 0
127 d = d.base
128 dens.add(d)
129 return dens
132def denoms(eq, *symbols):
133 """
134 Return (recursively) set of all denominators that appear in *eq*
135 that contain any symbol in *symbols*; if *symbols* are not
136 provided then all denominators will be returned.
138 Examples
139 ========
141 >>> from sympy.solvers.solvers import denoms
142 >>> from sympy.abc import x, y, z
144 >>> denoms(x/y)
145 {y}
147 >>> denoms(x/(y*z))
148 {y, z}
150 >>> denoms(3/x + y/z)
151 {x, z}
153 >>> denoms(x/2 + y/z)
154 {2, z}
156 If *symbols* are provided then only denominators containing
157 those symbols will be returned:
159 >>> denoms(1/x + 1/y + 1/z, y, z)
160 {y, z}
162 """
164 pot = preorder_traversal(eq)
165 dens = set()
166 for p in pot:
167 # Here p might be Tuple or Relational
168 # Expr subtrees (e.g. lhs and rhs) will be traversed after by pot
169 if not isinstance(p, Expr):
170 continue
171 den = denom(p)
172 if den is S.One:
173 continue
174 for d in Mul.make_args(den):
175 dens.add(d)
176 if not symbols:
177 return dens
178 elif len(symbols) == 1:
179 if iterable(symbols[0]):
180 symbols = symbols[0]
181 return {d for d in dens if any(s in d.free_symbols for s in symbols)}
184def checksol(f, symbol, sol=None, **flags):
185 """
186 Checks whether sol is a solution of equation f == 0.
188 Explanation
189 ===========
191 Input can be either a single symbol and corresponding value
192 or a dictionary of symbols and values. When given as a dictionary
193 and flag ``simplify=True``, the values in the dictionary will be
194 simplified. *f* can be a single equation or an iterable of equations.
195 A solution must satisfy all equations in *f* to be considered valid;
196 if a solution does not satisfy any equation, False is returned; if one or
197 more checks are inconclusive (and none are False) then None is returned.
199 Examples
200 ========
202 >>> from sympy import checksol, symbols
203 >>> x, y = symbols('x,y')
204 >>> checksol(x**4 - 1, x, 1)
205 True
206 >>> checksol(x**4 - 1, x, 0)
207 False
208 >>> checksol(x**2 + y**2 - 5**2, {x: 3, y: 4})
209 True
211 To check if an expression is zero using ``checksol()``, pass it
212 as *f* and send an empty dictionary for *symbol*:
214 >>> checksol(x**2 + x - x*(x + 1), {})
215 True
217 None is returned if ``checksol()`` could not conclude.
219 flags:
220 'numerical=True (default)'
221 do a fast numerical check if ``f`` has only one symbol.
222 'minimal=True (default is False)'
223 a very fast, minimal testing.
224 'warn=True (default is False)'
225 show a warning if checksol() could not conclude.
226 'simplify=True (default)'
227 simplify solution before substituting into function and
228 simplify the function before trying specific simplifications
229 'force=True (default is False)'
230 make positive all symbols without assumptions regarding sign.
232 """
233 from sympy.physics.units import Unit
235 minimal = flags.get('minimal', False)
237 if sol is not None:
238 sol = {symbol: sol}
239 elif isinstance(symbol, dict):
240 sol = symbol
241 else:
242 msg = 'Expecting (sym, val) or ({sym: val}, None) but got (%s, %s)'
243 raise ValueError(msg % (symbol, sol))
245 if iterable(f):
246 if not f:
247 raise ValueError('no functions to check')
248 rv = True
249 for fi in f:
250 check = checksol(fi, sol, **flags)
251 if check:
252 continue
253 if check is False:
254 return False
255 rv = None # don't return, wait to see if there's a False
256 return rv
258 f = _sympify(f)
260 if f.is_number:
261 return f.is_zero
263 if isinstance(f, Poly):
264 f = f.as_expr()
265 elif isinstance(f, (Eq, Ne)):
266 if f.rhs in (S.true, S.false):
267 f = f.reversed
268 B, E = f.args
269 if isinstance(B, BooleanAtom):
270 f = f.subs(sol)
271 if not f.is_Boolean:
272 return
273 else:
274 f = f.rewrite(Add, evaluate=False, deep=False)
276 if isinstance(f, BooleanAtom):
277 return bool(f)
278 elif not f.is_Relational and not f:
279 return True
281 illegal = set(_illegal)
282 if any(sympify(v).atoms() & illegal for k, v in sol.items()):
283 return False
285 attempt = -1
286 numerical = flags.get('numerical', True)
287 while 1:
288 attempt += 1
289 if attempt == 0:
290 val = f.subs(sol)
291 if isinstance(val, Mul):
292 val = val.as_independent(Unit)[0]
293 if val.atoms() & illegal:
294 return False
295 elif attempt == 1:
296 if not val.is_number:
297 if not val.is_constant(*list(sol.keys()), simplify=not minimal):
298 return False
299 # there are free symbols -- simple expansion might work
300 _, val = val.as_content_primitive()
301 val = _mexpand(val.as_numer_denom()[0], recursive=True)
302 elif attempt == 2:
303 if minimal:
304 return
305 if flags.get('simplify', True):
306 for k in sol:
307 sol[k] = simplify(sol[k])
308 # start over without the failed expanded form, possibly
309 # with a simplified solution
310 val = simplify(f.subs(sol))
311 if flags.get('force', True):
312 val, reps = posify(val)
313 # expansion may work now, so try again and check
314 exval = _mexpand(val, recursive=True)
315 if exval.is_number:
316 # we can decide now
317 val = exval
318 else:
319 # if there are no radicals and no functions then this can't be
320 # zero anymore -- can it?
321 pot = preorder_traversal(expand_mul(val))
322 seen = set()
323 saw_pow_func = False
324 for p in pot:
325 if p in seen:
326 continue
327 seen.add(p)
328 if p.is_Pow and not p.exp.is_Integer:
329 saw_pow_func = True
330 elif p.is_Function:
331 saw_pow_func = True
332 elif isinstance(p, UndefinedFunction):
333 saw_pow_func = True
334 if saw_pow_func:
335 break
336 if saw_pow_func is False:
337 return False
338 if flags.get('force', True):
339 # don't do a zero check with the positive assumptions in place
340 val = val.subs(reps)
341 nz = fuzzy_not(val.is_zero)
342 if nz is not None:
343 # issue 5673: nz may be True even when False
344 # so these are just hacks to keep a false positive
345 # from being returned
347 # HACK 1: LambertW (issue 5673)
348 if val.is_number and val.has(LambertW):
349 # don't eval this to verify solution since if we got here,
350 # numerical must be False
351 return None
353 # add other HACKs here if necessary, otherwise we assume
354 # the nz value is correct
355 return not nz
356 break
357 if val.is_Rational:
358 return val == 0
359 if numerical and val.is_number:
360 return (abs(val.n(18).n(12, chop=True)) < 1e-9) is S.true
362 if flags.get('warn', False):
363 warnings.warn("\n\tWarning: could not verify solution %s." % sol)
364 # returns None if it can't conclude
365 # TODO: improve solution testing
368def solve(f, *symbols, **flags):
369 r"""
370 Algebraically solves equations and systems of equations.
372 Explanation
373 ===========
375 Currently supported:
376 - polynomial
377 - transcendental
378 - piecewise combinations of the above
379 - systems of linear and polynomial equations
380 - systems containing relational expressions
381 - systems implied by undetermined coefficients
383 Examples
384 ========
386 The default output varies according to the input and might
387 be a list (possibly empty), a dictionary, a list of
388 dictionaries or tuples, or an expression involving relationals.
389 For specifics regarding different forms of output that may appear, see :ref:`solve_output`.
390 Let it suffice here to say that to obtain a uniform output from
391 `solve` use ``dict=True`` or ``set=True`` (see below).
393 >>> from sympy import solve, Poly, Eq, Matrix, Symbol
394 >>> from sympy.abc import x, y, z, a, b
396 The expressions that are passed can be Expr, Equality, or Poly
397 classes (or lists of the same); a Matrix is considered to be a
398 list of all the elements of the matrix:
400 >>> solve(x - 3, x)
401 [3]
402 >>> solve(Eq(x, 3), x)
403 [3]
404 >>> solve(Poly(x - 3), x)
405 [3]
406 >>> solve(Matrix([[x, x + y]]), x, y) == solve([x, x + y], x, y)
407 True
409 If no symbols are indicated to be of interest and the equation is
410 univariate, a list of values is returned; otherwise, the keys in
411 a dictionary will indicate which (of all the variables used in
412 the expression(s)) variables and solutions were found:
414 >>> solve(x**2 - 4)
415 [-2, 2]
416 >>> solve((x - a)*(y - b))
417 [{a: x}, {b: y}]
418 >>> solve([x - 3, y - 1])
419 {x: 3, y: 1}
420 >>> solve([x - 3, y**2 - 1])
421 [{x: 3, y: -1}, {x: 3, y: 1}]
423 If you pass symbols for which solutions are sought, the output will vary
424 depending on the number of symbols you passed, whether you are passing
425 a list of expressions or not, and whether a linear system was solved.
426 Uniform output is attained by using ``dict=True`` or ``set=True``.
428 >>> #### *** feel free to skip to the stars below *** ####
429 >>> from sympy import TableForm
430 >>> h = [None, ';|;'.join(['e', 's', 'solve(e, s)', 'solve(e, s, dict=True)',
431 ... 'solve(e, s, set=True)']).split(';')]
432 >>> t = []
433 >>> for e, s in [
434 ... (x - y, y),
435 ... (x - y, [x, y]),
436 ... (x**2 - y, [x, y]),
437 ... ([x - 3, y -1], [x, y]),
438 ... ]:
439 ... how = [{}, dict(dict=True), dict(set=True)]
440 ... res = [solve(e, s, **f) for f in how]
441 ... t.append([e, '|', s, '|'] + [res[0], '|', res[1], '|', res[2]])
442 ...
443 >>> # ******************************************************* #
444 >>> TableForm(t, headings=h, alignments="<")
445 e | s | solve(e, s) | solve(e, s, dict=True) | solve(e, s, set=True)
446 ---------------------------------------------------------------------------------------
447 x - y | y | [x] | [{y: x}] | ([y], {(x,)})
448 x - y | [x, y] | [(y, y)] | [{x: y}] | ([x, y], {(y, y)})
449 x**2 - y | [x, y] | [(x, x**2)] | [{y: x**2}] | ([x, y], {(x, x**2)})
450 [x - 3, y - 1] | [x, y] | {x: 3, y: 1} | [{x: 3, y: 1}] | ([x, y], {(3, 1)})
452 * If any equation does not depend on the symbol(s) given, it will be
453 eliminated from the equation set and an answer may be given
454 implicitly in terms of variables that were not of interest:
456 >>> solve([x - y, y - 3], x)
457 {x: y}
459 When you pass all but one of the free symbols, an attempt
460 is made to find a single solution based on the method of
461 undetermined coefficients. If it succeeds, a dictionary of values
462 is returned. If you want an algebraic solutions for one
463 or more of the symbols, pass the expression to be solved in a list:
465 >>> e = a*x + b - 2*x - 3
466 >>> solve(e, [a, b])
467 {a: 2, b: 3}
468 >>> solve([e], [a, b])
469 {a: -b/x + (2*x + 3)/x}
471 When there is no solution for any given symbol which will make all
472 expressions zero, the empty list is returned (or an empty set in
473 the tuple when ``set=True``):
475 >>> from sympy import sqrt
476 >>> solve(3, x)
477 []
478 >>> solve(x - 3, y)
479 []
480 >>> solve(sqrt(x) + 1, x, set=True)
481 ([x], set())
483 When an object other than a Symbol is given as a symbol, it is
484 isolated algebraically and an implicit solution may be obtained.
485 This is mostly provided as a convenience to save you from replacing
486 the object with a Symbol and solving for that Symbol. It will only
487 work if the specified object can be replaced with a Symbol using the
488 subs method:
490 >>> from sympy import exp, Function
491 >>> f = Function('f')
493 >>> solve(f(x) - x, f(x))
494 [x]
495 >>> solve(f(x).diff(x) - f(x) - x, f(x).diff(x))
496 [x + f(x)]
497 >>> solve(f(x).diff(x) - f(x) - x, f(x))
498 [-x + Derivative(f(x), x)]
499 >>> solve(x + exp(x)**2, exp(x), set=True)
500 ([exp(x)], {(-sqrt(-x),), (sqrt(-x),)})
502 >>> from sympy import Indexed, IndexedBase, Tuple
503 >>> A = IndexedBase('A')
504 >>> eqs = Tuple(A[1] + A[2] - 3, A[1] - A[2] + 1)
505 >>> solve(eqs, eqs.atoms(Indexed))
506 {A[1]: 1, A[2]: 2}
508 * To solve for a function within a derivative, use :func:`~.dsolve`.
510 To solve for a symbol implicitly, use implicit=True:
512 >>> solve(x + exp(x), x)
513 [-LambertW(1)]
514 >>> solve(x + exp(x), x, implicit=True)
515 [-exp(x)]
517 It is possible to solve for anything in an expression that can be
518 replaced with a symbol using :obj:`~sympy.core.basic.Basic.subs`:
520 >>> solve(x + 2 + sqrt(3), x + 2)
521 [-sqrt(3)]
522 >>> solve((x + 2 + sqrt(3), x + 4 + y), y, x + 2)
523 {y: -2 + sqrt(3), x + 2: -sqrt(3)}
525 * Nothing heroic is done in this implicit solving so you may end up
526 with a symbol still in the solution:
528 >>> eqs = (x*y + 3*y + sqrt(3), x + 4 + y)
529 >>> solve(eqs, y, x + 2)
530 {y: -sqrt(3)/(x + 3), x + 2: -2*x/(x + 3) - 6/(x + 3) + sqrt(3)/(x + 3)}
531 >>> solve(eqs, y*x, x)
532 {x: -y - 4, x*y: -3*y - sqrt(3)}
534 * If you attempt to solve for a number, remember that the number
535 you have obtained does not necessarily mean that the value is
536 equivalent to the expression obtained:
538 >>> solve(sqrt(2) - 1, 1)
539 [sqrt(2)]
540 >>> solve(x - y + 1, 1) # /!\ -1 is targeted, too
541 [x/(y - 1)]
542 >>> [_.subs(z, -1) for _ in solve((x - y + 1).subs(-1, z), 1)]
543 [-x + y]
545 **Additional Examples**
547 ``solve()`` with check=True (default) will run through the symbol tags to
548 eliminate unwanted solutions. If no assumptions are included, all possible
549 solutions will be returned:
551 >>> x = Symbol("x")
552 >>> solve(x**2 - 1)
553 [-1, 1]
555 By setting the ``positive`` flag, only one solution will be returned:
557 >>> pos = Symbol("pos", positive=True)
558 >>> solve(pos**2 - 1)
559 [1]
561 When the solutions are checked, those that make any denominator zero
562 are automatically excluded. If you do not want to exclude such solutions,
563 then use the check=False option:
565 >>> from sympy import sin, limit
566 >>> solve(sin(x)/x) # 0 is excluded
567 [pi]
569 If ``check=False``, then a solution to the numerator being zero is found
570 but the value of $x = 0$ is a spurious solution since $\sin(x)/x$ has the well
571 known limit (without discontinuity) of 1 at $x = 0$:
573 >>> solve(sin(x)/x, check=False)
574 [0, pi]
576 In the following case, however, the limit exists and is equal to the
577 value of $x = 0$ that is excluded when check=True:
579 >>> eq = x**2*(1/x - z**2/x)
580 >>> solve(eq, x)
581 []
582 >>> solve(eq, x, check=False)
583 [0]
584 >>> limit(eq, x, 0, '-')
585 0
586 >>> limit(eq, x, 0, '+')
587 0
589 **Solving Relationships**
591 When one or more expressions passed to ``solve`` is a relational,
592 a relational result is returned (and the ``dict`` and ``set`` flags
593 are ignored):
595 >>> solve(x < 3)
596 (-oo < x) & (x < 3)
597 >>> solve([x < 3, x**2 > 4], x)
598 ((-oo < x) & (x < -2)) | ((2 < x) & (x < 3))
599 >>> solve([x + y - 3, x > 3], x)
600 (3 < x) & (x < oo) & Eq(x, 3 - y)
602 Although checking of assumptions on symbols in relationals
603 is not done, setting assumptions will affect how certain
604 relationals might automatically simplify:
606 >>> solve(x**2 > 4)
607 ((-oo < x) & (x < -2)) | ((2 < x) & (x < oo))
609 >>> r = Symbol('r', real=True)
610 >>> solve(r**2 > 4)
611 (2 < r) | (r < -2)
613 There is currently no algorithm in SymPy that allows you to use
614 relationships to resolve more than one variable. So the following
615 does not determine that ``q < 0`` (and trying to solve for ``r``
616 and ``q`` will raise an error):
618 >>> from sympy import symbols
619 >>> r, q = symbols('r, q', real=True)
620 >>> solve([r + q - 3, r > 3], r)
621 (3 < r) & Eq(r, 3 - q)
623 You can directly call the routine that ``solve`` calls
624 when it encounters a relational: :func:`~.reduce_inequalities`.
625 It treats Expr like Equality.
627 >>> from sympy import reduce_inequalities
628 >>> reduce_inequalities([x**2 - 4])
629 Eq(x, -2) | Eq(x, 2)
631 If each relationship contains only one symbol of interest,
632 the expressions can be processed for multiple symbols:
634 >>> reduce_inequalities([0 <= x - 1, y < 3], [x, y])
635 (-oo < y) & (1 <= x) & (x < oo) & (y < 3)
637 But an error is raised if any relationship has more than one
638 symbol of interest:
640 >>> reduce_inequalities([0 <= x*y - 1, y < 3], [x, y])
641 Traceback (most recent call last):
642 ...
643 NotImplementedError:
644 inequality has more than one symbol of interest.
646 **Disabling High-Order Explicit Solutions**
648 When solving polynomial expressions, you might not want explicit solutions
649 (which can be quite long). If the expression is univariate, ``CRootOf``
650 instances will be returned instead:
652 >>> solve(x**3 - x + 1)
653 [-1/((-1/2 - sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)) -
654 (-1/2 - sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3,
655 -(-1/2 + sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3 -
656 1/((-1/2 + sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)),
657 -(3*sqrt(69)/2 + 27/2)**(1/3)/3 -
658 1/(3*sqrt(69)/2 + 27/2)**(1/3)]
659 >>> solve(x**3 - x + 1, cubics=False)
660 [CRootOf(x**3 - x + 1, 0),
661 CRootOf(x**3 - x + 1, 1),
662 CRootOf(x**3 - x + 1, 2)]
664 If the expression is multivariate, no solution might be returned:
666 >>> solve(x**3 - x + a, x, cubics=False)
667 []
669 Sometimes solutions will be obtained even when a flag is False because the
670 expression could be factored. In the following example, the equation can
671 be factored as the product of a linear and a quadratic factor so explicit
672 solutions (which did not require solving a cubic expression) are obtained:
674 >>> eq = x**3 + 3*x**2 + x - 1
675 >>> solve(eq, cubics=False)
676 [-1, -1 + sqrt(2), -sqrt(2) - 1]
678 **Solving Equations Involving Radicals**
680 Because of SymPy's use of the principle root, some solutions
681 to radical equations will be missed unless check=False:
683 >>> from sympy import root
684 >>> eq = root(x**3 - 3*x**2, 3) + 1 - x
685 >>> solve(eq)
686 []
687 >>> solve(eq, check=False)
688 [1/3]
690 In the above example, there is only a single solution to the
691 equation. Other expressions will yield spurious roots which
692 must be checked manually; roots which give a negative argument
693 to odd-powered radicals will also need special checking:
695 >>> from sympy import real_root, S
696 >>> eq = root(x, 3) - root(x, 5) + S(1)/7
697 >>> solve(eq) # this gives 2 solutions but misses a 3rd
698 [CRootOf(7*x**5 - 7*x**3 + 1, 1)**15,
699 CRootOf(7*x**5 - 7*x**3 + 1, 2)**15]
700 >>> sol = solve(eq, check=False)
701 >>> [abs(eq.subs(x,i).n(2)) for i in sol]
702 [0.48, 0.e-110, 0.e-110, 0.052, 0.052]
704 The first solution is negative so ``real_root`` must be used to see that it
705 satisfies the expression:
707 >>> abs(real_root(eq.subs(x, sol[0])).n(2))
708 0.e-110
710 If the roots of the equation are not real then more care will be
711 necessary to find the roots, especially for higher order equations.
712 Consider the following expression:
714 >>> expr = root(x, 3) - root(x, 5)
716 We will construct a known value for this expression at x = 3 by selecting
717 the 1-th root for each radical:
719 >>> expr1 = root(x, 3, 1) - root(x, 5, 1)
720 >>> v = expr1.subs(x, -3)
722 The ``solve`` function is unable to find any exact roots to this equation:
724 >>> eq = Eq(expr, v); eq1 = Eq(expr1, v)
725 >>> solve(eq, check=False), solve(eq1, check=False)
726 ([], [])
728 The function ``unrad``, however, can be used to get a form of the equation
729 for which numerical roots can be found:
731 >>> from sympy.solvers.solvers import unrad
732 >>> from sympy import nroots
733 >>> e, (p, cov) = unrad(eq)
734 >>> pvals = nroots(e)
735 >>> inversion = solve(cov, x)[0]
736 >>> xvals = [inversion.subs(p, i) for i in pvals]
738 Although ``eq`` or ``eq1`` could have been used to find ``xvals``, the
739 solution can only be verified with ``expr1``:
741 >>> z = expr - v
742 >>> [xi.n(chop=1e-9) for xi in xvals if abs(z.subs(x, xi).n()) < 1e-9]
743 []
744 >>> z1 = expr1 - v
745 >>> [xi.n(chop=1e-9) for xi in xvals if abs(z1.subs(x, xi).n()) < 1e-9]
746 [-3.0]
748 Parameters
749 ==========
751 f :
752 - a single Expr or Poly that must be zero
753 - an Equality
754 - a Relational expression
755 - a Boolean
756 - iterable of one or more of the above
758 symbols : (object(s) to solve for) specified as
759 - none given (other non-numeric objects will be used)
760 - single symbol
761 - denested list of symbols
762 (e.g., ``solve(f, x, y)``)
763 - ordered iterable of symbols
764 (e.g., ``solve(f, [x, y])``)
766 flags :
767 dict=True (default is False)
768 Return list (perhaps empty) of solution mappings.
769 set=True (default is False)
770 Return list of symbols and set of tuple(s) of solution(s).
771 exclude=[] (default)
772 Do not try to solve for any of the free symbols in exclude;
773 if expressions are given, the free symbols in them will
774 be extracted automatically.
775 check=True (default)
776 If False, do not do any testing of solutions. This can be
777 useful if you want to include solutions that make any
778 denominator zero.
779 numerical=True (default)
780 Do a fast numerical check if *f* has only one symbol.
781 minimal=True (default is False)
782 A very fast, minimal testing.
783 warn=True (default is False)
784 Show a warning if ``checksol()`` could not conclude.
785 simplify=True (default)
786 Simplify all but polynomials of order 3 or greater before
787 returning them and (if check is not False) use the
788 general simplify function on the solutions and the
789 expression obtained when they are substituted into the
790 function which should be zero.
791 force=True (default is False)
792 Make positive all symbols without assumptions regarding sign.
793 rational=True (default)
794 Recast Floats as Rational; if this option is not used, the
795 system containing Floats may fail to solve because of issues
796 with polys. If rational=None, Floats will be recast as
797 rationals but the answer will be recast as Floats. If the
798 flag is False then nothing will be done to the Floats.
799 manual=True (default is False)
800 Do not use the polys/matrix method to solve a system of
801 equations, solve them one at a time as you might "manually."
802 implicit=True (default is False)
803 Allows ``solve`` to return a solution for a pattern in terms of
804 other functions that contain that pattern; this is only
805 needed if the pattern is inside of some invertible function
806 like cos, exp, ect.
807 particular=True (default is False)
808 Instructs ``solve`` to try to find a particular solution to
809 a linear system with as many zeros as possible; this is very
810 expensive.
811 quick=True (default is False; ``particular`` must be True)
812 Selects a fast heuristic to find a solution with many zeros
813 whereas a value of False uses the very slow method guaranteed
814 to find the largest number of zeros possible.
815 cubics=True (default)
816 Return explicit solutions when cubic expressions are encountered.
817 When False, quartics and quintics are disabled, too.
818 quartics=True (default)
819 Return explicit solutions when quartic expressions are encountered.
820 When False, quintics are disabled, too.
821 quintics=True (default)
822 Return explicit solutions (if possible) when quintic expressions
823 are encountered.
825 See Also
826 ========
828 rsolve: For solving recurrence relationships
829 dsolve: For solving differential equations
831 """
832 from .inequalities import reduce_inequalities
834 # checking/recording flags
835 ###########################################################################
837 # set solver types explicitly; as soon as one is False
838 # all the rest will be False
839 hints = ('cubics', 'quartics', 'quintics')
840 default = True
841 for k in hints:
842 default = flags.setdefault(k, bool(flags.get(k, default)))
844 # allow solution to contain symbol if True:
845 implicit = flags.get('implicit', False)
847 # record desire to see warnings
848 warn = flags.get('warn', False)
850 # this flag will be needed for quick exits below, so record
851 # now -- but don't record `dict` yet since it might change
852 as_set = flags.get('set', False)
854 # keeping track of how f was passed
855 bare_f = not iterable(f)
857 # check flag usage for particular/quick which should only be used
858 # with systems of equations
859 if flags.get('quick', None) is not None:
860 if not flags.get('particular', None):
861 raise ValueError('when using `quick`, `particular` should be True')
862 if flags.get('particular', False) and bare_f:
863 raise ValueError(filldedent("""
864 The 'particular/quick' flag is usually used with systems of
865 equations. Either pass your equation in a list or
866 consider using a solver like `diophantine` if you are
867 looking for a solution in integers."""))
869 # sympify everything, creating list of expressions and list of symbols
870 ###########################################################################
872 def _sympified_list(w):
873 return list(map(sympify, w if iterable(w) else [w]))
874 f, symbols = (_sympified_list(w) for w in [f, symbols])
876 # preprocess symbol(s)
877 ###########################################################################
879 ordered_symbols = None # were the symbols in a well defined order?
880 if not symbols:
881 # get symbols from equations
882 symbols = set().union(*[fi.free_symbols for fi in f])
883 if len(symbols) < len(f):
884 for fi in f:
885 pot = preorder_traversal(fi)
886 for p in pot:
887 if isinstance(p, AppliedUndef):
888 if not as_set:
889 flags['dict'] = True # better show symbols
890 symbols.add(p)
891 pot.skip() # don't go any deeper
892 ordered_symbols = False
893 symbols = list(ordered(symbols)) # to make it canonical
894 else:
895 if len(symbols) == 1 and iterable(symbols[0]):
896 symbols = symbols[0]
897 ordered_symbols = symbols and is_sequence(symbols,
898 include=GeneratorType)
899 _symbols = list(uniq(symbols))
900 if len(_symbols) != len(symbols):
901 ordered_symbols = False
902 symbols = list(ordered(symbols))
903 else:
904 symbols = _symbols
906 # check for duplicates
907 if len(symbols) != len(set(symbols)):
908 raise ValueError('duplicate symbols given')
909 # remove those not of interest
910 exclude = flags.pop('exclude', set())
911 if exclude:
912 if isinstance(exclude, Expr):
913 exclude = [exclude]
914 exclude = set().union(*[e.free_symbols for e in sympify(exclude)])
915 symbols = [s for s in symbols if s not in exclude]
917 # preprocess equation(s)
918 ###########################################################################
920 # automatically ignore True values
921 if isinstance(f, list):
922 f = [s for s in f if s is not S.true]
924 # handle canonicalization of equation types
925 for i, fi in enumerate(f):
926 if isinstance(fi, (Eq, Ne)):
927 if 'ImmutableDenseMatrix' in [type(a).__name__ for a in fi.args]:
928 fi = fi.lhs - fi.rhs
929 else:
930 L, R = fi.args
931 if isinstance(R, BooleanAtom):
932 L, R = R, L
933 if isinstance(L, BooleanAtom):
934 if isinstance(fi, Ne):
935 L = ~L
936 if R.is_Relational:
937 fi = ~R if L is S.false else R
938 elif R.is_Symbol:
939 return L
940 elif R.is_Boolean and (~R).is_Symbol:
941 return ~L
942 else:
943 raise NotImplementedError(filldedent('''
944 Unanticipated argument of Eq when other arg
945 is True or False.
946 '''))
947 else:
948 fi = fi.rewrite(Add, evaluate=False, deep=False)
949 f[i] = fi
951 # *** dispatch and handle as a system of relationals
952 # **************************************************
953 if fi.is_Relational:
954 if len(symbols) != 1:
955 raise ValueError("can only solve for one symbol at a time")
956 if warn and symbols[0].assumptions0:
957 warnings.warn(filldedent("""
958 \tWarning: assumptions about variable '%s' are
959 not handled currently.""" % symbols[0]))
960 return reduce_inequalities(f, symbols=symbols)
962 # convert Poly to expression
963 if isinstance(fi, Poly):
964 f[i] = fi.as_expr()
966 # rewrite hyperbolics in terms of exp if they have symbols of
967 # interest
968 f[i] = f[i].replace(lambda w: isinstance(w, HyperbolicFunction) and \
969 w.has_free(*symbols), lambda w: w.rewrite(exp))
971 # if we have a Matrix, we need to iterate over its elements again
972 if f[i].is_Matrix:
973 bare_f = False
974 f.extend(list(f[i]))
975 f[i] = S.Zero
977 # if we can split it into real and imaginary parts then do so
978 freei = f[i].free_symbols
979 if freei and all(s.is_extended_real or s.is_imaginary for s in freei):
980 fr, fi = f[i].as_real_imag()
981 # accept as long as new re, im, arg or atan2 are not introduced
982 had = f[i].atoms(re, im, arg, atan2)
983 if fr and fi and fr != fi and not any(
984 i.atoms(re, im, arg, atan2) - had for i in (fr, fi)):
985 if bare_f:
986 bare_f = False
987 f[i: i + 1] = [fr, fi]
989 # real/imag handling -----------------------------
990 if any(isinstance(fi, (bool, BooleanAtom)) for fi in f):
991 if as_set:
992 return [], set()
993 return []
995 for i, fi in enumerate(f):
996 # Abs
997 while True:
998 was = fi
999 fi = fi.replace(Abs, lambda arg:
1000 separatevars(Abs(arg)).rewrite(Piecewise) if arg.has(*symbols)
1001 else Abs(arg))
1002 if was == fi:
1003 break
1005 for e in fi.find(Abs):
1006 if e.has(*symbols):
1007 raise NotImplementedError('solving %s when the argument '
1008 'is not real or imaginary.' % e)
1010 # arg
1011 fi = fi.replace(arg, lambda a: arg(a).rewrite(atan2).rewrite(atan))
1013 # save changes
1014 f[i] = fi
1016 # see if re(s) or im(s) appear
1017 freim = [fi for fi in f if fi.has(re, im)]
1018 if freim:
1019 irf = []
1020 for s in symbols:
1021 if s.is_real or s.is_imaginary:
1022 continue # neither re(x) nor im(x) will appear
1023 # if re(s) or im(s) appear, the auxiliary equation must be present
1024 if any(fi.has(re(s), im(s)) for fi in freim):
1025 irf.append((s, re(s) + S.ImaginaryUnit*im(s)))
1026 if irf:
1027 for s, rhs in irf:
1028 f = [fi.xreplace({s: rhs}) for fi in f] + [s - rhs]
1029 symbols.extend([re(s), im(s)])
1030 if bare_f:
1031 bare_f = False
1032 flags['dict'] = True
1033 # end of real/imag handling -----------------------------
1035 # we can solve for non-symbol entities by replacing them with Dummy symbols
1036 f, symbols, swap_sym = recast_to_symbols(f, symbols)
1037 # this set of symbols (perhaps recast) is needed below
1038 symset = set(symbols)
1040 # get rid of equations that have no symbols of interest; we don't
1041 # try to solve them because the user didn't ask and they might be
1042 # hard to solve; this means that solutions may be given in terms
1043 # of the eliminated equations e.g. solve((x-y, y-3), x) -> {x: y}
1044 newf = []
1045 for fi in f:
1046 # let the solver handle equations that..
1047 # - have no symbols but are expressions
1048 # - have symbols of interest
1049 # - have no symbols of interest but are constant
1050 # but when an expression is not constant and has no symbols of
1051 # interest, it can't change what we obtain for a solution from
1052 # the remaining equations so we don't include it; and if it's
1053 # zero it can be removed and if it's not zero, there is no
1054 # solution for the equation set as a whole
1055 #
1056 # The reason for doing this filtering is to allow an answer
1057 # to be obtained to queries like solve((x - y, y), x); without
1058 # this mod the return value is []
1059 ok = False
1060 if fi.free_symbols & symset:
1061 ok = True
1062 else:
1063 if fi.is_number:
1064 if fi.is_Number:
1065 if fi.is_zero:
1066 continue
1067 return []
1068 ok = True
1069 else:
1070 if fi.is_constant():
1071 ok = True
1072 if ok:
1073 newf.append(fi)
1074 if not newf:
1075 if as_set:
1076 return symbols, set()
1077 return []
1078 f = newf
1079 del newf
1081 # mask off any Object that we aren't going to invert: Derivative,
1082 # Integral, etc... so that solving for anything that they contain will
1083 # give an implicit solution
1084 seen = set()
1085 non_inverts = set()
1086 for fi in f:
1087 pot = preorder_traversal(fi)
1088 for p in pot:
1089 if not isinstance(p, Expr) or isinstance(p, Piecewise):
1090 pass
1091 elif (isinstance(p, bool) or
1092 not p.args or
1093 p in symset or
1094 p.is_Add or p.is_Mul or
1095 p.is_Pow and not implicit or
1096 p.is_Function and not implicit) and p.func not in (re, im):
1097 continue
1098 elif p not in seen:
1099 seen.add(p)
1100 if p.free_symbols & symset:
1101 non_inverts.add(p)
1102 else:
1103 continue
1104 pot.skip()
1105 del seen
1106 non_inverts = dict(list(zip(non_inverts, [Dummy() for _ in non_inverts])))
1107 f = [fi.subs(non_inverts) for fi in f]
1109 # Both xreplace and subs are needed below: xreplace to force substitution
1110 # inside Derivative, subs to handle non-straightforward substitutions
1111 non_inverts = [(v, k.xreplace(swap_sym).subs(swap_sym)) for k, v in non_inverts.items()]
1113 # rationalize Floats
1114 floats = False
1115 if flags.get('rational', True) is not False:
1116 for i, fi in enumerate(f):
1117 if fi.has(Float):
1118 floats = True
1119 f[i] = nsimplify(fi, rational=True)
1121 # capture any denominators before rewriting since
1122 # they may disappear after the rewrite, e.g. issue 14779
1123 flags['_denominators'] = _simple_dens(f[0], symbols)
1125 # Any embedded piecewise functions need to be brought out to the
1126 # top level so that the appropriate strategy gets selected.
1127 # However, this is necessary only if one of the piecewise
1128 # functions depends on one of the symbols we are solving for.
1129 def _has_piecewise(e):
1130 if e.is_Piecewise:
1131 return e.has(*symbols)
1132 return any(_has_piecewise(a) for a in e.args)
1133 for i, fi in enumerate(f):
1134 if _has_piecewise(fi):
1135 f[i] = piecewise_fold(fi)
1137 #
1138 # try to get a solution
1139 ###########################################################################
1140 if bare_f:
1141 solution = None
1142 if len(symbols) != 1:
1143 solution = _solve_undetermined(f[0], symbols, flags)
1144 if not solution:
1145 solution = _solve(f[0], *symbols, **flags)
1146 else:
1147 linear, solution = _solve_system(f, symbols, **flags)
1148 assert type(solution) is list
1149 assert not solution or type(solution[0]) is dict, solution
1150 #
1151 # postprocessing
1152 ###########################################################################
1153 # capture as_dict flag now (as_set already captured)
1154 as_dict = flags.get('dict', False)
1156 # define how solution will get unpacked
1157 tuple_format = lambda s: [tuple([i.get(x, x) for x in symbols]) for i in s]
1158 if as_dict or as_set:
1159 unpack = None
1160 elif bare_f:
1161 if len(symbols) == 1:
1162 unpack = lambda s: [i[symbols[0]] for i in s]
1163 elif len(solution) == 1 and len(solution[0]) == len(symbols):
1164 # undetermined linear coeffs solution
1165 unpack = lambda s: s[0]
1166 elif ordered_symbols:
1167 unpack = tuple_format
1168 else:
1169 unpack = lambda s: s
1170 else:
1171 if solution:
1172 if linear and len(solution) == 1:
1173 # if you want the tuple solution for the linear
1174 # case, use `set=True`
1175 unpack = lambda s: s[0]
1176 elif ordered_symbols:
1177 unpack = tuple_format
1178 else:
1179 unpack = lambda s: s
1180 else:
1181 unpack = None
1183 # Restore masked-off objects
1184 if non_inverts and type(solution) is list:
1185 solution = [{k: v.subs(non_inverts) for k, v in s.items()}
1186 for s in solution]
1188 # Restore original "symbols" if a dictionary is returned.
1189 # This is not necessary for
1190 # - the single univariate equation case
1191 # since the symbol will have been removed from the solution;
1192 # - the nonlinear poly_system since that only supports zero-dimensional
1193 # systems and those results come back as a list
1194 #
1195 # ** unless there were Derivatives with the symbols, but those were handled
1196 # above.
1197 if swap_sym:
1198 symbols = [swap_sym.get(k, k) for k in symbols]
1199 for i, sol in enumerate(solution):
1200 solution[i] = {swap_sym.get(k, k): v.subs(swap_sym)
1201 for k, v in sol.items()}
1203 # Get assumptions about symbols, to filter solutions.
1204 # Note that if assumptions about a solution can't be verified, it is still
1205 # returned.
1206 check = flags.get('check', True)
1208 # restore floats
1209 if floats and solution and flags.get('rational', None) is None:
1210 solution = nfloat(solution, exponent=False)
1211 # nfloat might reveal more duplicates
1212 solution = _remove_duplicate_solutions(solution)
1214 if check and solution: # assumption checking
1215 warn = flags.get('warn', False)
1216 got_None = [] # solutions for which one or more symbols gave None
1217 no_False = [] # solutions for which no symbols gave False
1218 for sol in solution:
1219 a_None = False
1220 for symb, val in sol.items():
1221 test = check_assumptions(val, **symb.assumptions0)
1222 if test:
1223 continue
1224 if test is False:
1225 break
1226 a_None = True
1227 else:
1228 no_False.append(sol)
1229 if a_None:
1230 got_None.append(sol)
1232 solution = no_False
1233 if warn and got_None:
1234 warnings.warn(filldedent("""
1235 \tWarning: assumptions concerning following solution(s)
1236 cannot be checked:""" + '\n\t' +
1237 ', '.join(str(s) for s in got_None)))
1239 #
1240 # done
1241 ###########################################################################
1243 if not solution:
1244 if as_set:
1245 return symbols, set()
1246 return []
1248 # make orderings canonical for list of dictionaries
1249 if not as_set: # for set, no point in ordering
1250 solution = [{k: s[k] for k in ordered(s)} for s in solution]
1251 solution.sort(key=default_sort_key)
1253 if not (as_set or as_dict):
1254 return unpack(solution)
1256 if as_dict:
1257 return solution
1259 # set output: (symbols, {t1, t2, ...}) from list of dictionaries;
1260 # include all symbols for those that like a verbose solution
1261 # and to resolve any differences in dictionary keys.
1262 #
1263 # The set results can easily be used to make a verbose dict as
1264 # k, v = solve(eqs, syms, set=True)
1265 # sol = [dict(zip(k,i)) for i in v]
1266 #
1267 if ordered_symbols:
1268 k = symbols # keep preferred order
1269 else:
1270 # just unify the symbols for which solutions were found
1271 k = list(ordered(set(flatten(tuple(i.keys()) for i in solution))))
1272 return k, {tuple([s.get(ki, ki) for ki in k]) for s in solution}
1275def _solve_undetermined(g, symbols, flags):
1276 """solve helper to return a list with one dict (solution) else None
1278 A direct call to solve_undetermined_coeffs is more flexible and
1279 can return both multiple solutions and handle more than one independent
1280 variable. Here, we have to be more cautious to keep from solving
1281 something that does not look like an undetermined coeffs system --
1282 to minimize the surprise factor since singularities that cancel are not
1283 prohibited in solve_undetermined_coeffs.
1284 """
1285 if g.free_symbols - set(symbols):
1286 sol = solve_undetermined_coeffs(g, symbols, **dict(flags, dict=True, set=None))
1287 if len(sol) == 1:
1288 return sol
1291def _solve(f, *symbols, **flags):
1292 """Return a checked solution for *f* in terms of one or more of the
1293 symbols in the form of a list of dictionaries.
1295 If no method is implemented to solve the equation, a NotImplementedError
1296 will be raised. In the case that conversion of an expression to a Poly
1297 gives None a ValueError will be raised.
1298 """
1300 not_impl_msg = "No algorithms are implemented to solve equation %s"
1302 if len(symbols) != 1:
1303 # look for solutions for desired symbols that are independent
1304 # of symbols already solved for, e.g. if we solve for x = y
1305 # then no symbol having x in its solution will be returned.
1307 # First solve for linear symbols (since that is easier and limits
1308 # solution size) and then proceed with symbols appearing
1309 # in a non-linear fashion. Ideally, if one is solving a single
1310 # expression for several symbols, they would have to be
1311 # appear in factors of an expression, but we do not here
1312 # attempt factorization. XXX perhaps handling a Mul
1313 # should come first in this routine whether there is
1314 # one or several symbols.
1315 nonlin_s = []
1316 got_s = set()
1317 rhs_s = set()
1318 result = []
1319 for s in symbols:
1320 xi, v = solve_linear(f, symbols=[s])
1321 if xi == s:
1322 # no need to check but we should simplify if desired
1323 if flags.get('simplify', True):
1324 v = simplify(v)
1325 vfree = v.free_symbols
1326 if vfree & got_s:
1327 # was linear, but has redundant relationship
1328 # e.g. x - y = 0 has y == x is redundant for x == y
1329 # so ignore
1330 continue
1331 rhs_s |= vfree
1332 got_s.add(xi)
1333 result.append({xi: v})
1334 elif xi: # there might be a non-linear solution if xi is not 0
1335 nonlin_s.append(s)
1336 if not nonlin_s:
1337 return result
1338 for s in nonlin_s:
1339 try:
1340 soln = _solve(f, s, **flags)
1341 for sol in soln:
1342 if sol[s].free_symbols & got_s:
1343 # depends on previously solved symbols: ignore
1344 continue
1345 got_s.add(s)
1346 result.append(sol)
1347 except NotImplementedError:
1348 continue
1349 if got_s:
1350 return result
1351 else:
1352 raise NotImplementedError(not_impl_msg % f)
1354 # solve f for a single variable
1356 symbol = symbols[0]
1358 # expand binomials only if it has the unknown symbol
1359 f = f.replace(lambda e: isinstance(e, binomial) and e.has(symbol),
1360 lambda e: expand_func(e))
1362 # checking will be done unless it is turned off before making a
1363 # recursive call; the variables `checkdens` and `check` are
1364 # captured here (for reference below) in case flag value changes
1365 flags['check'] = checkdens = check = flags.pop('check', True)
1367 # build up solutions if f is a Mul
1368 if f.is_Mul:
1369 result = set()
1370 for m in f.args:
1371 if m in {S.NegativeInfinity, S.ComplexInfinity, S.Infinity}:
1372 result = set()
1373 break
1374 soln = _vsolve(m, symbol, **flags)
1375 result.update(set(soln))
1376 result = [{symbol: v} for v in result]
1377 if check:
1378 # all solutions have been checked but now we must
1379 # check that the solutions do not set denominators
1380 # in any factor to zero
1381 dens = flags.get('_denominators', _simple_dens(f, symbols))
1382 result = [s for s in result if
1383 not any(checksol(den, s, **flags) for den in
1384 dens)]
1385 # set flags for quick exit at end; solutions for each
1386 # factor were already checked and simplified
1387 check = False
1388 flags['simplify'] = False
1390 elif f.is_Piecewise:
1391 result = set()
1392 for i, (expr, cond) in enumerate(f.args):
1393 if expr.is_zero:
1394 raise NotImplementedError(
1395 'solve cannot represent interval solutions')
1396 candidates = _vsolve(expr, symbol, **flags)
1397 # the explicit condition for this expr is the current cond
1398 # and none of the previous conditions
1399 args = [~c for _, c in f.args[:i]] + [cond]
1400 cond = And(*args)
1401 for candidate in candidates:
1402 if candidate in result:
1403 # an unconditional value was already there
1404 continue
1405 try:
1406 v = cond.subs(symbol, candidate)
1407 _eval_simplify = getattr(v, '_eval_simplify', None)
1408 if _eval_simplify is not None:
1409 # unconditionally take the simplification of v
1410 v = _eval_simplify(ratio=2, measure=lambda x: 1)
1411 except TypeError:
1412 # incompatible type with condition(s)
1413 continue
1414 if v == False:
1415 continue
1416 if v == True:
1417 result.add(candidate)
1418 else:
1419 result.add(Piecewise(
1420 (candidate, v),
1421 (S.NaN, True)))
1422 # solutions already checked and simplified
1423 # ****************************************
1424 return [{symbol: r} for r in result]
1425 else:
1426 # first see if it really depends on symbol and whether there
1427 # is only a linear solution
1428 f_num, sol = solve_linear(f, symbols=symbols)
1429 if f_num.is_zero or sol is S.NaN:
1430 return []
1431 elif f_num.is_Symbol:
1432 # no need to check but simplify if desired
1433 if flags.get('simplify', True):
1434 sol = simplify(sol)
1435 return [{f_num: sol}]
1437 poly = None
1438 # check for a single Add generator
1439 if not f_num.is_Add:
1440 add_args = [i for i in f_num.atoms(Add)
1441 if symbol in i.free_symbols]
1442 if len(add_args) == 1:
1443 gen = add_args[0]
1444 spart = gen.as_independent(symbol)[1].as_base_exp()[0]
1445 if spart == symbol:
1446 try:
1447 poly = Poly(f_num, spart)
1448 except PolynomialError:
1449 pass
1451 result = False # no solution was obtained
1452 msg = '' # there is no failure message
1454 # Poly is generally robust enough to convert anything to
1455 # a polynomial and tell us the different generators that it
1456 # contains, so we will inspect the generators identified by
1457 # polys to figure out what to do.
1459 # try to identify a single generator that will allow us to solve this
1460 # as a polynomial, followed (perhaps) by a change of variables if the
1461 # generator is not a symbol
1463 try:
1464 if poly is None:
1465 poly = Poly(f_num)
1466 if poly is None:
1467 raise ValueError('could not convert %s to Poly' % f_num)
1468 except GeneratorsNeeded:
1469 simplified_f = simplify(f_num)
1470 if simplified_f != f_num:
1471 return _solve(simplified_f, symbol, **flags)
1472 raise ValueError('expression appears to be a constant')
1474 gens = [g for g in poly.gens if g.has(symbol)]
1476 def _as_base_q(x):
1477 """Return (b**e, q) for x = b**(p*e/q) where p/q is the leading
1478 Rational of the exponent of x, e.g. exp(-2*x/3) -> (exp(x), 3)
1479 """
1480 b, e = x.as_base_exp()
1481 if e.is_Rational:
1482 return b, e.q
1483 if not e.is_Mul:
1484 return x, 1
1485 c, ee = e.as_coeff_Mul()
1486 if c.is_Rational and c is not S.One: # c could be a Float
1487 return b**ee, c.q
1488 return x, 1
1490 if len(gens) > 1:
1491 # If there is more than one generator, it could be that the
1492 # generators have the same base but different powers, e.g.
1493 # >>> Poly(exp(x) + 1/exp(x))
1494 # Poly(exp(-x) + exp(x), exp(-x), exp(x), domain='ZZ')
1495 #
1496 # If unrad was not disabled then there should be no rational
1497 # exponents appearing as in
1498 # >>> Poly(sqrt(x) + sqrt(sqrt(x)))
1499 # Poly(sqrt(x) + x**(1/4), sqrt(x), x**(1/4), domain='ZZ')
1501 bases, qs = list(zip(*[_as_base_q(g) for g in gens]))
1502 bases = set(bases)
1504 if len(bases) > 1 or not all(q == 1 for q in qs):
1505 funcs = {b for b in bases if b.is_Function}
1507 trig = {_ for _ in funcs if
1508 isinstance(_, TrigonometricFunction)}
1509 other = funcs - trig
1510 if not other and len(funcs.intersection(trig)) > 1:
1511 newf = None
1512 if f_num.is_Add and len(f_num.args) == 2:
1513 # check for sin(x)**p = cos(x)**p
1514 _args = f_num.args
1515 t = a, b = [i.atoms(Function).intersection(
1516 trig) for i in _args]
1517 if all(len(i) == 1 for i in t):
1518 a, b = [i.pop() for i in t]
1519 if isinstance(a, cos):
1520 a, b = b, a
1521 _args = _args[::-1]
1522 if isinstance(a, sin) and isinstance(b, cos
1523 ) and a.args[0] == b.args[0]:
1524 # sin(x) + cos(x) = 0 -> tan(x) + 1 = 0
1525 newf, _d = (TR2i(_args[0]/_args[1]) + 1
1526 ).as_numer_denom()
1527 if not _d.is_Number:
1528 newf = None
1529 if newf is None:
1530 newf = TR1(f_num).rewrite(tan)
1531 if newf != f_num:
1532 # don't check the rewritten form --check
1533 # solutions in the un-rewritten form below
1534 flags['check'] = False
1535 result = _solve(newf, symbol, **flags)
1536 flags['check'] = check
1538 # just a simple case - see if replacement of single function
1539 # clears all symbol-dependent functions, e.g.
1540 # log(x) - log(log(x) - 1) - 3 can be solved even though it has
1541 # two generators.
1543 if result is False and funcs:
1544 funcs = list(ordered(funcs)) # put shallowest function first
1545 f1 = funcs[0]
1546 t = Dummy('t')
1547 # perform the substitution
1548 ftry = f_num.subs(f1, t)
1550 # if no Functions left, we can proceed with usual solve
1551 if not ftry.has(symbol):
1552 cv_sols = _solve(ftry, t, **flags)
1553 cv_inv = list(ordered(_vsolve(t - f1, symbol, **flags)))[0]
1554 result = [{symbol: cv_inv.subs(sol)} for sol in cv_sols]
1556 if result is False:
1557 msg = 'multiple generators %s' % gens
1559 else:
1560 # e.g. case where gens are exp(x), exp(-x)
1561 u = bases.pop()
1562 t = Dummy('t')
1563 inv = _vsolve(u - t, symbol, **flags)
1564 if isinstance(u, (Pow, exp)):
1565 # this will be resolved by factor in _tsolve but we might
1566 # as well try a simple expansion here to get things in
1567 # order so something like the following will work now without
1568 # having to factor:
1569 #
1570 # >>> eq = (exp(I*(-x-2))+exp(I*(x+2)))
1571 # >>> eq.subs(exp(x),y) # fails
1572 # exp(I*(-x - 2)) + exp(I*(x + 2))
1573 # >>> eq.expand().subs(exp(x),y) # works
1574 # y**I*exp(2*I) + y**(-I)*exp(-2*I)
1575 def _expand(p):
1576 b, e = p.as_base_exp()
1577 e = expand_mul(e)
1578 return expand_power_exp(b**e)
1579 ftry = f_num.replace(
1580 lambda w: w.is_Pow or isinstance(w, exp),
1581 _expand).subs(u, t)
1582 if not ftry.has(symbol):
1583 soln = _solve(ftry, t, **flags)
1584 result = [{symbol: i.subs(s)} for i in inv for s in soln]
1586 elif len(gens) == 1:
1588 # There is only one generator that we are interested in, but
1589 # there may have been more than one generator identified by
1590 # polys (e.g. for symbols other than the one we are interested
1591 # in) so recast the poly in terms of our generator of interest.
1592 # Also use composite=True with f_num since Poly won't update
1593 # poly as documented in issue 8810.
1595 poly = Poly(f_num, gens[0], composite=True)
1597 # if we aren't on the tsolve-pass, use roots
1598 if not flags.pop('tsolve', False):
1599 soln = None
1600 deg = poly.degree()
1601 flags['tsolve'] = True
1602 hints = ('cubics', 'quartics', 'quintics')
1603 solvers = {h: flags.get(h) for h in hints}
1604 soln = roots(poly, **solvers)
1605 if sum(soln.values()) < deg:
1606 # e.g. roots(32*x**5 + 400*x**4 + 2032*x**3 +
1607 # 5000*x**2 + 6250*x + 3189) -> {}
1608 # so all_roots is used and RootOf instances are
1609 # returned *unless* the system is multivariate
1610 # or high-order EX domain.
1611 try:
1612 soln = poly.all_roots()
1613 except NotImplementedError:
1614 if not flags.get('incomplete', True):
1615 raise NotImplementedError(
1616 filldedent('''
1617 Neither high-order multivariate polynomials
1618 nor sorting of EX-domain polynomials is supported.
1619 If you want to see any results, pass keyword incomplete=True to
1620 solve; to see numerical values of roots
1621 for univariate expressions, use nroots.
1622 '''))
1623 else:
1624 pass
1625 else:
1626 soln = list(soln.keys())
1628 if soln is not None:
1629 u = poly.gen
1630 if u != symbol:
1631 try:
1632 t = Dummy('t')
1633 inv = _vsolve(u - t, symbol, **flags)
1634 soln = {i.subs(t, s) for i in inv for s in soln}
1635 except NotImplementedError:
1636 # perhaps _tsolve can handle f_num
1637 soln = None
1638 else:
1639 check = False # only dens need to be checked
1640 if soln is not None:
1641 if len(soln) > 2:
1642 # if the flag wasn't set then unset it since high-order
1643 # results are quite long. Perhaps one could base this
1644 # decision on a certain critical length of the
1645 # roots. In addition, wester test M2 has an expression
1646 # whose roots can be shown to be real with the
1647 # unsimplified form of the solution whereas only one of
1648 # the simplified forms appears to be real.
1649 flags['simplify'] = flags.get('simplify', False)
1650 if soln is not None:
1651 result = [{symbol: v} for v in soln]
1653 # fallback if above fails
1654 # -----------------------
1655 if result is False:
1656 # try unrad
1657 if flags.pop('_unrad', True):
1658 try:
1659 u = unrad(f_num, symbol)
1660 except (ValueError, NotImplementedError):
1661 u = False
1662 if u:
1663 eq, cov = u
1664 if cov:
1665 isym, ieq = cov
1666 inv = _vsolve(ieq, symbol, **flags)[0]
1667 rv = {inv.subs(xi) for xi in _solve(eq, isym, **flags)}
1668 else:
1669 try:
1670 rv = set(_vsolve(eq, symbol, **flags))
1671 except NotImplementedError:
1672 rv = None
1673 if rv is not None:
1674 result = [{symbol: v} for v in rv]
1675 # if the flag wasn't set then unset it since unrad results
1676 # can be quite long or of very high order
1677 flags['simplify'] = flags.get('simplify', False)
1678 else:
1679 pass # for coverage
1681 # try _tsolve
1682 if result is False:
1683 flags.pop('tsolve', None) # allow tsolve to be used on next pass
1684 try:
1685 soln = _tsolve(f_num, symbol, **flags)
1686 if soln is not None:
1687 result = [{symbol: v} for v in soln]
1688 except PolynomialError:
1689 pass
1690 # ----------- end of fallback ----------------------------
1692 if result is False:
1693 raise NotImplementedError('\n'.join([msg, not_impl_msg % f]))
1695 result = _remove_duplicate_solutions(result)
1697 if flags.get('simplify', True):
1698 result = [{k: d[k].simplify() for k in d} for d in result]
1699 # Simplification might reveal more duplicates
1700 result = _remove_duplicate_solutions(result)
1701 # we just simplified the solution so we now set the flag to
1702 # False so the simplification doesn't happen again in checksol()
1703 flags['simplify'] = False
1705 if checkdens:
1706 # reject any result that makes any denom. affirmatively 0;
1707 # if in doubt, keep it
1708 dens = _simple_dens(f, symbols)
1709 result = [r for r in result if
1710 not any(checksol(d, r, **flags)
1711 for d in dens)]
1712 if check:
1713 # keep only results if the check is not False
1714 result = [r for r in result if
1715 checksol(f_num, r, **flags) is not False]
1716 return result
1719def _remove_duplicate_solutions(solutions: list[dict[Expr, Expr]]
1720 ) -> list[dict[Expr, Expr]]:
1721 """Remove duplicates from a list of dicts"""
1722 solutions_set = set()
1723 solutions_new = []
1725 for sol in solutions:
1726 solset = frozenset(sol.items())
1727 if solset not in solutions_set:
1728 solutions_new.append(sol)
1729 solutions_set.add(solset)
1731 return solutions_new
1734def _solve_system(exprs, symbols, **flags):
1735 """return ``(linear, solution)`` where ``linear`` is True
1736 if the system was linear, else False; ``solution``
1737 is a list of dictionaries giving solutions for the symbols
1738 """
1739 if not exprs:
1740 return False, []
1742 if flags.pop('_split', True):
1743 # Split the system into connected components
1744 V = exprs
1745 symsset = set(symbols)
1746 exprsyms = {e: e.free_symbols & symsset for e in exprs}
1747 E = []
1748 sym_indices = {sym: i for i, sym in enumerate(symbols)}
1749 for n, e1 in enumerate(exprs):
1750 for e2 in exprs[:n]:
1751 # Equations are connected if they share a symbol
1752 if exprsyms[e1] & exprsyms[e2]:
1753 E.append((e1, e2))
1754 G = V, E
1755 subexprs = connected_components(G)
1756 if len(subexprs) > 1:
1757 subsols = []
1758 linear = True
1759 for subexpr in subexprs:
1760 subsyms = set()
1761 for e in subexpr:
1762 subsyms |= exprsyms[e]
1763 subsyms = sorted(subsyms, key = lambda x: sym_indices[x])
1764 flags['_split'] = False # skip split step
1765 _linear, subsol = _solve_system(subexpr, subsyms, **flags)
1766 if linear:
1767 linear = linear and _linear
1768 if not isinstance(subsol, list):
1769 subsol = [subsol]
1770 subsols.append(subsol)
1771 # Full solution is cartesion product of subsystems
1772 sols = []
1773 for soldicts in product(*subsols):
1774 sols.append(dict(item for sd in soldicts
1775 for item in sd.items()))
1776 return linear, sols
1778 polys = []
1779 dens = set()
1780 failed = []
1781 result = []
1782 solved_syms = []
1783 linear = True
1784 manual = flags.get('manual', False)
1785 checkdens = check = flags.get('check', True)
1787 for j, g in enumerate(exprs):
1788 dens.update(_simple_dens(g, symbols))
1789 i, d = _invert(g, *symbols)
1790 if d in symbols:
1791 if linear:
1792 linear = solve_linear(g, 0, [d])[0] == d
1793 g = d - i
1794 g = g.as_numer_denom()[0]
1795 if manual:
1796 failed.append(g)
1797 continue
1799 poly = g.as_poly(*symbols, extension=True)
1801 if poly is not None:
1802 polys.append(poly)
1803 else:
1804 failed.append(g)
1806 if polys:
1807 if all(p.is_linear for p in polys):
1808 n, m = len(polys), len(symbols)
1809 matrix = zeros(n, m + 1)
1811 for i, poly in enumerate(polys):
1812 for monom, coeff in poly.terms():
1813 try:
1814 j = monom.index(1)
1815 matrix[i, j] = coeff
1816 except ValueError:
1817 matrix[i, m] = -coeff
1819 # returns a dictionary ({symbols: values}) or None
1820 if flags.pop('particular', False):
1821 result = minsolve_linear_system(matrix, *symbols, **flags)
1822 else:
1823 result = solve_linear_system(matrix, *symbols, **flags)
1824 result = [result] if result else []
1825 if failed:
1826 if result:
1827 solved_syms = list(result[0].keys()) # there is only one result dict
1828 else:
1829 solved_syms = []
1830 # linear doesn't change
1831 else:
1832 linear = False
1833 if len(symbols) > len(polys):
1835 free = set().union(*[p.free_symbols for p in polys])
1836 free = list(ordered(free.intersection(symbols)))
1837 got_s = set()
1838 result = []
1839 for syms in subsets(free, len(polys)):
1840 try:
1841 # returns [], None or list of tuples
1842 res = solve_poly_system(polys, *syms)
1843 if res:
1844 for r in set(res):
1845 skip = False
1846 for r1 in r:
1847 if got_s and any(ss in r1.free_symbols
1848 for ss in got_s):
1849 # sol depends on previously
1850 # solved symbols: discard it
1851 skip = True
1852 if not skip:
1853 got_s.update(syms)
1854 result.append(dict(list(zip(syms, r))))
1855 except NotImplementedError:
1856 pass
1857 if got_s:
1858 solved_syms = list(got_s)
1859 else:
1860 raise NotImplementedError('no valid subset found')
1861 else:
1862 try:
1863 result = solve_poly_system(polys, *symbols)
1864 if result:
1865 solved_syms = symbols
1866 result = [dict(list(zip(solved_syms, r))) for r in set(result)]
1867 except NotImplementedError:
1868 failed.extend([g.as_expr() for g in polys])
1869 solved_syms = []
1871 # convert None or [] to [{}]
1872 result = result or [{}]
1874 if failed:
1875 linear = False
1876 # For each failed equation, see if we can solve for one of the
1877 # remaining symbols from that equation. If so, we update the
1878 # solution set and continue with the next failed equation,
1879 # repeating until we are done or we get an equation that can't
1880 # be solved.
1881 def _ok_syms(e, sort=False):
1882 rv = e.free_symbols & legal
1884 # Solve first for symbols that have lower degree in the equation.
1885 # Ideally we want to solve firstly for symbols that appear linearly
1886 # with rational coefficients e.g. if e = x*y + z then we should
1887 # solve for z first.
1888 def key(sym):
1889 ep = e.as_poly(sym)
1890 if ep is None:
1891 complexity = (S.Infinity, S.Infinity, S.Infinity)
1892 else:
1893 coeff_syms = ep.LC().free_symbols
1894 complexity = (ep.degree(), len(coeff_syms & rv), len(coeff_syms))
1895 return complexity + (default_sort_key(sym),)
1897 if sort:
1898 rv = sorted(rv, key=key)
1899 return rv
1901 legal = set(symbols) # what we are interested in
1902 # sort so equation with the fewest potential symbols is first
1903 u = Dummy() # used in solution checking
1904 for eq in ordered(failed, lambda _: len(_ok_syms(_))):
1905 newresult = []
1906 bad_results = []
1907 hit = False
1908 for r in result:
1909 got_s = set()
1910 # update eq with everything that is known so far
1911 eq2 = eq.subs(r)
1912 # if check is True then we see if it satisfies this
1913 # equation, otherwise we just accept it
1914 if check and r:
1915 b = checksol(u, u, eq2, minimal=True)
1916 if b is not None:
1917 # this solution is sufficient to know whether
1918 # it is valid or not so we either accept or
1919 # reject it, then continue
1920 if b:
1921 newresult.append(r)
1922 else:
1923 bad_results.append(r)
1924 continue
1925 # search for a symbol amongst those available that
1926 # can be solved for
1927 ok_syms = _ok_syms(eq2, sort=True)
1928 if not ok_syms:
1929 if r:
1930 newresult.append(r)
1931 break # skip as it's independent of desired symbols
1932 for s in ok_syms:
1933 try:
1934 soln = _vsolve(eq2, s, **flags)
1935 except NotImplementedError:
1936 continue
1937 # put each solution in r and append the now-expanded
1938 # result in the new result list; use copy since the
1939 # solution for s is being added in-place
1940 for sol in soln:
1941 if got_s and any(ss in sol.free_symbols for ss in got_s):
1942 # sol depends on previously solved symbols: discard it
1943 continue
1944 rnew = r.copy()
1945 for k, v in r.items():
1946 rnew[k] = v.subs(s, sol)
1947 # and add this new solution
1948 rnew[s] = sol
1949 # check that it is independent of previous solutions
1950 iset = set(rnew.items())
1951 for i in newresult:
1952 if len(i) < len(iset) and not set(i.items()) - iset:
1953 # this is a superset of a known solution that
1954 # is smaller
1955 break
1956 else:
1957 # keep it
1958 newresult.append(rnew)
1959 hit = True
1960 got_s.add(s)
1961 if not hit:
1962 raise NotImplementedError('could not solve %s' % eq2)
1963 else:
1964 result = newresult
1965 for b in bad_results:
1966 if b in result:
1967 result.remove(b)
1969 if not result:
1970 return False, []
1972 # rely on linear/polynomial system solvers to simplify
1973 # XXX the following tests show that the expressions
1974 # returned are not the same as they would be if simplify
1975 # were applied to this:
1976 # sympy/solvers/ode/tests/test_systems/test__classify_linear_system
1977 # sympy/solvers/tests/test_solvers/test_issue_4886
1978 # so the docs should be updated to reflect that or else
1979 # the following should be `bool(failed) or not linear`
1980 default_simplify = bool(failed)
1981 if flags.get('simplify', default_simplify):
1982 for r in result:
1983 for k in r:
1984 r[k] = simplify(r[k])
1985 flags['simplify'] = False # don't need to do so in checksol now
1987 if checkdens:
1988 result = [r for r in result
1989 if not any(checksol(d, r, **flags) for d in dens)]
1991 if check and not linear:
1992 result = [r for r in result
1993 if not any(checksol(e, r, **flags) is False for e in exprs)]
1995 result = [r for r in result if r]
1996 return linear, result
1999def solve_linear(lhs, rhs=0, symbols=[], exclude=[]):
2000 r"""
2001 Return a tuple derived from ``f = lhs - rhs`` that is one of
2002 the following: ``(0, 1)``, ``(0, 0)``, ``(symbol, solution)``, ``(n, d)``.
2004 Explanation
2005 ===========
2007 ``(0, 1)`` meaning that ``f`` is independent of the symbols in *symbols*
2008 that are not in *exclude*.
2010 ``(0, 0)`` meaning that there is no solution to the equation amongst the
2011 symbols given. If the first element of the tuple is not zero, then the
2012 function is guaranteed to be dependent on a symbol in *symbols*.
2014 ``(symbol, solution)`` where symbol appears linearly in the numerator of
2015 ``f``, is in *symbols* (if given), and is not in *exclude* (if given). No
2016 simplification is done to ``f`` other than a ``mul=True`` expansion, so the
2017 solution will correspond strictly to a unique solution.
2019 ``(n, d)`` where ``n`` and ``d`` are the numerator and denominator of ``f``
2020 when the numerator was not linear in any symbol of interest; ``n`` will
2021 never be a symbol unless a solution for that symbol was found (in which case
2022 the second element is the solution, not the denominator).
2024 Examples
2025 ========
2027 >>> from sympy import cancel, Pow
2029 ``f`` is independent of the symbols in *symbols* that are not in
2030 *exclude*:
2032 >>> from sympy import cos, sin, solve_linear
2033 >>> from sympy.abc import x, y, z
2034 >>> eq = y*cos(x)**2 + y*sin(x)**2 - y # = y*(1 - 1) = 0
2035 >>> solve_linear(eq)
2036 (0, 1)
2037 >>> eq = cos(x)**2 + sin(x)**2 # = 1
2038 >>> solve_linear(eq)
2039 (0, 1)
2040 >>> solve_linear(x, exclude=[x])
2041 (0, 1)
2043 The variable ``x`` appears as a linear variable in each of the
2044 following:
2046 >>> solve_linear(x + y**2)
2047 (x, -y**2)
2048 >>> solve_linear(1/x - y**2)
2049 (x, y**(-2))
2051 When not linear in ``x`` or ``y`` then the numerator and denominator are
2052 returned:
2054 >>> solve_linear(x**2/y**2 - 3)
2055 (x**2 - 3*y**2, y**2)
2057 If the numerator of the expression is a symbol, then ``(0, 0)`` is
2058 returned if the solution for that symbol would have set any
2059 denominator to 0:
2061 >>> eq = 1/(1/x - 2)
2062 >>> eq.as_numer_denom()
2063 (x, 1 - 2*x)
2064 >>> solve_linear(eq)
2065 (0, 0)
2067 But automatic rewriting may cause a symbol in the denominator to
2068 appear in the numerator so a solution will be returned:
2070 >>> (1/x)**-1
2071 x
2072 >>> solve_linear((1/x)**-1)
2073 (x, 0)
2075 Use an unevaluated expression to avoid this:
2077 >>> solve_linear(Pow(1/x, -1, evaluate=False))
2078 (0, 0)
2080 If ``x`` is allowed to cancel in the following expression, then it
2081 appears to be linear in ``x``, but this sort of cancellation is not
2082 done by ``solve_linear`` so the solution will always satisfy the
2083 original expression without causing a division by zero error.
2085 >>> eq = x**2*(1/x - z**2/x)
2086 >>> solve_linear(cancel(eq))
2087 (x, 0)
2088 >>> solve_linear(eq)
2089 (x**2*(1 - z**2), x)
2091 A list of symbols for which a solution is desired may be given:
2093 >>> solve_linear(x + y + z, symbols=[y])
2094 (y, -x - z)
2096 A list of symbols to ignore may also be given:
2098 >>> solve_linear(x + y + z, exclude=[x])
2099 (y, -x - z)
2101 (A solution for ``y`` is obtained because it is the first variable
2102 from the canonically sorted list of symbols that had a linear
2103 solution.)
2105 """
2106 if isinstance(lhs, Eq):
2107 if rhs:
2108 raise ValueError(filldedent('''
2109 If lhs is an Equality, rhs must be 0 but was %s''' % rhs))
2110 rhs = lhs.rhs
2111 lhs = lhs.lhs
2112 dens = None
2113 eq = lhs - rhs
2114 n, d = eq.as_numer_denom()
2115 if not n:
2116 return S.Zero, S.One
2118 free = n.free_symbols
2119 if not symbols:
2120 symbols = free
2121 else:
2122 bad = [s for s in symbols if not s.is_Symbol]
2123 if bad:
2124 if len(bad) == 1:
2125 bad = bad[0]
2126 if len(symbols) == 1:
2127 eg = 'solve(%s, %s)' % (eq, symbols[0])
2128 else:
2129 eg = 'solve(%s, *%s)' % (eq, list(symbols))
2130 raise ValueError(filldedent('''
2131 solve_linear only handles symbols, not %s. To isolate
2132 non-symbols use solve, e.g. >>> %s <<<.
2133 ''' % (bad, eg)))
2134 symbols = free.intersection(symbols)
2135 symbols = symbols.difference(exclude)
2136 if not symbols:
2137 return S.Zero, S.One
2139 # derivatives are easy to do but tricky to analyze to see if they
2140 # are going to disallow a linear solution, so for simplicity we
2141 # just evaluate the ones that have the symbols of interest
2142 derivs = defaultdict(list)
2143 for der in n.atoms(Derivative):
2144 csym = der.free_symbols & symbols
2145 for c in csym:
2146 derivs[c].append(der)
2148 all_zero = True
2149 for xi in sorted(symbols, key=default_sort_key): # canonical order
2150 # if there are derivatives in this var, calculate them now
2151 if isinstance(derivs[xi], list):
2152 derivs[xi] = {der: der.doit() for der in derivs[xi]}
2153 newn = n.subs(derivs[xi])
2154 dnewn_dxi = newn.diff(xi)
2155 # dnewn_dxi can be nonzero if it survives differentation by any
2156 # of its free symbols
2157 free = dnewn_dxi.free_symbols
2158 if dnewn_dxi and (not free or any(dnewn_dxi.diff(s) for s in free) or free == symbols):
2159 all_zero = False
2160 if dnewn_dxi is S.NaN:
2161 break
2162 if xi not in dnewn_dxi.free_symbols:
2163 vi = -1/dnewn_dxi*(newn.subs(xi, 0))
2164 if dens is None:
2165 dens = _simple_dens(eq, symbols)
2166 if not any(checksol(di, {xi: vi}, minimal=True) is True
2167 for di in dens):
2168 # simplify any trivial integral
2169 irep = [(i, i.doit()) for i in vi.atoms(Integral) if
2170 i.function.is_number]
2171 # do a slight bit of simplification
2172 vi = expand_mul(vi.subs(irep))
2173 return xi, vi
2174 if all_zero:
2175 return S.Zero, S.One
2176 if n.is_Symbol: # no solution for this symbol was found
2177 return S.Zero, S.Zero
2178 return n, d
2181def minsolve_linear_system(system, *symbols, **flags):
2182 r"""
2183 Find a particular solution to a linear system.
2185 Explanation
2186 ===========
2188 In particular, try to find a solution with the minimal possible number
2189 of non-zero variables using a naive algorithm with exponential complexity.
2190 If ``quick=True``, a heuristic is used.
2192 """
2193 quick = flags.get('quick', False)
2194 # Check if there are any non-zero solutions at all
2195 s0 = solve_linear_system(system, *symbols, **flags)
2196 if not s0 or all(v == 0 for v in s0.values()):
2197 return s0
2198 if quick:
2199 # We just solve the system and try to heuristically find a nice
2200 # solution.
2201 s = solve_linear_system(system, *symbols)
2202 def update(determined, solution):
2203 delete = []
2204 for k, v in solution.items():
2205 solution[k] = v.subs(determined)
2206 if not solution[k].free_symbols:
2207 delete.append(k)
2208 determined[k] = solution[k]
2209 for k in delete:
2210 del solution[k]
2211 determined = {}
2212 update(determined, s)
2213 while s:
2214 # NOTE sort by default_sort_key to get deterministic result
2215 k = max((k for k in s.values()),
2216 key=lambda x: (len(x.free_symbols), default_sort_key(x)))
2217 kfree = k.free_symbols
2218 x = next(reversed(list(ordered(kfree))))
2219 if len(kfree) != 1:
2220 determined[x] = S.Zero
2221 else:
2222 val = _vsolve(k, x, check=False)[0]
2223 if not val and not any(v.subs(x, val) for v in s.values()):
2224 determined[x] = S.One
2225 else:
2226 determined[x] = val
2227 update(determined, s)
2228 return determined
2229 else:
2230 # We try to select n variables which we want to be non-zero.
2231 # All others will be assumed zero. We try to solve the modified system.
2232 # If there is a non-trivial solution, just set the free variables to
2233 # one. If we do this for increasing n, trying all combinations of
2234 # variables, we will find an optimal solution.
2235 # We speed up slightly by starting at one less than the number of
2236 # variables the quick method manages.
2237 N = len(symbols)
2238 bestsol = minsolve_linear_system(system, *symbols, quick=True)
2239 n0 = len([x for x in bestsol.values() if x != 0])
2240 for n in range(n0 - 1, 1, -1):
2241 debugf('minsolve: %s', n)
2242 thissol = None
2243 for nonzeros in combinations(range(N), n):
2244 subm = Matrix([system.col(i).T for i in nonzeros] + [system.col(-1).T]).T
2245 s = solve_linear_system(subm, *[symbols[i] for i in nonzeros])
2246 if s and not all(v == 0 for v in s.values()):
2247 subs = [(symbols[v], S.One) for v in nonzeros]
2248 for k, v in s.items():
2249 s[k] = v.subs(subs)
2250 for sym in symbols:
2251 if sym not in s:
2252 if symbols.index(sym) in nonzeros:
2253 s[sym] = S.One
2254 else:
2255 s[sym] = S.Zero
2256 thissol = s
2257 break
2258 if thissol is None:
2259 break
2260 bestsol = thissol
2261 return bestsol
2264def solve_linear_system(system, *symbols, **flags):
2265 r"""
2266 Solve system of $N$ linear equations with $M$ variables, which means
2267 both under- and overdetermined systems are supported.
2269 Explanation
2270 ===========
2272 The possible number of solutions is zero, one, or infinite. Respectively,
2273 this procedure will return None or a dictionary with solutions. In the
2274 case of underdetermined systems, all arbitrary parameters are skipped.
2275 This may cause a situation in which an empty dictionary is returned.
2276 In that case, all symbols can be assigned arbitrary values.
2278 Input to this function is a $N\times M + 1$ matrix, which means it has
2279 to be in augmented form. If you prefer to enter $N$ equations and $M$
2280 unknowns then use ``solve(Neqs, *Msymbols)`` instead. Note: a local
2281 copy of the matrix is made by this routine so the matrix that is
2282 passed will not be modified.
2284 The algorithm used here is fraction-free Gaussian elimination,
2285 which results, after elimination, in an upper-triangular matrix.
2286 Then solutions are found using back-substitution. This approach
2287 is more efficient and compact than the Gauss-Jordan method.
2289 Examples
2290 ========
2292 >>> from sympy import Matrix, solve_linear_system
2293 >>> from sympy.abc import x, y
2295 Solve the following system::
2297 x + 4 y == 2
2298 -2 x + y == 14
2300 >>> system = Matrix(( (1, 4, 2), (-2, 1, 14)))
2301 >>> solve_linear_system(system, x, y)
2302 {x: -6, y: 2}
2304 A degenerate system returns an empty dictionary:
2306 >>> system = Matrix(( (0,0,0), (0,0,0) ))
2307 >>> solve_linear_system(system, x, y)
2308 {}
2310 """
2311 assert system.shape[1] == len(symbols) + 1
2313 # This is just a wrapper for solve_lin_sys
2314 eqs = list(system * Matrix(symbols + (-1,)))
2315 eqs, ring = sympy_eqs_to_ring(eqs, symbols)
2316 sol = solve_lin_sys(eqs, ring, _raw=False)
2317 if sol is not None:
2318 sol = {sym:val for sym, val in sol.items() if sym != val}
2319 return sol
2322def solve_undetermined_coeffs(equ, coeffs, *syms, **flags):
2323 r"""
2324 Solve a system of equations in $k$ parameters that is formed by
2325 matching coefficients in variables ``coeffs`` that are on
2326 factors dependent on the remaining variables (or those given
2327 explicitly by ``syms``.
2329 Explanation
2330 ===========
2332 The result of this function is a dictionary with symbolic values of those
2333 parameters with respect to coefficients in $q$ -- empty if there
2334 is no solution or coefficients do not appear in the equation -- else
2335 None (if the system was not recognized). If there is more than one
2336 solution, the solutions are passed as a list. The output can be modified using
2337 the same semantics as for `solve` since the flags that are passed are sent
2338 directly to `solve` so, for example the flag ``dict=True`` will always return a list
2339 of solutions as dictionaries.
2341 This function accepts both Equality and Expr class instances.
2342 The solving process is most efficient when symbols are specified
2343 in addition to parameters to be determined, but an attempt to
2344 determine them (if absent) will be made. If an expected solution is not
2345 obtained (and symbols were not specified) try specifying them.
2347 Examples
2348 ========
2350 >>> from sympy import Eq, solve_undetermined_coeffs
2351 >>> from sympy.abc import a, b, c, h, p, k, x, y
2353 >>> solve_undetermined_coeffs(Eq(a*x + a + b, x/2), [a, b], x)
2354 {a: 1/2, b: -1/2}
2355 >>> solve_undetermined_coeffs(a - 2, [a])
2356 {a: 2}
2358 The equation can be nonlinear in the symbols:
2360 >>> X, Y, Z = y, x**y, y*x**y
2361 >>> eq = a*X + b*Y + c*Z - X - 2*Y - 3*Z
2362 >>> coeffs = a, b, c
2363 >>> syms = x, y
2364 >>> solve_undetermined_coeffs(eq, coeffs, syms)
2365 {a: 1, b: 2, c: 3}
2367 And the system can be nonlinear in coefficients, too, but if
2368 there is only a single solution, it will be returned as a
2369 dictionary:
2371 >>> eq = a*x**2 + b*x + c - ((x - h)**2 + 4*p*k)/4/p
2372 >>> solve_undetermined_coeffs(eq, (h, p, k), x)
2373 {h: -b/(2*a), k: (4*a*c - b**2)/(4*a), p: 1/(4*a)}
2375 Multiple solutions are always returned in a list:
2377 >>> solve_undetermined_coeffs(a**2*x + b - x, [a, b], x)
2378 [{a: -1, b: 0}, {a: 1, b: 0}]
2380 Using flag ``dict=True`` (in keeping with semantics in :func:`~.solve`)
2381 will force the result to always be a list with any solutions
2382 as elements in that list.
2384 >>> solve_undetermined_coeffs(a*x - 2*x, [a], dict=True)
2385 [{a: 2}]
2386 """
2387 if not (coeffs and all(i.is_Symbol for i in coeffs)):
2388 raise ValueError('must provide symbols for coeffs')
2390 if isinstance(equ, Eq):
2391 eq = equ.lhs - equ.rhs
2392 else:
2393 eq = equ
2395 ceq = cancel(eq)
2396 xeq = _mexpand(ceq.as_numer_denom()[0], recursive=True)
2398 free = xeq.free_symbols
2399 coeffs = free & set(coeffs)
2400 if not coeffs:
2401 return ([], {}) if flags.get('set', None) else [] # solve(0, x) -> []
2403 if not syms:
2404 # e.g. A*exp(x) + B - (exp(x) + y) separated into parts that
2405 # don't/do depend on coeffs gives
2406 # -(exp(x) + y), A*exp(x) + B
2407 # then see what symbols are common to both
2408 # {x} = {x, A, B} - {x, y}
2409 ind, dep = xeq.as_independent(*coeffs, as_Add=True)
2410 dfree = dep.free_symbols
2411 syms = dfree & ind.free_symbols
2412 if not syms:
2413 # but if the system looks like (a + b)*x + b - c
2414 # then {} = {a, b, x} - c
2415 # so calculate {x} = {a, b, x} - {a, b}
2416 syms = dfree - set(coeffs)
2417 if not syms:
2418 syms = [Dummy()]
2419 else:
2420 if len(syms) == 1 and iterable(syms[0]):
2421 syms = syms[0]
2422 e, s, _ = recast_to_symbols([xeq], syms)
2423 xeq = e[0]
2424 syms = s
2426 # find the functional forms in which symbols appear
2428 gens = set(xeq.as_coefficients_dict(*syms).keys()) - {1}
2429 cset = set(coeffs)
2430 if any(g.has_xfree(cset) for g in gens):
2431 return # a generator contained a coefficient symbol
2433 # make sure we are working with symbols for generators
2435 e, gens, _ = recast_to_symbols([xeq], list(gens))
2436 xeq = e[0]
2438 # collect coefficients in front of generators
2440 system = list(collect(xeq, gens, evaluate=False).values())
2442 # get a solution
2444 soln = solve(system, coeffs, **flags)
2446 # unpack unless told otherwise if length is 1
2448 settings = flags.get('dict', None) or flags.get('set', None)
2449 if type(soln) is dict or settings or len(soln) != 1:
2450 return soln
2451 return soln[0]
2454def solve_linear_system_LU(matrix, syms):
2455 """
2456 Solves the augmented matrix system using ``LUsolve`` and returns a
2457 dictionary in which solutions are keyed to the symbols of *syms* as ordered.
2459 Explanation
2460 ===========
2462 The matrix must be invertible.
2464 Examples
2465 ========
2467 >>> from sympy import Matrix, solve_linear_system_LU
2468 >>> from sympy.abc import x, y, z
2470 >>> solve_linear_system_LU(Matrix([
2471 ... [1, 2, 0, 1],
2472 ... [3, 2, 2, 1],
2473 ... [2, 0, 0, 1]]), [x, y, z])
2474 {x: 1/2, y: 1/4, z: -1/2}
2476 See Also
2477 ========
2479 LUsolve
2481 """
2482 if matrix.rows != matrix.cols - 1:
2483 raise ValueError("Rows should be equal to columns - 1")
2484 A = matrix[:matrix.rows, :matrix.rows]
2485 b = matrix[:, matrix.cols - 1:]
2486 soln = A.LUsolve(b)
2487 solutions = {}
2488 for i in range(soln.rows):
2489 solutions[syms[i]] = soln[i, 0]
2490 return solutions
2493def det_perm(M):
2494 """
2495 Return the determinant of *M* by using permutations to select factors.
2497 Explanation
2498 ===========
2500 For sizes larger than 8 the number of permutations becomes prohibitively
2501 large, or if there are no symbols in the matrix, it is better to use the
2502 standard determinant routines (e.g., ``M.det()``.)
2504 See Also
2505 ========
2507 det_minor
2508 det_quick
2510 """
2511 args = []
2512 s = True
2513 n = M.rows
2514 list_ = M.flat()
2515 for perm in generate_bell(n):
2516 fac = []
2517 idx = 0
2518 for j in perm:
2519 fac.append(list_[idx + j])
2520 idx += n
2521 term = Mul(*fac) # disaster with unevaluated Mul -- takes forever for n=7
2522 args.append(term if s else -term)
2523 s = not s
2524 return Add(*args)
2527def det_minor(M):
2528 """
2529 Return the ``det(M)`` computed from minors without
2530 introducing new nesting in products.
2532 See Also
2533 ========
2535 det_perm
2536 det_quick
2538 """
2539 n = M.rows
2540 if n == 2:
2541 return M[0, 0]*M[1, 1] - M[1, 0]*M[0, 1]
2542 else:
2543 return sum([(1, -1)[i % 2]*Add(*[M[0, i]*d for d in
2544 Add.make_args(det_minor(M.minor_submatrix(0, i)))])
2545 if M[0, i] else S.Zero for i in range(n)])
2548def det_quick(M, method=None):
2549 """
2550 Return ``det(M)`` assuming that either
2551 there are lots of zeros or the size of the matrix
2552 is small. If this assumption is not met, then the normal
2553 Matrix.det function will be used with method = ``method``.
2555 See Also
2556 ========
2558 det_minor
2559 det_perm
2561 """
2562 if any(i.has(Symbol) for i in M):
2563 if M.rows < 8 and all(i.has(Symbol) for i in M):
2564 return det_perm(M)
2565 return det_minor(M)
2566 else:
2567 return M.det(method=method) if method else M.det()
2570def inv_quick(M):
2571 """Return the inverse of ``M``, assuming that either
2572 there are lots of zeros or the size of the matrix
2573 is small.
2574 """
2575 if not all(i.is_Number for i in M):
2576 if not any(i.is_Number for i in M):
2577 det = lambda _: det_perm(_)
2578 else:
2579 det = lambda _: det_minor(_)
2580 else:
2581 return M.inv()
2582 n = M.rows
2583 d = det(M)
2584 if d == S.Zero:
2585 raise NonInvertibleMatrixError("Matrix det == 0; not invertible")
2586 ret = zeros(n)
2587 s1 = -1
2588 for i in range(n):
2589 s = s1 = -s1
2590 for j in range(n):
2591 di = det(M.minor_submatrix(i, j))
2592 ret[j, i] = s*di/d
2593 s = -s
2594 return ret
2597# these are functions that have multiple inverse values per period
2598multi_inverses = {
2599 sin: lambda x: (asin(x), S.Pi - asin(x)),
2600 cos: lambda x: (acos(x), 2*S.Pi - acos(x)),
2601}
2604def _vsolve(e, s, **flags):
2605 """return list of scalar values for the solution of e for symbol s"""
2606 return [i[s] for i in _solve(e, s, **flags)]
2609def _tsolve(eq, sym, **flags):
2610 """
2611 Helper for ``_solve`` that solves a transcendental equation with respect
2612 to the given symbol. Various equations containing powers and logarithms,
2613 can be solved.
2615 There is currently no guarantee that all solutions will be returned or
2616 that a real solution will be favored over a complex one.
2618 Either a list of potential solutions will be returned or None will be
2619 returned (in the case that no method was known to get a solution
2620 for the equation). All other errors (like the inability to cast an
2621 expression as a Poly) are unhandled.
2623 Examples
2624 ========
2626 >>> from sympy import log, ordered
2627 >>> from sympy.solvers.solvers import _tsolve as tsolve
2628 >>> from sympy.abc import x
2630 >>> list(ordered(tsolve(3**(2*x + 5) - 4, x)))
2631 [-5/2 + log(2)/log(3), (-5*log(3)/2 + log(2) + I*pi)/log(3)]
2633 >>> tsolve(log(x) + 2*x, x)
2634 [LambertW(2)/2]
2636 """
2637 if 'tsolve_saw' not in flags:
2638 flags['tsolve_saw'] = []
2639 if eq in flags['tsolve_saw']:
2640 return None
2641 else:
2642 flags['tsolve_saw'].append(eq)
2644 rhs, lhs = _invert(eq, sym)
2646 if lhs == sym:
2647 return [rhs]
2648 try:
2649 if lhs.is_Add:
2650 # it's time to try factoring; powdenest is used
2651 # to try get powers in standard form for better factoring
2652 f = factor(powdenest(lhs - rhs))
2653 if f.is_Mul:
2654 return _vsolve(f, sym, **flags)
2655 if rhs:
2656 f = logcombine(lhs, force=flags.get('force', True))
2657 if f.count(log) != lhs.count(log):
2658 if isinstance(f, log):
2659 return _vsolve(f.args[0] - exp(rhs), sym, **flags)
2660 return _tsolve(f - rhs, sym, **flags)
2662 elif lhs.is_Pow:
2663 if lhs.exp.is_Integer:
2664 if lhs - rhs != eq:
2665 return _vsolve(lhs - rhs, sym, **flags)
2667 if sym not in lhs.exp.free_symbols:
2668 return _vsolve(lhs.base - rhs**(1/lhs.exp), sym, **flags)
2670 # _tsolve calls this with Dummy before passing the actual number in.
2671 if any(t.is_Dummy for t in rhs.free_symbols):
2672 raise NotImplementedError # _tsolve will call here again...
2674 # a ** g(x) == 0
2675 if not rhs:
2676 # f(x)**g(x) only has solutions where f(x) == 0 and g(x) != 0 at
2677 # the same place
2678 sol_base = _vsolve(lhs.base, sym, **flags)
2679 return [s for s in sol_base if lhs.exp.subs(sym, s) != 0] # XXX use checksol here?
2681 # a ** g(x) == b
2682 if not lhs.base.has(sym):
2683 if lhs.base == 0:
2684 return _vsolve(lhs.exp, sym, **flags) if rhs != 0 else []
2686 # Gets most solutions...
2687 if lhs.base == rhs.as_base_exp()[0]:
2688 # handles case when bases are equal
2689 sol = _vsolve(lhs.exp - rhs.as_base_exp()[1], sym, **flags)
2690 else:
2691 # handles cases when bases are not equal and exp
2692 # may or may not be equal
2693 f = exp(log(lhs.base)*lhs.exp) - exp(log(rhs))
2694 sol = _vsolve(f, sym, **flags)
2696 # Check for duplicate solutions
2697 def equal(expr1, expr2):
2698 _ = Dummy()
2699 eq = checksol(expr1 - _, _, expr2)
2700 if eq is None:
2701 if nsimplify(expr1) != nsimplify(expr2):
2702 return False
2703 # they might be coincidentally the same
2704 # so check more rigorously
2705 eq = expr1.equals(expr2) # XXX expensive but necessary?
2706 return eq
2708 # Guess a rational exponent
2709 e_rat = nsimplify(log(abs(rhs))/log(abs(lhs.base)))
2710 e_rat = simplify(posify(e_rat)[0])
2711 n, d = fraction(e_rat)
2712 if expand(lhs.base**n - rhs**d) == 0:
2713 sol = [s for s in sol if not equal(lhs.exp.subs(sym, s), e_rat)]
2714 sol.extend(_vsolve(lhs.exp - e_rat, sym, **flags))
2716 return list(set(sol))
2718 # f(x) ** g(x) == c
2719 else:
2720 sol = []
2721 logform = lhs.exp*log(lhs.base) - log(rhs)
2722 if logform != lhs - rhs:
2723 try:
2724 sol.extend(_vsolve(logform, sym, **flags))
2725 except NotImplementedError:
2726 pass
2728 # Collect possible solutions and check with substitution later.
2729 check = []
2730 if rhs == 1:
2731 # f(x) ** g(x) = 1 -- g(x)=0 or f(x)=+-1
2732 check.extend(_vsolve(lhs.exp, sym, **flags))
2733 check.extend(_vsolve(lhs.base - 1, sym, **flags))
2734 check.extend(_vsolve(lhs.base + 1, sym, **flags))
2735 elif rhs.is_Rational:
2736 for d in (i for i in divisors(abs(rhs.p)) if i != 1):
2737 e, t = integer_log(rhs.p, d)
2738 if not t:
2739 continue # rhs.p != d**b
2740 for s in divisors(abs(rhs.q)):
2741 if s**e== rhs.q:
2742 r = Rational(d, s)
2743 check.extend(_vsolve(lhs.base - r, sym, **flags))
2744 check.extend(_vsolve(lhs.base + r, sym, **flags))
2745 check.extend(_vsolve(lhs.exp - e, sym, **flags))
2746 elif rhs.is_irrational:
2747 b_l, e_l = lhs.base.as_base_exp()
2748 n, d = (e_l*lhs.exp).as_numer_denom()
2749 b, e = sqrtdenest(rhs).as_base_exp()
2750 check = [sqrtdenest(i) for i in (_vsolve(lhs.base - b, sym, **flags))]
2751 check.extend([sqrtdenest(i) for i in (_vsolve(lhs.exp - e, sym, **flags))])
2752 if e_l*d != 1:
2753 check.extend(_vsolve(b_l**n - rhs**(e_l*d), sym, **flags))
2754 for s in check:
2755 ok = checksol(eq, sym, s)
2756 if ok is None:
2757 ok = eq.subs(sym, s).equals(0)
2758 if ok:
2759 sol.append(s)
2760 return list(set(sol))
2762 elif lhs.is_Function and len(lhs.args) == 1:
2763 if lhs.func in multi_inverses:
2764 # sin(x) = 1/3 -> x - asin(1/3) & x - (pi - asin(1/3))
2765 soln = []
2766 for i in multi_inverses[type(lhs)](rhs):
2767 soln.extend(_vsolve(lhs.args[0] - i, sym, **flags))
2768 return list(set(soln))
2769 elif lhs.func == LambertW:
2770 return _vsolve(lhs.args[0] - rhs*exp(rhs), sym, **flags)
2772 rewrite = lhs.rewrite(exp)
2773 if rewrite != lhs:
2774 return _vsolve(rewrite - rhs, sym, **flags)
2775 except NotImplementedError:
2776 pass
2778 # maybe it is a lambert pattern
2779 if flags.pop('bivariate', True):
2780 # lambert forms may need some help being recognized, e.g. changing
2781 # 2**(3*x) + x**3*log(2)**3 + 3*x**2*log(2)**2 + 3*x*log(2) + 1
2782 # to 2**(3*x) + (x*log(2) + 1)**3
2784 # make generator in log have exponent of 1
2785 logs = eq.atoms(log)
2786 spow = min(
2787 {i.exp for j in logs for i in j.atoms(Pow)
2788 if i.base == sym} or {1})
2789 if spow != 1:
2790 p = sym**spow
2791 u = Dummy('bivariate-cov')
2792 ueq = eq.subs(p, u)
2793 if not ueq.has_free(sym):
2794 sol = _vsolve(ueq, u, **flags)
2795 inv = _vsolve(p - u, sym)
2796 return [i.subs(u, s) for i in inv for s in sol]
2798 g = _filtered_gens(eq.as_poly(), sym)
2799 up_or_log = set()
2800 for gi in g:
2801 if isinstance(gi, (exp, log)) or (gi.is_Pow and gi.base == S.Exp1):
2802 up_or_log.add(gi)
2803 elif gi.is_Pow:
2804 gisimp = powdenest(expand_power_exp(gi))
2805 if gisimp.is_Pow and sym in gisimp.exp.free_symbols:
2806 up_or_log.add(gi)
2807 eq_down = expand_log(expand_power_exp(eq)).subs(
2808 dict(list(zip(up_or_log, [0]*len(up_or_log)))))
2809 eq = expand_power_exp(factor(eq_down, deep=True) + (eq - eq_down))
2810 rhs, lhs = _invert(eq, sym)
2811 if lhs.has(sym):
2812 try:
2813 poly = lhs.as_poly()
2814 g = _filtered_gens(poly, sym)
2815 _eq = lhs - rhs
2816 sols = _solve_lambert(_eq, sym, g)
2817 # use a simplified form if it satisfies eq
2818 # and has fewer operations
2819 for n, s in enumerate(sols):
2820 ns = nsimplify(s)
2821 if ns != s and ns.count_ops() <= s.count_ops():
2822 ok = checksol(_eq, sym, ns)
2823 if ok is None:
2824 ok = _eq.subs(sym, ns).equals(0)
2825 if ok:
2826 sols[n] = ns
2827 return sols
2828 except NotImplementedError:
2829 # maybe it's a convoluted function
2830 if len(g) == 2:
2831 try:
2832 gpu = bivariate_type(lhs - rhs, *g)
2833 if gpu is None:
2834 raise NotImplementedError
2835 g, p, u = gpu
2836 flags['bivariate'] = False
2837 inversion = _tsolve(g - u, sym, **flags)
2838 if inversion:
2839 sol = _vsolve(p, u, **flags)
2840 return list({i.subs(u, s)
2841 for i in inversion for s in sol})
2842 except NotImplementedError:
2843 pass
2844 else:
2845 pass
2847 if flags.pop('force', True):
2848 flags['force'] = False
2849 pos, reps = posify(lhs - rhs)
2850 if rhs == S.ComplexInfinity:
2851 return []
2852 for u, s in reps.items():
2853 if s == sym:
2854 break
2855 else:
2856 u = sym
2857 if pos.has(u):
2858 try:
2859 soln = _vsolve(pos, u, **flags)
2860 return [s.subs(reps) for s in soln]
2861 except NotImplementedError:
2862 pass
2863 else:
2864 pass # here for coverage
2866 return # here for coverage
2869# TODO: option for calculating J numerically
2871@conserve_mpmath_dps
2872def nsolve(*args, dict=False, **kwargs):
2873 r"""
2874 Solve a nonlinear equation system numerically: ``nsolve(f, [args,] x0,
2875 modules=['mpmath'], **kwargs)``.
2877 Explanation
2878 ===========
2880 ``f`` is a vector function of symbolic expressions representing the system.
2881 *args* are the variables. If there is only one variable, this argument can
2882 be omitted. ``x0`` is a starting vector close to a solution.
2884 Use the modules keyword to specify which modules should be used to
2885 evaluate the function and the Jacobian matrix. Make sure to use a module
2886 that supports matrices. For more information on the syntax, please see the
2887 docstring of ``lambdify``.
2889 If the keyword arguments contain ``dict=True`` (default is False) ``nsolve``
2890 will return a list (perhaps empty) of solution mappings. This might be
2891 especially useful if you want to use ``nsolve`` as a fallback to solve since
2892 using the dict argument for both methods produces return values of
2893 consistent type structure. Please note: to keep this consistent with
2894 ``solve``, the solution will be returned in a list even though ``nsolve``
2895 (currently at least) only finds one solution at a time.
2897 Overdetermined systems are supported.
2899 Examples
2900 ========
2902 >>> from sympy import Symbol, nsolve
2903 >>> import mpmath
2904 >>> mpmath.mp.dps = 15
2905 >>> x1 = Symbol('x1')
2906 >>> x2 = Symbol('x2')
2907 >>> f1 = 3 * x1**2 - 2 * x2**2 - 1
2908 >>> f2 = x1**2 - 2 * x1 + x2**2 + 2 * x2 - 8
2909 >>> print(nsolve((f1, f2), (x1, x2), (-1, 1)))
2910 Matrix([[-1.19287309935246], [1.27844411169911]])
2912 For one-dimensional functions the syntax is simplified:
2914 >>> from sympy import sin, nsolve
2915 >>> from sympy.abc import x
2916 >>> nsolve(sin(x), x, 2)
2917 3.14159265358979
2918 >>> nsolve(sin(x), 2)
2919 3.14159265358979
2921 To solve with higher precision than the default, use the prec argument:
2923 >>> from sympy import cos
2924 >>> nsolve(cos(x) - x, 1)
2925 0.739085133215161
2926 >>> nsolve(cos(x) - x, 1, prec=50)
2927 0.73908513321516064165531208767387340401341175890076
2928 >>> cos(_)
2929 0.73908513321516064165531208767387340401341175890076
2931 To solve for complex roots of real functions, a nonreal initial point
2932 must be specified:
2934 >>> from sympy import I
2935 >>> nsolve(x**2 + 2, I)
2936 1.4142135623731*I
2938 ``mpmath.findroot`` is used and you can find their more extensive
2939 documentation, especially concerning keyword parameters and
2940 available solvers. Note, however, that functions which are very
2941 steep near the root, the verification of the solution may fail. In
2942 this case you should use the flag ``verify=False`` and
2943 independently verify the solution.
2945 >>> from sympy import cos, cosh
2946 >>> f = cos(x)*cosh(x) - 1
2947 >>> nsolve(f, 3.14*100)
2948 Traceback (most recent call last):
2949 ...
2950 ValueError: Could not find root within given tolerance. (1.39267e+230 > 2.1684e-19)
2951 >>> ans = nsolve(f, 3.14*100, verify=False); ans
2952 312.588469032184
2953 >>> f.subs(x, ans).n(2)
2954 2.1e+121
2955 >>> (f/f.diff(x)).subs(x, ans).n(2)
2956 7.4e-15
2958 One might safely skip the verification if bounds of the root are known
2959 and a bisection method is used:
2961 >>> bounds = lambda i: (3.14*i, 3.14*(i + 1))
2962 >>> nsolve(f, bounds(100), solver='bisect', verify=False)
2963 315.730061685774
2965 Alternatively, a function may be better behaved when the
2966 denominator is ignored. Since this is not always the case, however,
2967 the decision of what function to use is left to the discretion of
2968 the user.
2970 >>> eq = x**2/(1 - x)/(1 - 2*x)**2 - 100
2971 >>> nsolve(eq, 0.46)
2972 Traceback (most recent call last):
2973 ...
2974 ValueError: Could not find root within given tolerance. (10000 > 2.1684e-19)
2975 Try another starting point or tweak arguments.
2976 >>> nsolve(eq.as_numer_denom()[0], 0.46)
2977 0.46792545969349058
2979 """
2980 # there are several other SymPy functions that use method= so
2981 # guard against that here
2982 if 'method' in kwargs:
2983 raise ValueError(filldedent('''
2984 Keyword "method" should not be used in this context. When using
2985 some mpmath solvers directly, the keyword "method" is
2986 used, but when using nsolve (and findroot) the keyword to use is
2987 "solver".'''))
2989 if 'prec' in kwargs:
2990 import mpmath
2991 mpmath.mp.dps = kwargs.pop('prec')
2993 # keyword argument to return result as a dictionary
2994 as_dict = dict
2995 from builtins import dict # to unhide the builtin
2997 # interpret arguments
2998 if len(args) == 3:
2999 f = args[0]
3000 fargs = args[1]
3001 x0 = args[2]
3002 if iterable(fargs) and iterable(x0):
3003 if len(x0) != len(fargs):
3004 raise TypeError('nsolve expected exactly %i guess vectors, got %i'
3005 % (len(fargs), len(x0)))
3006 elif len(args) == 2:
3007 f = args[0]
3008 fargs = None
3009 x0 = args[1]
3010 if iterable(f):
3011 raise TypeError('nsolve expected 3 arguments, got 2')
3012 elif len(args) < 2:
3013 raise TypeError('nsolve expected at least 2 arguments, got %i'
3014 % len(args))
3015 else:
3016 raise TypeError('nsolve expected at most 3 arguments, got %i'
3017 % len(args))
3018 modules = kwargs.get('modules', ['mpmath'])
3019 if iterable(f):
3020 f = list(f)
3021 for i, fi in enumerate(f):
3022 if isinstance(fi, Eq):
3023 f[i] = fi.lhs - fi.rhs
3024 f = Matrix(f).T
3025 if iterable(x0):
3026 x0 = list(x0)
3027 if not isinstance(f, Matrix):
3028 # assume it's a SymPy expression
3029 if isinstance(f, Eq):
3030 f = f.lhs - f.rhs
3031 syms = f.free_symbols
3032 if fargs is None:
3033 fargs = syms.copy().pop()
3034 if not (len(syms) == 1 and (fargs in syms or fargs[0] in syms)):
3035 raise ValueError(filldedent('''
3036 expected a one-dimensional and numerical function'''))
3038 # the function is much better behaved if there is no denominator
3039 # but sending the numerator is left to the user since sometimes
3040 # the function is better behaved when the denominator is present
3041 # e.g., issue 11768
3043 f = lambdify(fargs, f, modules)
3044 x = sympify(findroot(f, x0, **kwargs))
3045 if as_dict:
3046 return [{fargs: x}]
3047 return x
3049 if len(fargs) > f.cols:
3050 raise NotImplementedError(filldedent('''
3051 need at least as many equations as variables'''))
3052 verbose = kwargs.get('verbose', False)
3053 if verbose:
3054 print('f(x):')
3055 print(f)
3056 # derive Jacobian
3057 J = f.jacobian(fargs)
3058 if verbose:
3059 print('J(x):')
3060 print(J)
3061 # create functions
3062 f = lambdify(fargs, f.T, modules)
3063 J = lambdify(fargs, J, modules)
3064 # solve the system numerically
3065 x = findroot(f, x0, J=J, **kwargs)
3066 if as_dict:
3067 return [dict(zip(fargs, [sympify(xi) for xi in x]))]
3068 return Matrix(x)
3071def _invert(eq, *symbols, **kwargs):
3072 """
3073 Return tuple (i, d) where ``i`` is independent of *symbols* and ``d``
3074 contains symbols.
3076 Explanation
3077 ===========
3079 ``i`` and ``d`` are obtained after recursively using algebraic inversion
3080 until an uninvertible ``d`` remains. If there are no free symbols then
3081 ``d`` will be zero. Some (but not necessarily all) solutions to the
3082 expression ``i - d`` will be related to the solutions of the original
3083 expression.
3085 Examples
3086 ========
3088 >>> from sympy.solvers.solvers import _invert as invert
3089 >>> from sympy import sqrt, cos
3090 >>> from sympy.abc import x, y
3091 >>> invert(x - 3)
3092 (3, x)
3093 >>> invert(3)
3094 (3, 0)
3095 >>> invert(2*cos(x) - 1)
3096 (1/2, cos(x))
3097 >>> invert(sqrt(x) - 3)
3098 (3, sqrt(x))
3099 >>> invert(sqrt(x) + y, x)
3100 (-y, sqrt(x))
3101 >>> invert(sqrt(x) + y, y)
3102 (-sqrt(x), y)
3103 >>> invert(sqrt(x) + y, x, y)
3104 (0, sqrt(x) + y)
3106 If there is more than one symbol in a power's base and the exponent
3107 is not an Integer, then the principal root will be used for the
3108 inversion:
3110 >>> invert(sqrt(x + y) - 2)
3111 (4, x + y)
3112 >>> invert(sqrt(x + y) - 2)
3113 (4, x + y)
3115 If the exponent is an Integer, setting ``integer_power`` to True
3116 will force the principal root to be selected:
3118 >>> invert(x**2 - 4, integer_power=True)
3119 (2, x)
3121 """
3122 eq = sympify(eq)
3123 if eq.args:
3124 # make sure we are working with flat eq
3125 eq = eq.func(*eq.args)
3126 free = eq.free_symbols
3127 if not symbols:
3128 symbols = free
3129 if not free & set(symbols):
3130 return eq, S.Zero
3132 dointpow = bool(kwargs.get('integer_power', False))
3134 lhs = eq
3135 rhs = S.Zero
3136 while True:
3137 was = lhs
3138 while True:
3139 indep, dep = lhs.as_independent(*symbols)
3141 # dep + indep == rhs
3142 if lhs.is_Add:
3143 # this indicates we have done it all
3144 if indep.is_zero:
3145 break
3147 lhs = dep
3148 rhs -= indep
3150 # dep * indep == rhs
3151 else:
3152 # this indicates we have done it all
3153 if indep is S.One:
3154 break
3156 lhs = dep
3157 rhs /= indep
3159 # collect like-terms in symbols
3160 if lhs.is_Add:
3161 terms = {}
3162 for a in lhs.args:
3163 i, d = a.as_independent(*symbols)
3164 terms.setdefault(d, []).append(i)
3165 if any(len(v) > 1 for v in terms.values()):
3166 args = []
3167 for d, i in terms.items():
3168 if len(i) > 1:
3169 args.append(Add(*i)*d)
3170 else:
3171 args.append(i[0]*d)
3172 lhs = Add(*args)
3174 # if it's a two-term Add with rhs = 0 and two powers we can get the
3175 # dependent terms together, e.g. 3*f(x) + 2*g(x) -> f(x)/g(x) = -2/3
3176 if lhs.is_Add and not rhs and len(lhs.args) == 2 and \
3177 not lhs.is_polynomial(*symbols):
3178 a, b = ordered(lhs.args)
3179 ai, ad = a.as_independent(*symbols)
3180 bi, bd = b.as_independent(*symbols)
3181 if any(_ispow(i) for i in (ad, bd)):
3182 a_base, a_exp = ad.as_base_exp()
3183 b_base, b_exp = bd.as_base_exp()
3184 if a_base == b_base:
3185 # a = -b
3186 lhs = powsimp(powdenest(ad/bd))
3187 rhs = -bi/ai
3188 else:
3189 rat = ad/bd
3190 _lhs = powsimp(ad/bd)
3191 if _lhs != rat:
3192 lhs = _lhs
3193 rhs = -bi/ai
3194 elif ai == -bi:
3195 if isinstance(ad, Function) and ad.func == bd.func:
3196 if len(ad.args) == len(bd.args) == 1:
3197 lhs = ad.args[0] - bd.args[0]
3198 elif len(ad.args) == len(bd.args):
3199 # should be able to solve
3200 # f(x, y) - f(2 - x, 0) == 0 -> x == 1
3201 raise NotImplementedError(
3202 'equal function with more than 1 argument')
3203 else:
3204 raise ValueError(
3205 'function with different numbers of args')
3207 elif lhs.is_Mul and any(_ispow(a) for a in lhs.args):
3208 lhs = powsimp(powdenest(lhs))
3210 if lhs.is_Function:
3211 if hasattr(lhs, 'inverse') and lhs.inverse() is not None and len(lhs.args) == 1:
3212 # -1
3213 # f(x) = g -> x = f (g)
3214 #
3215 # /!\ inverse should not be defined if there are multiple values
3216 # for the function -- these are handled in _tsolve
3217 #
3218 rhs = lhs.inverse()(rhs)
3219 lhs = lhs.args[0]
3220 elif isinstance(lhs, atan2):
3221 y, x = lhs.args
3222 lhs = 2*atan(y/(sqrt(x**2 + y**2) + x))
3223 elif lhs.func == rhs.func:
3224 if len(lhs.args) == len(rhs.args) == 1:
3225 lhs = lhs.args[0]
3226 rhs = rhs.args[0]
3227 elif len(lhs.args) == len(rhs.args):
3228 # should be able to solve
3229 # f(x, y) == f(2, 3) -> x == 2
3230 # f(x, x + y) == f(2, 3) -> x == 2
3231 raise NotImplementedError(
3232 'equal function with more than 1 argument')
3233 else:
3234 raise ValueError(
3235 'function with different numbers of args')
3238 if rhs and lhs.is_Pow and lhs.exp.is_Integer and lhs.exp < 0:
3239 lhs = 1/lhs
3240 rhs = 1/rhs
3242 # base**a = b -> base = b**(1/a) if
3243 # a is an Integer and dointpow=True (this gives real branch of root)
3244 # a is not an Integer and the equation is multivariate and the
3245 # base has more than 1 symbol in it
3246 # The rationale for this is that right now the multi-system solvers
3247 # doesn't try to resolve generators to see, for example, if the whole
3248 # system is written in terms of sqrt(x + y) so it will just fail, so we
3249 # do that step here.
3250 if lhs.is_Pow and (
3251 lhs.exp.is_Integer and dointpow or not lhs.exp.is_Integer and
3252 len(symbols) > 1 and len(lhs.base.free_symbols & set(symbols)) > 1):
3253 rhs = rhs**(1/lhs.exp)
3254 lhs = lhs.base
3256 if lhs == was:
3257 break
3258 return rhs, lhs
3261def unrad(eq, *syms, **flags):
3262 """
3263 Remove radicals with symbolic arguments and return (eq, cov),
3264 None, or raise an error.
3266 Explanation
3267 ===========
3269 None is returned if there are no radicals to remove.
3271 NotImplementedError is raised if there are radicals and they cannot be
3272 removed or if the relationship between the original symbols and the
3273 change of variable needed to rewrite the system as a polynomial cannot
3274 be solved.
3276 Otherwise the tuple, ``(eq, cov)``, is returned where:
3278 *eq*, ``cov``
3279 *eq* is an equation without radicals (in the symbol(s) of
3280 interest) whose solutions are a superset of the solutions to the
3281 original expression. *eq* might be rewritten in terms of a new
3282 variable; the relationship to the original variables is given by
3283 ``cov`` which is a list containing ``v`` and ``v**p - b`` where
3284 ``p`` is the power needed to clear the radical and ``b`` is the
3285 radical now expressed as a polynomial in the symbols of interest.
3286 For example, for sqrt(2 - x) the tuple would be
3287 ``(c, c**2 - 2 + x)``. The solutions of *eq* will contain
3288 solutions to the original equation (if there are any).
3290 *syms*
3291 An iterable of symbols which, if provided, will limit the focus of
3292 radical removal: only radicals with one or more of the symbols of
3293 interest will be cleared. All free symbols are used if *syms* is not
3294 set.
3296 *flags* are used internally for communication during recursive calls.
3297 Two options are also recognized:
3299 ``take``, when defined, is interpreted as a single-argument function
3300 that returns True if a given Pow should be handled.
3302 Radicals can be removed from an expression if:
3304 * All bases of the radicals are the same; a change of variables is
3305 done in this case.
3306 * If all radicals appear in one term of the expression.
3307 * There are only four terms with sqrt() factors or there are less than
3308 four terms having sqrt() factors.
3309 * There are only two terms with radicals.
3311 Examples
3312 ========
3314 >>> from sympy.solvers.solvers import unrad
3315 >>> from sympy.abc import x
3316 >>> from sympy import sqrt, Rational, root
3318 >>> unrad(sqrt(x)*x**Rational(1, 3) + 2)
3319 (x**5 - 64, [])
3320 >>> unrad(sqrt(x) + root(x + 1, 3))
3321 (-x**3 + x**2 + 2*x + 1, [])
3322 >>> eq = sqrt(x) + root(x, 3) - 2
3323 >>> unrad(eq)
3324 (_p**3 + _p**2 - 2, [_p, _p**6 - x])
3326 """
3328 uflags = {"check": False, "simplify": False}
3330 def _cov(p, e):
3331 if cov:
3332 # XXX - uncovered
3333 oldp, olde = cov
3334 if Poly(e, p).degree(p) in (1, 2):
3335 cov[:] = [p, olde.subs(oldp, _vsolve(e, p, **uflags)[0])]
3336 else:
3337 raise NotImplementedError
3338 else:
3339 cov[:] = [p, e]
3341 def _canonical(eq, cov):
3342 if cov:
3343 # change symbol to vanilla so no solutions are eliminated
3344 p, e = cov
3345 rep = {p: Dummy(p.name)}
3346 eq = eq.xreplace(rep)
3347 cov = [p.xreplace(rep), e.xreplace(rep)]
3349 # remove constants and powers of factors since these don't change
3350 # the location of the root; XXX should factor or factor_terms be used?
3351 eq = factor_terms(_mexpand(eq.as_numer_denom()[0], recursive=True), clear=True)
3352 if eq.is_Mul:
3353 args = []
3354 for f in eq.args:
3355 if f.is_number:
3356 continue
3357 if f.is_Pow:
3358 args.append(f.base)
3359 else:
3360 args.append(f)
3361 eq = Mul(*args) # leave as Mul for more efficient solving
3363 # make the sign canonical
3364 margs = list(Mul.make_args(eq))
3365 changed = False
3366 for i, m in enumerate(margs):
3367 if m.could_extract_minus_sign():
3368 margs[i] = -m
3369 changed = True
3370 if changed:
3371 eq = Mul(*margs, evaluate=False)
3373 return eq, cov
3375 def _Q(pow):
3376 # return leading Rational of denominator of Pow's exponent
3377 c = pow.as_base_exp()[1].as_coeff_Mul()[0]
3378 if not c.is_Rational:
3379 return S.One
3380 return c.q
3382 # define the _take method that will determine whether a term is of interest
3383 def _take(d):
3384 # return True if coefficient of any factor's exponent's den is not 1
3385 for pow in Mul.make_args(d):
3386 if not pow.is_Pow:
3387 continue
3388 if _Q(pow) == 1:
3389 continue
3390 if pow.free_symbols & syms:
3391 return True
3392 return False
3393 _take = flags.setdefault('_take', _take)
3395 if isinstance(eq, Eq):
3396 eq = eq.lhs - eq.rhs # XXX legacy Eq as Eqn support
3397 elif not isinstance(eq, Expr):
3398 return
3400 cov, nwas, rpt = [flags.setdefault(k, v) for k, v in
3401 sorted({"cov": [], "n": None, "rpt": 0}.items())]
3403 # preconditioning
3404 eq = powdenest(factor_terms(eq, radical=True, clear=True))
3405 eq = eq.as_numer_denom()[0]
3406 eq = _mexpand(eq, recursive=True)
3407 if eq.is_number:
3408 return
3410 # see if there are radicals in symbols of interest
3411 syms = set(syms) or eq.free_symbols # _take uses this
3412 poly = eq.as_poly()
3413 gens = [g for g in poly.gens if _take(g)]
3414 if not gens:
3415 return
3417 # recast poly in terms of eigen-gens
3418 poly = eq.as_poly(*gens)
3420 # not a polynomial e.g. 1 + sqrt(x)*exp(sqrt(x)) with gen sqrt(x)
3421 if poly is None:
3422 return
3424 # - an exponent has a symbol of interest (don't handle)
3425 if any(g.exp.has(*syms) for g in gens):
3426 return
3428 def _rads_bases_lcm(poly):
3429 # if all the bases are the same or all the radicals are in one
3430 # term, `lcm` will be the lcm of the denominators of the
3431 # exponents of the radicals
3432 lcm = 1
3433 rads = set()
3434 bases = set()
3435 for g in poly.gens:
3436 q = _Q(g)
3437 if q != 1:
3438 rads.add(g)
3439 lcm = ilcm(lcm, q)
3440 bases.add(g.base)
3441 return rads, bases, lcm
3442 rads, bases, lcm = _rads_bases_lcm(poly)
3444 covsym = Dummy('p', nonnegative=True)
3446 # only keep in syms symbols that actually appear in radicals;
3447 # and update gens
3448 newsyms = set()
3449 for r in rads:
3450 newsyms.update(syms & r.free_symbols)
3451 if newsyms != syms:
3452 syms = newsyms
3453 # get terms together that have common generators
3454 drad = dict(zip(rads, range(len(rads))))
3455 rterms = {(): []}
3456 args = Add.make_args(poly.as_expr())
3457 for t in args:
3458 if _take(t):
3459 common = set(t.as_poly().gens).intersection(rads)
3460 key = tuple(sorted([drad[i] for i in common]))
3461 else:
3462 key = ()
3463 rterms.setdefault(key, []).append(t)
3464 others = Add(*rterms.pop(()))
3465 rterms = [Add(*rterms[k]) for k in rterms.keys()]
3467 # the output will depend on the order terms are processed, so
3468 # make it canonical quickly
3469 rterms = list(reversed(list(ordered(rterms))))
3471 ok = False # we don't have a solution yet
3472 depth = sqrt_depth(eq)
3474 if len(rterms) == 1 and not (rterms[0].is_Add and lcm > 2):
3475 eq = rterms[0]**lcm - ((-others)**lcm)
3476 ok = True
3477 else:
3478 if len(rterms) == 1 and rterms[0].is_Add:
3479 rterms = list(rterms[0].args)
3480 if len(bases) == 1:
3481 b = bases.pop()
3482 if len(syms) > 1:
3483 x = b.free_symbols
3484 else:
3485 x = syms
3486 x = list(ordered(x))[0]
3487 try:
3488 inv = _vsolve(covsym**lcm - b, x, **uflags)
3489 if not inv:
3490 raise NotImplementedError
3491 eq = poly.as_expr().subs(b, covsym**lcm).subs(x, inv[0])
3492 _cov(covsym, covsym**lcm - b)
3493 return _canonical(eq, cov)
3494 except NotImplementedError:
3495 pass
3497 if len(rterms) == 2:
3498 if not others:
3499 eq = rterms[0]**lcm - (-rterms[1])**lcm
3500 ok = True
3501 elif not log(lcm, 2).is_Integer:
3502 # the lcm-is-power-of-two case is handled below
3503 r0, r1 = rterms
3504 if flags.get('_reverse', False):
3505 r1, r0 = r0, r1
3506 i0 = _rads0, _bases0, lcm0 = _rads_bases_lcm(r0.as_poly())
3507 i1 = _rads1, _bases1, lcm1 = _rads_bases_lcm(r1.as_poly())
3508 for reverse in range(2):
3509 if reverse:
3510 i0, i1 = i1, i0
3511 r0, r1 = r1, r0
3512 _rads1, _, lcm1 = i1
3513 _rads1 = Mul(*_rads1)
3514 t1 = _rads1**lcm1
3515 c = covsym**lcm1 - t1
3516 for x in syms:
3517 try:
3518 sol = _vsolve(c, x, **uflags)
3519 if not sol:
3520 raise NotImplementedError
3521 neweq = r0.subs(x, sol[0]) + covsym*r1/_rads1 + \
3522 others
3523 tmp = unrad(neweq, covsym)
3524 if tmp:
3525 eq, newcov = tmp
3526 if newcov:
3527 newp, newc = newcov
3528 _cov(newp, c.subs(covsym,
3529 _vsolve(newc, covsym, **uflags)[0]))
3530 else:
3531 _cov(covsym, c)
3532 else:
3533 eq = neweq
3534 _cov(covsym, c)
3535 ok = True
3536 break
3537 except NotImplementedError:
3538 if reverse:
3539 raise NotImplementedError(
3540 'no successful change of variable found')
3541 else:
3542 pass
3543 if ok:
3544 break
3545 elif len(rterms) == 3:
3546 # two cube roots and another with order less than 5
3547 # (so an analytical solution can be found) or a base
3548 # that matches one of the cube root bases
3549 info = [_rads_bases_lcm(i.as_poly()) for i in rterms]
3550 RAD = 0
3551 BASES = 1
3552 LCM = 2
3553 if info[0][LCM] != 3:
3554 info.append(info.pop(0))
3555 rterms.append(rterms.pop(0))
3556 elif info[1][LCM] != 3:
3557 info.append(info.pop(1))
3558 rterms.append(rterms.pop(1))
3559 if info[0][LCM] == info[1][LCM] == 3:
3560 if info[1][BASES] != info[2][BASES]:
3561 info[0], info[1] = info[1], info[0]
3562 rterms[0], rterms[1] = rterms[1], rterms[0]
3563 if info[1][BASES] == info[2][BASES]:
3564 eq = rterms[0]**3 + (rterms[1] + rterms[2] + others)**3
3565 ok = True
3566 elif info[2][LCM] < 5:
3567 # a*root(A, 3) + b*root(B, 3) + others = c
3568 a, b, c, d, A, B = [Dummy(i) for i in 'abcdAB']
3569 # zz represents the unraded expression into which the
3570 # specifics for this case are substituted
3571 zz = (c - d)*(A**3*a**9 + 3*A**2*B*a**6*b**3 -
3572 3*A**2*a**6*c**3 + 9*A**2*a**6*c**2*d - 9*A**2*a**6*c*d**2 +
3573 3*A**2*a**6*d**3 + 3*A*B**2*a**3*b**6 + 21*A*B*a**3*b**3*c**3 -
3574 63*A*B*a**3*b**3*c**2*d + 63*A*B*a**3*b**3*c*d**2 -
3575 21*A*B*a**3*b**3*d**3 + 3*A*a**3*c**6 - 18*A*a**3*c**5*d +
3576 45*A*a**3*c**4*d**2 - 60*A*a**3*c**3*d**3 + 45*A*a**3*c**2*d**4 -
3577 18*A*a**3*c*d**5 + 3*A*a**3*d**6 + B**3*b**9 - 3*B**2*b**6*c**3 +
3578 9*B**2*b**6*c**2*d - 9*B**2*b**6*c*d**2 + 3*B**2*b**6*d**3 +
3579 3*B*b**3*c**6 - 18*B*b**3*c**5*d + 45*B*b**3*c**4*d**2 -
3580 60*B*b**3*c**3*d**3 + 45*B*b**3*c**2*d**4 - 18*B*b**3*c*d**5 +
3581 3*B*b**3*d**6 - c**9 + 9*c**8*d - 36*c**7*d**2 + 84*c**6*d**3 -
3582 126*c**5*d**4 + 126*c**4*d**5 - 84*c**3*d**6 + 36*c**2*d**7 -
3583 9*c*d**8 + d**9)
3584 def _t(i):
3585 b = Mul(*info[i][RAD])
3586 return cancel(rterms[i]/b), Mul(*info[i][BASES])
3587 aa, AA = _t(0)
3588 bb, BB = _t(1)
3589 cc = -rterms[2]
3590 dd = others
3591 eq = zz.xreplace(dict(zip(
3592 (a, A, b, B, c, d),
3593 (aa, AA, bb, BB, cc, dd))))
3594 ok = True
3595 # handle power-of-2 cases
3596 if not ok:
3597 if log(lcm, 2).is_Integer and (not others and
3598 len(rterms) == 4 or len(rterms) < 4):
3599 def _norm2(a, b):
3600 return a**2 + b**2 + 2*a*b
3602 if len(rterms) == 4:
3603 # (r0+r1)**2 - (r2+r3)**2
3604 r0, r1, r2, r3 = rterms
3605 eq = _norm2(r0, r1) - _norm2(r2, r3)
3606 ok = True
3607 elif len(rterms) == 3:
3608 # (r1+r2)**2 - (r0+others)**2
3609 r0, r1, r2 = rterms
3610 eq = _norm2(r1, r2) - _norm2(r0, others)
3611 ok = True
3612 elif len(rterms) == 2:
3613 # r0**2 - (r1+others)**2
3614 r0, r1 = rterms
3615 eq = r0**2 - _norm2(r1, others)
3616 ok = True
3618 new_depth = sqrt_depth(eq) if ok else depth
3619 rpt += 1 # XXX how many repeats with others unchanging is enough?
3620 if not ok or (
3621 nwas is not None and len(rterms) == nwas and
3622 new_depth is not None and new_depth == depth and
3623 rpt > 3):
3624 raise NotImplementedError('Cannot remove all radicals')
3626 flags.update({"cov": cov, "n": len(rterms), "rpt": rpt})
3627 neq = unrad(eq, *syms, **flags)
3628 if neq:
3629 eq, cov = neq
3630 eq, cov = _canonical(eq, cov)
3631 return eq, cov
3634# delayed imports
3635from sympy.solvers.bivariate import (
3636 bivariate_type, _solve_lambert, _filtered_gens)