-
-
Notifications
You must be signed in to change notification settings - Fork 4.3k
/
solvers.py
3640 lines (3140 loc) · 133 KB
/
solvers.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
"""
This module contain solvers for all kinds of equations:
- algebraic or transcendental, use solve()
- recurrence, use rsolve()
- differential, use dsolve()
- nonlinear (numerically), use nsolve()
(you will need a good starting point)
"""
from __future__ import annotations
from sympy.core import (S, Add, Symbol, Dummy, Expr, Mul)
from sympy.core.assumptions import check_assumptions
from sympy.core.exprtools import factor_terms
from sympy.core.function import (expand_mul, expand_log, Derivative,
AppliedUndef, UndefinedFunction, nfloat,
Function, expand_power_exp, _mexpand, expand,
expand_func)
from sympy.core.logic import fuzzy_not
from sympy.core.numbers import ilcm, Float, Rational, _illegal
from sympy.core.power import integer_log, Pow
from sympy.core.relational import Eq, Ne
from sympy.core.sorting import ordered, default_sort_key
from sympy.core.sympify import sympify, _sympify
from sympy.core.traversal import preorder_traversal
from sympy.logic.boolalg import And, BooleanAtom
from sympy.functions import (log, exp, LambertW, cos, sin, tan, acos, asin, atan,
Abs, re, im, arg, sqrt, atan2)
from sympy.functions.combinatorial.factorials import binomial
from sympy.functions.elementary.hyperbolic import HyperbolicFunction
from sympy.functions.elementary.piecewise import piecewise_fold, Piecewise
from sympy.functions.elementary.trigonometric import TrigonometricFunction
from sympy.integrals.integrals import Integral
from sympy.ntheory.factor_ import divisors
from sympy.simplify import (simplify, collect, powsimp, posify, # type: ignore
powdenest, nsimplify, denom, logcombine, sqrtdenest, fraction,
separatevars)
from sympy.simplify.sqrtdenest import sqrt_depth
from sympy.simplify.fu import TR1, TR2i
from sympy.strategies.rl import rebuild
from sympy.matrices.common import NonInvertibleMatrixError
from sympy.matrices import Matrix, zeros
from sympy.polys import roots, cancel, factor, Poly
from sympy.polys.polyerrors import GeneratorsNeeded, PolynomialError
from sympy.polys.solvers import sympy_eqs_to_ring, solve_lin_sys
from sympy.utilities.lambdify import lambdify
from sympy.utilities.misc import filldedent, debugf
from sympy.utilities.iterables import (connected_components,
generate_bell, uniq, iterable, is_sequence, subsets, flatten)
from sympy.utilities.decorator import conserve_mpmath_dps
from mpmath import findroot
from sympy.solvers.polysys import solve_poly_system
from types import GeneratorType
from collections import defaultdict
from itertools import combinations, product
import warnings
def recast_to_symbols(eqs, symbols):
"""
Return (e, s, d) where e and s are versions of *eqs* and
*symbols* in which any non-Symbol objects in *symbols* have
been replaced with generic Dummy symbols and d is a dictionary
that can be used to restore the original expressions.
Examples
========
>>> from sympy.solvers.solvers import recast_to_symbols
>>> from sympy import symbols, Function
>>> x, y = symbols('x y')
>>> fx = Function('f')(x)
>>> eqs, syms = [fx + 1, x, y], [fx, y]
>>> e, s, d = recast_to_symbols(eqs, syms); (e, s, d)
([_X0 + 1, x, y], [_X0, y], {_X0: f(x)})
The original equations and symbols can be restored using d:
>>> assert [i.xreplace(d) for i in eqs] == eqs
>>> assert [d.get(i, i) for i in s] == syms
"""
if not iterable(eqs) and iterable(symbols):
raise ValueError('Both eqs and symbols must be iterable')
orig = list(symbols)
symbols = list(ordered(symbols))
swap_sym = {}
i = 0
for j, s in enumerate(symbols):
if not isinstance(s, Symbol) and s not in swap_sym:
swap_sym[s] = Dummy('X%d' % i)
i += 1
new_f = []
for i in eqs:
isubs = getattr(i, 'subs', None)
if isubs is not None:
new_f.append(isubs(swap_sym))
else:
new_f.append(i)
restore = {v: k for k, v in swap_sym.items()}
return new_f, [swap_sym.get(i, i) for i in orig], restore
def _ispow(e):
"""Return True if e is a Pow or is exp."""
return isinstance(e, Expr) and (e.is_Pow or isinstance(e, exp))
def _simple_dens(f, symbols):
# when checking if a denominator is zero, we can just check the
# base of powers with nonzero exponents since if the base is zero
# the power will be zero, too. To keep it simple and fast, we
# limit simplification to exponents that are Numbers
dens = set()
for d in denoms(f, symbols):
if d.is_Pow and d.exp.is_Number:
if d.exp.is_zero:
continue # foo**0 is never 0
d = d.base
dens.add(d)
return dens
def denoms(eq, *symbols):
"""
Return (recursively) set of all denominators that appear in *eq*
that contain any symbol in *symbols*; if *symbols* are not
provided then all denominators will be returned.
Examples
========
>>> from sympy.solvers.solvers import denoms
>>> from sympy.abc import x, y, z
>>> denoms(x/y)
{y}
>>> denoms(x/(y*z))
{y, z}
>>> denoms(3/x + y/z)
{x, z}
>>> denoms(x/2 + y/z)
{2, z}
If *symbols* are provided then only denominators containing
those symbols will be returned:
>>> denoms(1/x + 1/y + 1/z, y, z)
{y, z}
"""
pot = preorder_traversal(eq)
dens = set()
for p in pot:
# Here p might be Tuple or Relational
# Expr subtrees (e.g. lhs and rhs) will be traversed after by pot
if not isinstance(p, Expr):
continue
den = denom(p)
if den is S.One:
continue
for d in Mul.make_args(den):
dens.add(d)
if not symbols:
return dens
elif len(symbols) == 1:
if iterable(symbols[0]):
symbols = symbols[0]
return {d for d in dens if any(s in d.free_symbols for s in symbols)}
def checksol(f, symbol, sol=None, **flags):
"""
Checks whether sol is a solution of equation f == 0.
Explanation
===========
Input can be either a single symbol and corresponding value
or a dictionary of symbols and values. When given as a dictionary
and flag ``simplify=True``, the values in the dictionary will be
simplified. *f* can be a single equation or an iterable of equations.
A solution must satisfy all equations in *f* to be considered valid;
if a solution does not satisfy any equation, False is returned; if one or
more checks are inconclusive (and none are False) then None is returned.
Examples
========
>>> from sympy import checksol, symbols
>>> x, y = symbols('x,y')
>>> checksol(x**4 - 1, x, 1)
True
>>> checksol(x**4 - 1, x, 0)
False
>>> checksol(x**2 + y**2 - 5**2, {x: 3, y: 4})
True
To check if an expression is zero using ``checksol()``, pass it
as *f* and send an empty dictionary for *symbol*:
>>> checksol(x**2 + x - x*(x + 1), {})
True
None is returned if ``checksol()`` could not conclude.
flags:
'numerical=True (default)'
do a fast numerical check if ``f`` has only one symbol.
'minimal=True (default is False)'
a very fast, minimal testing.
'warn=True (default is False)'
show a warning if checksol() could not conclude.
'simplify=True (default)'
simplify solution before substituting into function and
simplify the function before trying specific simplifications
'force=True (default is False)'
make positive all symbols without assumptions regarding sign.
"""
from sympy.physics.units import Unit
minimal = flags.get('minimal', False)
if sol is not None:
sol = {symbol: sol}
elif isinstance(symbol, dict):
sol = symbol
else:
msg = 'Expecting (sym, val) or ({sym: val}, None) but got (%s, %s)'
raise ValueError(msg % (symbol, sol))
if iterable(f):
if not f:
raise ValueError('no functions to check')
rv = True
for fi in f:
check = checksol(fi, sol, **flags)
if check:
continue
if check is False:
return False
rv = None # don't return, wait to see if there's a False
return rv
f = _sympify(f)
if f.is_number:
return f.is_zero
if isinstance(f, Poly):
f = f.as_expr()
elif isinstance(f, (Eq, Ne)):
if f.rhs in (S.true, S.false):
f = f.reversed
B, E = f.args
if isinstance(B, BooleanAtom):
f = f.subs(sol)
if not f.is_Boolean:
return
else:
f = f.rewrite(Add, evaluate=False, deep=False)
if isinstance(f, BooleanAtom):
return bool(f)
elif not f.is_Relational and not f:
return True
illegal = set(_illegal)
if any(sympify(v).atoms() & illegal for k, v in sol.items()):
return False
attempt = -1
numerical = flags.get('numerical', True)
while 1:
attempt += 1
if attempt == 0:
val = f.subs(sol)
if isinstance(val, Mul):
val = val.as_independent(Unit)[0]
if val.atoms() & illegal:
return False
elif attempt == 1:
if not val.is_number:
if not val.is_constant(*list(sol.keys()), simplify=not minimal):
return False
# there are free symbols -- simple expansion might work
_, val = val.as_content_primitive()
val = _mexpand(val.as_numer_denom()[0], recursive=True)
elif attempt == 2:
if minimal:
return
if flags.get('simplify', True):
for k in sol:
sol[k] = simplify(sol[k])
# start over without the failed expanded form, possibly
# with a simplified solution
val = simplify(f.subs(sol))
if flags.get('force', True):
val, reps = posify(val)
# expansion may work now, so try again and check
exval = _mexpand(val, recursive=True)
if exval.is_number:
# we can decide now
val = exval
else:
# if there are no radicals and no functions then this can't be
# zero anymore -- can it?
pot = preorder_traversal(expand_mul(val))
seen = set()
saw_pow_func = False
for p in pot:
if p in seen:
continue
seen.add(p)
if p.is_Pow and not p.exp.is_Integer:
saw_pow_func = True
elif p.is_Function:
saw_pow_func = True
elif isinstance(p, UndefinedFunction):
saw_pow_func = True
if saw_pow_func:
break
if saw_pow_func is False:
return False
if flags.get('force', True):
# don't do a zero check with the positive assumptions in place
val = val.subs(reps)
nz = fuzzy_not(val.is_zero)
if nz is not None:
# issue 5673: nz may be True even when False
# so these are just hacks to keep a false positive
# from being returned
# HACK 1: LambertW (issue 5673)
if val.is_number and val.has(LambertW):
# don't eval this to verify solution since if we got here,
# numerical must be False
return None
# add other HACKs here if necessary, otherwise we assume
# the nz value is correct
return not nz
break
if val.is_Rational:
return val == 0
if numerical and val.is_number:
return (abs(val.n(18).n(12, chop=True)) < 1e-9) is S.true
if flags.get('warn', False):
warnings.warn("\n\tWarning: could not verify solution %s." % sol)
# returns None if it can't conclude
# TODO: improve solution testing
def solve(f, *symbols, **flags):
r"""
Algebraically solves equations and systems of equations.
Explanation
===========
Currently supported:
- polynomial
- transcendental
- piecewise combinations of the above
- systems of linear and polynomial equations
- systems containing relational expressions
- systems implied by undetermined coefficients
Examples
========
The default output varies according to the input and might
be a list (possibly empty), a dictionary, a list of
dictionaries or tuples, or an expression involving relationals.
For specifics regarding different forms of output that may appear, see :ref:`solve_output`.
Let it suffice here to say that to obtain a uniform output from
`solve` use ``dict=True`` or ``set=True`` (see below).
>>> from sympy import solve, Poly, Eq, Matrix, Symbol
>>> from sympy.abc import x, y, z, a, b
The expressions that are passed can be Expr, Equality, or Poly
classes (or lists of the same); a Matrix is considered to be a
list of all the elements of the matrix:
>>> solve(x - 3, x)
[3]
>>> solve(Eq(x, 3), x)
[3]
>>> solve(Poly(x - 3), x)
[3]
>>> solve(Matrix([[x, x + y]]), x, y) == solve([x, x + y], x, y)
True
If no symbols are indicated to be of interest and the equation is
univariate, a list of values is returned; otherwise, the keys in
a dictionary will indicate which (of all the variables used in
the expression(s)) variables and solutions were found:
>>> solve(x**2 - 4)
[-2, 2]
>>> solve((x - a)*(y - b))
[{a: x}, {b: y}]
>>> solve([x - 3, y - 1])
{x: 3, y: 1}
>>> solve([x - 3, y**2 - 1])
[{x: 3, y: -1}, {x: 3, y: 1}]
If you pass symbols for which solutions are sought, the output will vary
depending on the number of symbols you passed, whether you are passing
a list of expressions or not, and whether a linear system was solved.
Uniform output is attained by using ``dict=True`` or ``set=True``.
>>> #### *** feel free to skip to the stars below *** ####
>>> from sympy import TableForm
>>> h = [None, ';|;'.join(['e', 's', 'solve(e, s)', 'solve(e, s, dict=True)',
... 'solve(e, s, set=True)']).split(';')]
>>> t = []
>>> for e, s in [
... (x - y, y),
... (x - y, [x, y]),
... (x**2 - y, [x, y]),
... ([x - 3, y -1], [x, y]),
... ]:
... how = [{}, dict(dict=True), dict(set=True)]
... res = [solve(e, s, **f) for f in how]
... t.append([e, '|', s, '|'] + [res[0], '|', res[1], '|', res[2]])
...
>>> # ******************************************************* #
>>> TableForm(t, headings=h, alignments="<")
e | s | solve(e, s) | solve(e, s, dict=True) | solve(e, s, set=True)
---------------------------------------------------------------------------------------
x - y | y | [x] | [{y: x}] | ([y], {(x,)})
x - y | [x, y] | [(y, y)] | [{x: y}] | ([x, y], {(y, y)})
x**2 - y | [x, y] | [(x, x**2)] | [{y: x**2}] | ([x, y], {(x, x**2)})
[x - 3, y - 1] | [x, y] | {x: 3, y: 1} | [{x: 3, y: 1}] | ([x, y], {(3, 1)})
* If any equation does not depend on the symbol(s) given, it will be
eliminated from the equation set and an answer may be given
implicitly in terms of variables that were not of interest:
>>> solve([x - y, y - 3], x)
{x: y}
When you pass all but one of the free symbols, an attempt
is made to find a single solution based on the method of
undetermined coefficients. If it succeeds, a dictionary of values
is returned. If you want an algebraic solutions for one
or more of the symbols, pass the expression to be solved in a list:
>>> e = a*x + b - 2*x - 3
>>> solve(e, [a, b])
{a: 2, b: 3}
>>> solve([e], [a, b])
{a: -b/x + (2*x + 3)/x}
When there is no solution for any given symbol which will make all
expressions zero, the empty list is returned (or an empty set in
the tuple when ``set=True``):
>>> from sympy import sqrt
>>> solve(3, x)
[]
>>> solve(x - 3, y)
[]
>>> solve(sqrt(x) + 1, x, set=True)
([x], set())
When an object other than a Symbol is given as a symbol, it is
isolated algebraically and an implicit solution may be obtained.
This is mostly provided as a convenience to save you from replacing
the object with a Symbol and solving for that Symbol. It will only
work if the specified object can be replaced with a Symbol using the
subs method:
>>> from sympy import exp, Function
>>> f = Function('f')
>>> solve(f(x) - x, f(x))
[x]
>>> solve(f(x).diff(x) - f(x) - x, f(x).diff(x))
[x + f(x)]
>>> solve(f(x).diff(x) - f(x) - x, f(x))
[-x + Derivative(f(x), x)]
>>> solve(x + exp(x)**2, exp(x), set=True)
([exp(x)], {(-sqrt(-x),), (sqrt(-x),)})
>>> from sympy import Indexed, IndexedBase, Tuple
>>> A = IndexedBase('A')
>>> eqs = Tuple(A[1] + A[2] - 3, A[1] - A[2] + 1)
>>> solve(eqs, eqs.atoms(Indexed))
{A[1]: 1, A[2]: 2}
* To solve for a function within a derivative, use :func:`~.dsolve`.
To solve for a symbol implicitly, use implicit=True:
>>> solve(x + exp(x), x)
[-LambertW(1)]
>>> solve(x + exp(x), x, implicit=True)
[-exp(x)]
It is possible to solve for anything in an expression that can be
replaced with a symbol using :obj:`~sympy.core.basic.Basic.subs`:
>>> solve(x + 2 + sqrt(3), x + 2)
[-sqrt(3)]
>>> solve((x + 2 + sqrt(3), x + 4 + y), y, x + 2)
{y: -2 + sqrt(3), x + 2: -sqrt(3)}
* Nothing heroic is done in this implicit solving so you may end up
with a symbol still in the solution:
>>> eqs = (x*y + 3*y + sqrt(3), x + 4 + y)
>>> solve(eqs, y, x + 2)
{y: -sqrt(3)/(x + 3), x + 2: -2*x/(x + 3) - 6/(x + 3) + sqrt(3)/(x + 3)}
>>> solve(eqs, y*x, x)
{x: -y - 4, x*y: -3*y - sqrt(3)}
* If you attempt to solve for a number, remember that the number
you have obtained does not necessarily mean that the value is
equivalent to the expression obtained:
>>> solve(sqrt(2) - 1, 1)
[sqrt(2)]
>>> solve(x - y + 1, 1) # /!\ -1 is targeted, too
[x/(y - 1)]
>>> [_.subs(z, -1) for _ in solve((x - y + 1).subs(-1, z), 1)]
[-x + y]
**Additional Examples**
``solve()`` with check=True (default) will run through the symbol tags to
eliminate unwanted solutions. If no assumptions are included, all possible
solutions will be returned:
>>> x = Symbol("x")
>>> solve(x**2 - 1)
[-1, 1]
By setting the ``positive`` flag, only one solution will be returned:
>>> pos = Symbol("pos", positive=True)
>>> solve(pos**2 - 1)
[1]
When the solutions are checked, those that make any denominator zero
are automatically excluded. If you do not want to exclude such solutions,
then use the check=False option:
>>> from sympy import sin, limit
>>> solve(sin(x)/x) # 0 is excluded
[pi]
If ``check=False``, then a solution to the numerator being zero is found
but the value of $x = 0$ is a spurious solution since $\sin(x)/x$ has the well
known limit (without discontinuity) of 1 at $x = 0$:
>>> solve(sin(x)/x, check=False)
[0, pi]
In the following case, however, the limit exists and is equal to the
value of $x = 0$ that is excluded when check=True:
>>> eq = x**2*(1/x - z**2/x)
>>> solve(eq, x)
[]
>>> solve(eq, x, check=False)
[0]
>>> limit(eq, x, 0, '-')
0
>>> limit(eq, x, 0, '+')
0
**Solving Relationships**
When one or more expressions passed to ``solve`` is a relational,
a relational result is returned (and the ``dict`` and ``set`` flags
are ignored):
>>> solve(x < 3)
(-oo < x) & (x < 3)
>>> solve([x < 3, x**2 > 4], x)
((-oo < x) & (x < -2)) | ((2 < x) & (x < 3))
>>> solve([x + y - 3, x > 3], x)
(3 < x) & (x < oo) & Eq(x, 3 - y)
Although checking of assumptions on symbols in relationals
is not done, setting assumptions will affect how certain
relationals might automatically simplify:
>>> solve(x**2 > 4)
((-oo < x) & (x < -2)) | ((2 < x) & (x < oo))
>>> r = Symbol('r', real=True)
>>> solve(r**2 > 4)
(2 < r) | (r < -2)
There is currently no algorithm in SymPy that allows you to use
relationships to resolve more than one variable. So the following
does not determine that ``q < 0`` (and trying to solve for ``r``
and ``q`` will raise an error):
>>> from sympy import symbols
>>> r, q = symbols('r, q', real=True)
>>> solve([r + q - 3, r > 3], r)
(3 < r) & Eq(r, 3 - q)
You can directly call the routine that ``solve`` calls
when it encounters a relational: :func:`~.reduce_inequalities`.
It treats Expr like Equality.
>>> from sympy import reduce_inequalities
>>> reduce_inequalities([x**2 - 4])
Eq(x, -2) | Eq(x, 2)
If each relationship contains only one symbol of interest,
the expressions can be processed for multiple symbols:
>>> reduce_inequalities([0 <= x - 1, y < 3], [x, y])
(-oo < y) & (1 <= x) & (x < oo) & (y < 3)
But an error is raised if any relationship has more than one
symbol of interest:
>>> reduce_inequalities([0 <= x*y - 1, y < 3], [x, y])
Traceback (most recent call last):
...
NotImplementedError:
inequality has more than one symbol of interest.
**Disabling High-Order Explicit Solutions**
When solving polynomial expressions, you might not want explicit solutions
(which can be quite long). If the expression is univariate, ``CRootOf``
instances will be returned instead:
>>> solve(x**3 - x + 1)
[-1/((-1/2 - sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)) -
(-1/2 - sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3,
-(-1/2 + sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)/3 -
1/((-1/2 + sqrt(3)*I/2)*(3*sqrt(69)/2 + 27/2)**(1/3)),
-(3*sqrt(69)/2 + 27/2)**(1/3)/3 -
1/(3*sqrt(69)/2 + 27/2)**(1/3)]
>>> solve(x**3 - x + 1, cubics=False)
[CRootOf(x**3 - x + 1, 0),
CRootOf(x**3 - x + 1, 1),
CRootOf(x**3 - x + 1, 2)]
If the expression is multivariate, no solution might be returned:
>>> solve(x**3 - x + a, x, cubics=False)
[]
Sometimes solutions will be obtained even when a flag is False because the
expression could be factored. In the following example, the equation can
be factored as the product of a linear and a quadratic factor so explicit
solutions (which did not require solving a cubic expression) are obtained:
>>> eq = x**3 + 3*x**2 + x - 1
>>> solve(eq, cubics=False)
[-1, -1 + sqrt(2), -sqrt(2) - 1]
**Solving Equations Involving Radicals**
Because of SymPy's use of the principle root, some solutions
to radical equations will be missed unless check=False:
>>> from sympy import root
>>> eq = root(x**3 - 3*x**2, 3) + 1 - x
>>> solve(eq)
[]
>>> solve(eq, check=False)
[1/3]
In the above example, there is only a single solution to the
equation. Other expressions will yield spurious roots which
must be checked manually; roots which give a negative argument
to odd-powered radicals will also need special checking:
>>> from sympy import real_root, S
>>> eq = root(x, 3) - root(x, 5) + S(1)/7
>>> solve(eq) # this gives 2 solutions but misses a 3rd
[CRootOf(7*x**5 - 7*x**3 + 1, 1)**15,
CRootOf(7*x**5 - 7*x**3 + 1, 2)**15]
>>> sol = solve(eq, check=False)
>>> [abs(eq.subs(x,i).n(2)) for i in sol]
[0.48, 0.e-110, 0.e-110, 0.052, 0.052]
The first solution is negative so ``real_root`` must be used to see that it
satisfies the expression:
>>> abs(real_root(eq.subs(x, sol[0])).n(2))
0.e-110
If the roots of the equation are not real then more care will be
necessary to find the roots, especially for higher order equations.
Consider the following expression:
>>> expr = root(x, 3) - root(x, 5)
We will construct a known value for this expression at x = 3 by selecting
the 1-th root for each radical:
>>> expr1 = root(x, 3, 1) - root(x, 5, 1)
>>> v = expr1.subs(x, -3)
The ``solve`` function is unable to find any exact roots to this equation:
>>> eq = Eq(expr, v); eq1 = Eq(expr1, v)
>>> solve(eq, check=False), solve(eq1, check=False)
([], [])
The function ``unrad``, however, can be used to get a form of the equation
for which numerical roots can be found:
>>> from sympy.solvers.solvers import unrad
>>> from sympy import nroots
>>> e, (p, cov) = unrad(eq)
>>> pvals = nroots(e)
>>> inversion = solve(cov, x)[0]
>>> xvals = [inversion.subs(p, i) for i in pvals]
Although ``eq`` or ``eq1`` could have been used to find ``xvals``, the
solution can only be verified with ``expr1``:
>>> z = expr - v
>>> [xi.n(chop=1e-9) for xi in xvals if abs(z.subs(x, xi).n()) < 1e-9]
[]
>>> z1 = expr1 - v
>>> [xi.n(chop=1e-9) for xi in xvals if abs(z1.subs(x, xi).n()) < 1e-9]
[-3.0]
Parameters
==========
f :
- a single Expr or Poly that must be zero
- an Equality
- a Relational expression
- a Boolean
- iterable of one or more of the above
symbols : (object(s) to solve for) specified as
- none given (other non-numeric objects will be used)
- single symbol
- denested list of symbols
(e.g., ``solve(f, x, y)``)
- ordered iterable of symbols
(e.g., ``solve(f, [x, y])``)
flags :
dict=True (default is False)
Return list (perhaps empty) of solution mappings.
set=True (default is False)
Return list of symbols and set of tuple(s) of solution(s).
exclude=[] (default)
Do not try to solve for any of the free symbols in exclude;
if expressions are given, the free symbols in them will
be extracted automatically.
check=True (default)
If False, do not do any testing of solutions. This can be
useful if you want to include solutions that make any
denominator zero.
numerical=True (default)
Do a fast numerical check if *f* has only one symbol.
minimal=True (default is False)
A very fast, minimal testing.
warn=True (default is False)
Show a warning if ``checksol()`` could not conclude.
simplify=True (default)
Simplify all but polynomials of order 3 or greater before
returning them and (if check is not False) use the
general simplify function on the solutions and the
expression obtained when they are substituted into the
function which should be zero.
force=True (default is False)
Make positive all symbols without assumptions regarding sign.
rational=True (default)
Recast Floats as Rational; if this option is not used, the
system containing Floats may fail to solve because of issues
with polys. If rational=None, Floats will be recast as
rationals but the answer will be recast as Floats. If the
flag is False then nothing will be done to the Floats.
manual=True (default is False)
Do not use the polys/matrix method to solve a system of
equations, solve them one at a time as you might "manually."
implicit=True (default is False)
Allows ``solve`` to return a solution for a pattern in terms of
other functions that contain that pattern; this is only
needed if the pattern is inside of some invertible function
like cos, exp, ect.
particular=True (default is False)
Instructs ``solve`` to try to find a particular solution to
a linear system with as many zeros as possible; this is very
expensive.
quick=True (default is False; ``particular`` must be True)
Selects a fast heuristic to find a solution with many zeros
whereas a value of False uses the very slow method guaranteed
to find the largest number of zeros possible.
cubics=True (default)
Return explicit solutions when cubic expressions are encountered.
When False, quartics and quintics are disabled, too.
quartics=True (default)
Return explicit solutions when quartic expressions are encountered.
When False, quintics are disabled, too.
quintics=True (default)
Return explicit solutions (if possible) when quintic expressions
are encountered.
See Also
========
rsolve: For solving recurrence relationships
dsolve: For solving differential equations
"""
from .inequalities import reduce_inequalities
# checking/recording flags
###########################################################################
# set solver types explicitly; as soon as one is False
# all the rest will be False
hints = ('cubics', 'quartics', 'quintics')
default = True
for k in hints:
default = flags.setdefault(k, bool(flags.get(k, default)))
# allow solution to contain symbol if True:
implicit = flags.get('implicit', False)
# record desire to see warnings
warn = flags.get('warn', False)
# this flag will be needed for quick exits below, so record
# now -- but don't record `dict` yet since it might change
as_set = flags.get('set', False)
# keeping track of how f was passed
bare_f = not iterable(f)
# check flag usage for particular/quick which should only be used
# with systems of equations
if flags.get('quick', None) is not None:
if not flags.get('particular', None):
raise ValueError('when using `quick`, `particular` should be True')
if flags.get('particular', False) and bare_f:
raise ValueError(filldedent("""
The 'particular/quick' flag is usually used with systems of
equations. Either pass your equation in a list or
consider using a solver like `diophantine` if you are
looking for a solution in integers."""))
# sympify everything, creating list of expressions and list of symbols
###########################################################################
def _sympified_list(w):
return list(map(sympify, w if iterable(w) else [w]))
f, symbols = (_sympified_list(w) for w in [f, symbols])
# preprocess symbol(s)
###########################################################################
ordered_symbols = None # were the symbols in a well defined order?
if not symbols:
# get symbols from equations
symbols = set().union(*[fi.free_symbols for fi in f])
if len(symbols) < len(f):
for fi in f:
pot = preorder_traversal(fi)
for p in pot:
if isinstance(p, AppliedUndef):
if not as_set:
flags['dict'] = True # better show symbols
symbols.add(p)
pot.skip() # don't go any deeper
ordered_symbols = False
symbols = list(ordered(symbols)) # to make it canonical
else:
if len(symbols) == 1 and iterable(symbols[0]):
symbols = symbols[0]
ordered_symbols = symbols and is_sequence(symbols,
include=GeneratorType)
_symbols = list(uniq(symbols))
if len(_symbols) != len(symbols):
ordered_symbols = False
symbols = list(ordered(symbols))
else:
symbols = _symbols
# check for duplicates
if len(symbols) != len(set(symbols)):
raise ValueError('duplicate symbols given')
# remove those not of interest
exclude = flags.pop('exclude', set())
if exclude:
if isinstance(exclude, Expr):
exclude = [exclude]
exclude = set().union(*[e.free_symbols for e in sympify(exclude)])
symbols = [s for s in symbols if s not in exclude]
# preprocess equation(s)
###########################################################################
# automatically ignore True values
if isinstance(f, list):
f = [s for s in f if s is not S.true]
# handle canonicalization of equation types
for i, fi in enumerate(f):
if isinstance(fi, (Eq, Ne)):
if 'ImmutableDenseMatrix' in [type(a).__name__ for a in fi.args]:
fi = fi.lhs - fi.rhs
else:
L, R = fi.args
if isinstance(R, BooleanAtom):
L, R = R, L
if isinstance(L, BooleanAtom):
if isinstance(fi, Ne):
L = ~L
if R.is_Relational:
fi = ~R if L is S.false else R
elif R.is_Symbol:
return L
elif R.is_Boolean and (~R).is_Symbol:
return ~L
else:
raise NotImplementedError(filldedent('''
Unanticipated argument of Eq when other arg
is True or False.
'''))
else:
fi = fi.rewrite(Add, evaluate=False, deep=False)
f[i] = fi
# *** dispatch and handle as a system of relationals
# **************************************************
if fi.is_Relational:
if len(symbols) != 1:
raise ValueError("can only solve for one symbol at a time")
if warn and symbols[0].assumptions0:
warnings.warn(filldedent("""
\tWarning: assumptions about variable '%s' are
not handled currently.""" % symbols[0]))
return reduce_inequalities(f, symbols=symbols)
# convert Poly to expression
if isinstance(fi, Poly):
f[i] = fi.as_expr()
# rewrite hyperbolics in terms of exp if they have symbols of
# interest
f[i] = f[i].replace(lambda w: isinstance(w, HyperbolicFunction) and \
w.has_free(*symbols), lambda w: w.rewrite(exp))
# if we have a Matrix, we need to iterate over its elements again
if f[i].is_Matrix:
bare_f = False
f.extend(list(f[i]))
f[i] = S.Zero
# if we can split it into real and imaginary parts then do so
freei = f[i].free_symbols
if freei and all(s.is_extended_real or s.is_imaginary for s in freei):
fr, fi = f[i].as_real_imag()
# accept as long as new re, im, arg or atan2 are not introduced
had = f[i].atoms(re, im, arg, atan2)
if fr and fi and fr != fi and not any(
i.atoms(re, im, arg, atan2) - had for i in (fr, fi)):
if bare_f:
bare_f = False
f[i: i + 1] = [fr, fi]
# real/imag handling -----------------------------
if any(isinstance(fi, (bool, BooleanAtom)) for fi in f):
if as_set:
return [], set()
return []
for i, fi in enumerate(f):
# Abs
while True:
was = fi
fi = fi.replace(Abs, lambda arg: