On Mon, Nov 24, 2008 at 1:35 AM, <[EMAIL PROTECTED]> wrote:
>
> From: Fabian Seoane <[EMAIL PROTECTED]>
>
> For this, a guess_solve_strategy method was added
>
> Before of this, solve() parsed the expression searching for a polynomial
> and if this parse failed, it just called tsolve.
>
> Now the expression is parsed, and while parsing the best algorithm is
> selected for that particular expression, i.e. not only it is conceptually
> clearer and easier to extend, but there is no need to parse the entire
> expression to decide what strategy to use.
>
> New features
> ============
> - Support for more expressions, mostly expressions that can
> be transformed to polynomial expressions using simple changes of
> variables. These incluide:
> - x**(p1/q1) + x**Rational(p2/q2) + ...
> - x**-n1 + x**-n2 + ...
> - P(x)/Q(x) (rational functions)
> - tsolve now returns a tuple of expressions instead only one expression.
> - It is easy to add new types of functions. See docstring of
> solve and guess_solve_strategy for more info
> - Test for all new types of supported expressions
>
> Performace
> ==========
> - If it is a polynomial, it will partially parse the expression twice,
> since it will parse the polynomial inside guess_solve_strategy to
> determine the solve strategy and then it will be parsed again in
> solve.
> - In other cases, it can be faster or it can be slower. For trascebdental
> functions it will generally be faster, since there is no need to
> parse the expression completely
>
> Random notes
> ============
> - The change of variable could be computed during the parse of
> guess_solve_strategy, and so we could avoid to partially parse
> the expression twice. I didn't implement this because it would
> make the code more complex, it is nice right now that functions
> guess_solve_strategy and solve have a clear separate entity.
> Plus, parsing the expression twice is not a bit issue since
> the parsing is fast, O(n) typically, where n is the number of arguments
> for the expression, and most of the time is spended in simplify
> and in subs
> ---
> sympy/solvers/solvers.py | 192
> +++++++++++++++++++++++++++++++----
> sympy/solvers/tests/test_solvers.py | 135 ++++++++++++++++++++-----
> 2 files changed, 279 insertions(+), 48 deletions(-)
>
> diff --git a/sympy/solvers/solvers.py b/sympy/solvers/solvers.py
> index 06d9a8f..6a9daf4 100644
> --- a/sympy/solvers/solvers.py
> +++ b/sympy/solvers/solvers.py
> @@ -14,10 +14,12 @@
> """
>
> from sympy.core.sympify import sympify
> -from sympy.core.basic import Basic, S, C
> +from sympy.core.basic import Basic, S, C, Mul, Add
> +from sympy.core.power import Pow
> from sympy.core.symbol import Symbol, Wild
> from sympy.core.relational import Equality
> from sympy.core.function import Derivative, diff
> +from sympy.core.numbers import ilcm
>
> from sympy.functions import sqrt, log, exp, LambertW
> from sympy.simplify import simplify, collect
> @@ -30,6 +32,86 @@ from sympy.solvers.numeric import newton
>
> from sympy.solvers.polysys import solve_poly_system
>
> +# Codes for guess solve strategy
> +GS_POLY = 0
> +GS_RATIONAL = 1
> +GS_POLY_CV_1 = 2 # can be converted to a polynomial equation via the change
> of variable y -> x**n
> +GS_POLY_CV_2 = 3 # can be converted to a polynomial equation multiplying on
> both sides by x**m
> + # for example, x + 1/x == 0. Multiplying by x yields x**2 +
> x == 0
> +GS_RATIONAL_CV_1 = 4 # can be converted to a rational equation via the
> change of variable y -> x**n
> +GS_TRASCENDENTAL = 5
> +
> +def guess_solve_strategy(expr, symbol):
> + """
> + Tries to guess what approach should be used to solve a specific equation
> +
> + Returns
> + =======
> + - -1: could not guess
> + - integer > 0: code representing certain type of equation. See GS_*
> fields
> + on this module for a complete list
> +
> + Examples
> + ========
> + >>> from sympy import Symbol, Rational
> + >>> x = Symbol('x')
> + >>> guess_solve_strategy(x**2 + 1, x)
> + 0
> + >>> guess_solve_strategy(x**Rational(1,2) + 1, x)
> + 2
> + """
> + eq_type = -1
> + if expr.is_Add:
> + items = expr.args
> + for item in items:
> + if item.is_Number or item.is_Symbol:
> + eq_type = max(eq_type, GS_POLY)
> + elif item.is_Mul:
> + for arg in item.args:
> + eq_type = max(guess_solve_strategy(arg, symbol), eq_type)
> + elif item.is_Pow and item.base.has(symbol):
> + if item.exp.is_Integer:
> + if item.exp > 0:
> + eq_type = max(eq_type, GS_POLY)
> + else:
> + eq_type = max(eq_type, GS_POLY_CV_2)
> + elif item.exp.is_Rational:
> + eq_type = max(eq_type, GS_POLY_CV_1)
> + elif item.is_Function:
> + return GS_TRASCENDENTAL
> +
> + elif expr.is_Mul:
> + # check for rational functions
> + num, denom = expr.as_numer_denom()
> + if denom != 1 and denom.has(symbol):
> + #we have a quotient
> + m = max(guess_solve_strategy(num, symbol),
> guess_solve_strategy(denom, symbol))
> + if m == GS_POLY:
> + return GS_RATIONAL
> + elif m == GS_POLY_CV_1:
> + return GS_RATIONAL_CV_1
> + else:
> + raise NotImplementedError
> + else:
> + return max(map(guess_solve_strategy, expr.args, [symbol]))
> +
> + elif expr.is_Symbol:
> + return GS_POLY
> +
> + elif expr.is_Pow:
> + if expr.exp.has(symbol):
> + return GS_TRASCENDENTAL
> + elif expr.exp.is_Number and expr.base.has(symbol):
> + if expr.exp.is_Integer:
> + eq_type = max(eq_type, GS_POLY)
> + else:
> + eq_type = max(eq_type, GS_POLY_CV_1)
> +
> + elif expr.is_Function and expr.has(symbol):
> + return GS_TRASCENDENTAL
> +
> + return eq_type
> +
> def solve(f, *symbols, **flags):
> """Solves equations and systems of equations.
>
> @@ -81,13 +163,87 @@ def solve(f, *symbols, **flags):
> f = f.lhs - f.rhs
>
> if len(symbols) == 1:
> - poly = f.as_poly(*symbols)
> + symbol = symbols[0]
> +
> + strategy = guess_solve_strategy(f, symbol)
>
> - if poly is not None:
> + if strategy == GS_POLY:
> + poly = f.as_poly( symbol )
> + assert poly is not None
> result = roots(poly, cubics=True, quartics=True).keys()
> +
> + elif strategy == GS_RATIONAL:
> + P, Q = f.as_numer_denom()
> + #TODO: check for Q != 0
> + return solve(P, symbol, **flags)
> +
> + elif strategy == GS_POLY_CV_1:
> + # we must search for a suitable change of variable
> + # collect exponents
> + exponents_denom = list()
> + args = list(f.args)
> + if isinstance(f, Add):
> + for arg in args:
> + if isinstance(arg, Pow):
> + exponents_denom.append(arg.exp.q)
> + elif isinstance(arg, Mul):
> + for mul_arg in arg.args:
> + if isinstance(mul_arg, Pow):
> + exponents_denom.append(mul_arg.exp.q)
> + elif isinstance(f, Mul):
> + for mul_arg in args:
> + if isinstance(mul_arg, Pow):
> + exponents_denom.append(mul_arg.exp.q)
> +
> + assert len(exponents_denom) > 0
> + if len(exponents_denom) == 1:
> + m = exponents_denom[0]
> + else:
> + # get the GCD of the denominators
> + m = ilcm(*exponents_denom)
> + # x -> y**m.
> + # we assume positive for simplification purposes
> + t = Symbol('t', positive=True, dummy=True)
> + f_ = f.subs(symbol, t**m)
> + if guess_solve_strategy(f_, t) != GS_POLY:
> + raise TypeError("Could not convert to a polynomial
> equation: %s" % f_)
> + cv_sols = solve(f_, t)
> + result = list()
> + for sol in cv_sols:
> + result.append(sol**(S.One/m))
> +
> + elif strategy == GS_POLY_CV_2:
> + m = 0
> + args = list(f.args)
> + if isinstance(f, Add):
> + for arg in args:
> + if isinstance(arg, Pow):
> + m = min(m, arg.exp)
> + elif isinstance(arg, Mul):
> + for mul_arg in arg.args:
> + if isinstance(mul_arg, Pow):
> + m = min(m, mul_arg.exp)
> + elif isinstance(f, Mul):
> + for mul_arg in args:
> + if isinstance(mul_arg, Pow):
> + m = min(m, mul_arg.exp)
> + f1 = simplify(f*symbol**(-m))
> + result = solve(f1, symbol)
> + # TODO: we might have introduced unwanted solutions
> + # when multiplied by x**-m
> +
> + elif strategy == GS_TRASCENDENTAL:
> + #a, b = f.as_numer_denom()
> + # Let's throw away the denominator for now. When we have
> robust
> + # assumptions, it should be checked, that for the solution,
> + # b!=0.
> + result = tsolve(f, *symbols)
> + elif strategy == -1:
> + raise Exception('Could not parse expression %s' % f)
> else:
> - result = [tsolve(f, *symbols)]
> + raise NotImplementedError("No algorithms where implemented
> to solve equation %s" % f)
> else:
> + #TODO: if we do it the other way round we can avoid one nesting
> raise NotImplementedError('multivariate equation')
>
> if flags.get('simplified', True):
> @@ -530,10 +686,10 @@ def tsolve(eq, sym):
> >>> x = Symbol('x')
>
> >>> tsolve(3**(2*x+5)-4, x)
> - (-5*log(3) + log(4))/(2*log(3))
> + [(-5*log(3) + log(4))/(2*log(3))]
>
> >>> tsolve(log(x) + 2*x, x)
> - 1/2*LambertW(2)
> + [1/2*LambertW(2)]
>
> """
> if patterns is None:
> @@ -549,11 +705,11 @@ def tsolve(eq, sym):
> r = Wild('r')
> m = eq2.match((a*x+b)*r)
> if m and m[a]:
> - return (-b/a).subs(m).subs(x, sym)
> + return [(-b/a).subs(m).subs(x, sym)]
> for p, sol in patterns:
> m = eq2.match(p)
> if m:
> - return sol.subs(m).subs(x, sym)
> + return [sol.subs(m).subs(x, sym)]
>
> # let's also try to inverse the equation
> lhs = eq
> @@ -587,7 +743,7 @@ def tsolve(eq, sym):
> lhs = lhs.args[0]
>
> sol = solve(lhs-rhs, sym)
> - return sol[0]
> + return sol
>
> elif lhs.is_Add:
> # just a simple case - we do variable substitution for first function,
> @@ -610,18 +766,12 @@ def tsolve(eq, sym):
> # if no Functions left, we can proceed with usual solve
> if not (lhs_.is_Function or
> any(term.is_Function for term in lhs_.args)):
> -
> - # FIXME at present solve cannot solve x + 1/x = y
> - # FIXME so we do this:
> - numer, denom = lhs_.as_numer_denom()
> - sol = solve(numer-rhs*denom, t)
> - #sol = solve(lhs_-rhs, t)
> - sol = sol[0]
> -
> - sol = tsolve(sol-f1, sym)
> - return sol
> -
> -
> + cv_sols = solve(lhs_ - rhs, t)
> + cv_inv = solve( t - f1, sym )[0]
> + sols = list()
> + for sol in cv_sols:
> + sols.append(cv_inv.subs(t, sol))
> + return sols
>
>
> raise ValueError("unable to solve the equation")
> diff --git a/sympy/solvers/tests/test_solvers.py
> b/sympy/solvers/tests/test_solvers.py
> index 65c7f36..b5c09e8 100644
> --- a/sympy/solvers/tests/test_solvers.py
> +++ b/sympy/solvers/tests/test_solvers.py
> @@ -1,20 +1,68 @@
> from sympy import Matrix, Symbol, solve, exp, log, cos, acos, Rational, Eq, \
> sqrt, oo, LambertW, pi, I, sin, asin, Function, diff, Derivative, \
> - symbols, S, raises
> + symbols, S, raises, sympify, var, simplify
> from sympy.solvers import solve_linear_system,
> solve_linear_system_LU,dsolve,\
> tsolve, deriv_degree
>
> +from sympy.solvers.solvers import guess_solve_strategy, GS_POLY,
> GS_POLY_CV_1, GS_TRASCENDENTAL, \
> + GS_RATIONAL, GS_RATIONAL_CV_1
> +
> from sympy.utilities.pytest import XFAIL
>
> -def test_solve():
> +def test_guess_strategy():
> + """
> + See solvers._guess_solve_strategy
> + """
> + x, y = symbols('xy')
> +
> + # polynomial equations
> + assert guess_solve_strategy( x, x ) == GS_POLY
> + assert guess_solve_strategy( 2*x, x ) == GS_POLY
> + assert guess_solve_strategy( x + sqrt(2), x) == GS_POLY
> + assert guess_solve_strategy( x + 2**Rational(1,4), x) == GS_POLY
> + assert guess_solve_strategy( x**2 + 1, x ) == GS_POLY
> + assert guess_solve_strategy( x**2 - 1, x ) == GS_POLY
> + assert guess_solve_strategy( x*y + y, x ) == GS_POLY
> + assert guess_solve_strategy( x*exp(y) + y, x) == GS_POLY
> + assert guess_solve_strategy( (x - y**3)/(y**2*(1 - y**2)**(1/2)), x) ==
> GS_POLY
> +
> + # polynomial equations via a change of variable
> + assert guess_solve_strategy( x**Rational(1,2) + 1, x ) == GS_POLY_CV_1
> + assert guess_solve_strategy( x**Rational(1,3) + x**Rational(1,2) + 1, x
> ) == GS_POLY_CV_1
> +
> + # polynomial equation multiplying both sides by x**n
> + assert guess_solve_strategy( x + 1/x + y, x )
> +
> + # rational functions
> + assert guess_solve_strategy( (x+1)/(x**2 + 2), x) == GS_RATIONAL
> + assert guess_solve_strategy( (x - y**3)/(y**2*(1 - y**2)**(1/2)), y) ==
> GS_RATIONAL
> +
> + # rational functions via the change of variable y -> x**n
> + assert guess_solve_strategy( (x**Rational(1,2) + 1)/(x**Rational(1,3) +
> x**Rational(1,2) + 1), x ) \
> + == GS_RATIONAL_CV_1
> +
> + #trascendental functions
> + assert guess_solve_strategy( exp(x) + 1, x ) == GS_TRASCENDENTAL
> + assert guess_solve_strategy( 2*cos(x)-y, x ) == GS_TRASCENDENTAL
> + assert guess_solve_strategy( exp(x) + exp(-x) - y, x ) ==
> GS_TRASCENDENTAL
> +
> +def test_solve_polynomial():
> x, y = map(Symbol, 'xy')
>
> +
> assert solve(3*x-2, x) == [Rational(2,3)]
> assert solve(Eq(3*x, 2), x) == [Rational(2,3)]
>
> assert solve(x**2-1, x) in [[-1, 1], [1, -1]]
> assert solve(Eq(x**2, 1), x) in [[-1, 1], [1, -1]]
>
> + assert solve( x - y**3, x) == [y**3]
> + assert solve( x - y**3, y) == [
> + (-x**Rational(1,3))/2 + I*sqrt(3)*x**Rational(1,3)/2,
> + x**Rational(1,3),
> + (-x**Rational(1,3))/2 - I*sqrt(3)*x**Rational(1,3)/2
> + ]
> +
> a11,a12,a21,a22,b1,b2 = symbols('a11','a12','a21','a22','b1','b2')
>
> assert solve([a11*x + a12*y - b1, a21*x + a22*y - b2], x, y) == \
> @@ -27,9 +75,40 @@ def test_solve():
> assert solve((x-y, x+y), (x, y)) == solution
> assert solve((x-y, x+y), [x, y]) == solution
>
> + assert solve( x**3 - 15*x - 4, x) == [-2 + 3**Rational(1,2),
> + 4,
> + -2 - 3**Rational(1,2) ]
> +
> raises(TypeError, "solve(x**2-pi, pi)")
> raises(ValueError, "solve(x**2-pi)")
>
> +def test_solve_polynomial_cv_1():
> + """
> + Test for solving on equations that can be converted to a polynomial
> equation
> + using the change of variable y -> x**Rational(p, q)
> + """
> +
> + x = Symbol('x')
> +
> + assert solve( x**Rational(1,2) - 1, x) == [1]
> + assert solve( x**Rational(1,2) - 2, x) == [sqrt(2)]
> +
> +def test_solve_polynomial_cv_2():
> + """
> + Test for solving on equations that can be converted to a polynomial
> equation
> + multiplying both sides of the equation by x**m
> + """
> +
> + x = Symbol('x')
> +
> + assert solve( x + 1/x - 1, x) == [Rational(1,2) + I*sqrt(3)/2,
> Rational(1,2) - I*sqrt(3)/2]
> +
> +def test_solve_rational():
> + x = Symbol('x')
> + y = Symbol('y')
> +
> + solve( ( x - y**3 )/( (y**2)*sqrt(1 - y**2) ), x) == [x**Rational(1,3)]
> +
> def test_linear_system():
> x, y, z, t, n = map(Symbol, 'xyztn')
>
> @@ -101,35 +180,37 @@ def test_tsolve():
> assert solve(cos(x)-y, x) == [acos(y)]
> assert solve(2*cos(x)-y,x)== [acos(y/2)]
> # XXX in the following test, log(2*y + 2*...) should -> log(2) + log(y
> +...)
> - assert solve(exp(x)+exp(-x)-y,x) == [-log(4) + log(2*y + 2*(-4 +
> y**2)**Rational(1,2))]
> - assert tsolve(exp(x)-3, x) == log(3)
> - assert tsolve(Eq(exp(x), 3), x) == log(3)
> - assert tsolve(log(x)-3, x) == exp(3)
> - assert tsolve(sqrt(3*x)-4, x) == Rational(16,3)
> - assert tsolve(3**(x+2), x) == -oo
> - assert tsolve(3**(2-x), x) == oo
> - assert tsolve(4*3**(5*x+2)-7, x) ==
> (log(Rational(7,4))-2*log(3))/(5*log(3))
> - assert tsolve(x+2**x, x) == -LambertW(log(2))/log(2)
> + assert solve(exp(x)+exp(-x)-y,x) == [-log(4) + log(2*y + 2*(-4 +
> y**2)**Rational(1,2)),
> + -log(4) + log(2*y - 2*(-4 +
> y**2)**Rational(1,2))]
> + assert tsolve(exp(x)-3, x) == [log(3)]
> + assert tsolve(Eq(exp(x), 3), x) == [log(3)]
> + assert tsolve(log(x)-3, x) == [exp(3)]
> + assert tsolve(sqrt(3*x)-4, x) == [Rational(16,3)]
> + assert tsolve(3**(x+2), x) == [-oo]
> + assert tsolve(3**(2-x), x) == [oo]
> + assert tsolve(4*3**(5*x+2)-7, x) ==
> [(log(Rational(7,4))-2*log(3))/(5*log(3))]
> + assert tsolve(x+2**x, x) == [-LambertW(log(2))/log(2)]
> assert tsolve(3*x+5+2**(-5*x+3), x) == \
> - -Rational(5,3) +
> LambertW(-10240*2**Rational(1,3)*log(2)/3)/(5*log(2))
> + [-Rational(5,3) +
> LambertW(-10240*2**Rational(1,3)*log(2)/3)/(5*log(2))]
> assert tsolve(5*x-1+3*exp(2-7*x), x) == \
> - Rational(1,5) + LambertW(-21*exp(Rational(3,5))/5)/7
> + [Rational(1,5) + LambertW(-21*exp(Rational(3,5))/5)/7]
> assert tsolve(2*x+5+log(3*x-2), x) == \
> - Rational(2,3) + LambertW(2*exp(-Rational(19,3))/3)/2
> - assert tsolve(3*x+log(4*x), x) == LambertW(Rational(3,4))/3
> - assert tsolve((2*x+8)*(8+exp(x)), x) == -4
> - assert tsolve(2*exp(3*x+4)-3, x) == -Rational(4,3)+log(Rational(3,2))/3
> - assert tsolve(2*log(3*x+4)-3, x) == (exp(Rational(3,2))-4)/3
> - assert tsolve(exp(x)+1, x) == pi*I
> - assert tsolve(x**2 - 2**x, x) == 2
> - assert tsolve(x**3 - 3**x, x) == -3/log(3)*LambertW(-log(3)/3)
> + [Rational(2,3) + LambertW(2*exp(-Rational(19,3))/3)/2]
> + assert tsolve(3*x+log(4*x), x) == [LambertW(Rational(3,4))/3]
> + assert tsolve((2*x+8)*(8+exp(x)), x) == [-4]
> + assert tsolve(2*exp(3*x+4)-3, x) == [-Rational(4,3)+log(Rational(3,2))/3]
> + assert tsolve(2*log(3*x+4)-3, x) == [(exp(Rational(3,2))-4)/3]
> + assert tsolve(exp(x)+1, x) == [pi*I]
> + assert tsolve(x**2 - 2**x, x) == [2]
> + assert tsolve(x**3 - 3**x, x) == [-3/log(3)*LambertW(-log(3)/3)]
> assert tsolve(2*(3*x+4)**5 - 6*7**(3*x+9), x) == \
> - Rational(-4,3) -
> 5/log(7)/3*LambertW(-7*2**Rational(4,5)*6**Rational(1,5)*log(7)/10)
> + [Rational(-4,3) -
> 5/log(7)/3*LambertW(-7*2**Rational(4,5)*6**Rational(1,5)*log(7)/10)]
>
> - assert tsolve(z*cos(x)-y, x) == acos(y/z)
> - assert tsolve(z*cos(2*x)-y, x) == acos(y/z)/2
> - assert tsolve(z*cos(sin(x))-y, x) == asin(acos(y/z))
> + assert tsolve(z*cos(x)-y, x) == [acos(y/z)]
> + assert tsolve(z*cos(2*x)-y, x) == [acos(y/z)/2]
> + assert tsolve(z*cos(sin(x))-y, x) == [asin(acos(y/z))]
>
> - assert tsolve(z*cos(x), x) == acos(0)
> + assert tsolve(z*cos(x), x) == [acos(0)]
>
> - assert tsolve(exp(x)+exp(-x)-y, x)== log(y/2 + Rational(1,2)*(-4 +
> y**2)**Rational(1,2))
> + assert tsolve(exp(x)+exp(-x)-y, x)== [log(y/2 + Rational(1,2)*(-4 +
> y**2)**Rational(1,2)),
> + log(y/2 - Rational(1,2)*(-4 +
> y**2)**Rational(1,2))]
All patches are +1.
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy-patches" group.
To post to this group, send email to sympy-patches@googlegroups.com
To unsubscribe from this group, send email to [EMAIL PROTECTED]
For more options, visit this group at
http://groups.google.com/group/sympy-patches?hl=en
-~----------~----~----~----~------~----~------~--~---
