Applications of Algorithms for Solving Equations and Inequalities in Mathematica

Adam Strzebonski

Introduction

Computational problems arising in many areas of a computer algebra system can be formulated as questions about solution sets of systems of equations and inequalities. The ability to answer such questions is therefore a very important factor contributing to the computational strength of the system. Mathematica contains an implementation of algorithms for solving many kinds of systems of equations and inequalities over various domains. Several of these algorithms have been invented or improved by Volker Weispfenning. The algorithms are used directly in the Mathematica functionality for solving systems of equations and inequalities over complex, real and integer domains, as well as for finding solution instances, quantifier elimination, and exact optimization. They are also used behind the scenes in many other parts of the system. The assumption mechanism uses decision algorithms for proving inference. Assumptions are used in simplification, computation of limits, series, integrals and integral transforms. Real system solving is used in symbolic and numeric integration over regions. It allows to find explicit integration bounds from implicit region descriptions. Mathematica 5.1 contains system-wide support for piecewise functions. Simplification and computation of limits, series, derivatives and integrals of piecewise functions, solving differential equations involving piecewise functions as well as visualization of piecewise functions all use equation and inequality solving algorithms for investigating geometry of piecewise conditions.

Solving Systems of Equations and Inequalities

Mathematica provides functionality for finding explicit descriptions of solution sets of equations and inequalities.

Solving of this system involves using a heuristic for solving transcendental systems, as well as linear equation solving, real system solving using CAD, and binary quadratic Diophantine equation solving.

Reduce[
      Sin[π x + π y] == (1/2) && 
        Cos[π x - π y] == ((Sqrt[3])/2) && -72 x + 36 (x^2) + 
            36 y - 144 x y + 144 (y^2) == 25 && 
        (x^2) + (y^2) ≤ (10^6), {x, y}, Reals]

Solving of this system involves using a heuristic for solving transcendental systems and Groebner basis computation.

Reduce[
      Sinh[x] == a && (E^x) == (b^2) + (a^3) + 1, 
      x]

Methods

Complex and real systems of linear equations are solved using linear algebra methods. Complex polynomial systems are solved using Groebner basis methods. To solve real algebraic systems and mixed real-complex polynomial systems the cylindrical algebraic decomposition (CAD) algorithm is used. There is no general method for solving Diophantine systems, but there are algorithms for solving many large special cases. Mathematica contains implementation of several such algorithms. For systems involving transcendental functions heuristics are used for decomposing the systems into algebraic systems composed with simple transcendental systems that can be solved using inverse image information for transcendental functions.

Finding Solution Instances, Decision Problems and Assumptions

This uses Loos-Weispfenning linear quantifier elimination algorithm to eliminate x, then the quadratic case of Weispfenning's virtual substitution algorithm to eliminate y, then a simplified strict inequality version of CAD to get a solution for z, t, and u.

The quadratic case of Weispfenning's virtual substitution algorithm is implemented in the development version of Mathematica which I used at the conference. Mathematica 5.1, which was used to generate outputs in this notebook, uses here a combination of linear quantifier elimination and strict inequality CAD.

FindInstance[
      2 x + (y^2) z - t u y ≥ 2 && 
        (z^4) + (t^4) + (u^4) == (y^2) - t z u y + 1, {x, y, z, t, u}, 
      Reals]

Algorithms for solving decision problems are used in the assumption mechanism.

expr = (Sqrt[x (((x - (y^2) - 1))^2) (((z - 2 x))^4)])

Refine[expr, x ≥ 0]

Refine[expr, 0 < x < y < (z/2)]

To compute this limit we need to know that |x|<1.

Limit[(x^n), n -> ∞]

Limit[(x^n), n -> ∞, 
      Assumptions -> Abs[x] < 1]

The assumption mechanism uses decision algorithms to prove that  Re[x]^4 + Im[x]^4 < (1/2) implies |x|<1.

Limit[(x^n), n -> ∞, 
      Assumptions -> Re[x]^4 + Im[x]^4 < (1/2)]

Assumptions are used to determine validity of special function transformations. This function can be simplified if -(π/2) < Re[x] < (π/2).

FunctionExpand[EllipticF[x,  1]]

FunctionExpand[EllipticF[x, 1], (x^2) < (E/2)]

Methods

The general algorithms for solving systems of equations and inequalities can be used to decide existence of solutions or to find a single solution instance, and most of them can be significantly simplified when used for this purpose. The case of real polynomial systems is of special importance here, since real inequalities are the most commonly used assumptions. There are several special case algorithms applicable to real systems.

For linear systems Mathematica uses a method based on the simplex algorithm. The Loos-Weispfenning linear quantifier elimination algorithm is used for linear systems with non-rational coefficients and for systems which contain some linear variables. In the latter case the algorithm is used to eliminate linear variables with constant coefficients, as long as there are any, and the resulting nonlinear system is handled by CAD. The quadratic case of Weispfenning's virtual substitution algorithm is used for sytems with at least four variables. This choice of default methods was made based on experiments, as a most often, but not always, fastest method combination. Other method combinations can be specified with options.

When a real polynomial system contains only strict inequalities a simplified version of CAD is used. Two other methods, by Aubry, Rouillier, and Safey El Din for real systems containing only equations and by G.X. Zeng and X.N. Zeng for systems consisting of a single strict inequality, are implemented. While there are examples for which these two methods are significantly faster than CAD-based methods, they seem to be slower for randomly selected systems. It is not clear how to characterize systems for which these methods are faster, and therefore they are used only if requested through a non-default option setting.

Full Dimensional Real Solving, Integration and Visualization

When computing, symbolically or numerically, integrals over regions or integrals of piecewise defined functions, or when plotting piecewise functions, it is very advantageous (or necessary in the symbolic integration case) to be able to obtain explicit functional bounds for variable values within regions.

This computes the volume of the intersection of a ball and a double cone.

(∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞) Boole[
                (x^2) + (y^2) + (z^2) ≤ 1 && 
                  3 (x^2) + 3 (y^2) ≤ 
                    (z^2)] \[DifferentialD]z \[DifferentialD]y \[DifferentialD]x  // Timing

Full dimensional CAD is used to divide the set into cylindrical regions with boundaries given by algebraic functions.

Simplify`PiecewiseIntegrate; 
    PiecewiseDump`$ShowIntegrals = 
      True; (∫_(-∞)^∞)(∫_(-∞)^∞)(∫_(-∞)^∞) Boole[
                (x^2) + (y^2) + (z^2) ≤ 1 && 
                  3 (x^2) + 3 (y^2) ≤ 
                    (z^2)] \[DifferentialD]z \[DifferentialD]y \[DifferentialD]x  // Timing

PiecewiseDump`$ShowIntegrals = False;

Here CAD is used to get a parametrized cylindrical description of the integration region. The result is a piecewise function of the parameter.

(∫_0^1)(∫_0^1) 
            x (E^y) Boole[
              a (x^2) + (y^2) < 
                1] \[DifferentialD]y \[DifferentialD]x

This uses transcendental inequality solver to partition the integration region into countably many intervals. Symbolic summation is then used to get the answer. This method is valid only if the integral is convergent. A real decision algorithm is used to prove that the integrand is nonnegative.

(∫_(-∞)^∞) Boole[
              (1/2) < Sin[x] < (2/3)]/(x^2) \[DifferentialD]x

Here numeric integration uses full dimensional CAD to integrate a piecewise function.

NIntegrate[(Piecewise[{{(y^4), (x^4) (y^2) + (y^2) ≤ 1}, {(1/y), True}}]),{x, -10, 10},{y, -10, 10}]

Without symbolic preprocessing the numeric integration fails.

f[x_?NumericQ, y_?NumericQ]:=(Piecewise[{{(y^4), (x^4) (y^2) + (y^2) ≤ 1}, {(1/y), True}}])

NIntegrate[f[x, y], {x, -10, 10}, {y, -10, 10}]

The following graphics examples need the development version of Mathematica. Symbolic preprocessing of Plot3D inputs is not implemented in Mathematica 5.1. The graphics that follow have been produced with the development version and copied as bitmaps.

Here Plot3D uses full dimensional CAD to decompose the domain of a piecewise function into parts on which the piecewise conditions have constant values.

Plot3D[
      Boole[(((x + (7/4)))^2) + (((y - (9/4)))^2) < 
            (1/9) || (((x - (9/4)))^2) + (((y - (9/4)))^2) < 
            (1/9) || ((x^2) + (y^2) < 25 && (((x - 2))^2) + (((y - 2))^2) > 
                1 && (((x + 2))^2) + (((y - 2))^2) > 
                1 && ((x^2) + (((y + 1))^2) > 5 || 
                  (x^2) + (((y - 1))^2) < 9))], {x, -5, 
        5}, {y, -5, 5}, BoundaryStyle -> None, 
      ColorFunction -> (Hue[#(3/2) + (1/4)] &)]

Without the symbolic preprocessing Plot3D cannot identifiy the case boundaries that well.

Plot3D[
      Boole[(((x + (7/4)))^2) + (((y - (9/4)))^2) < 
            (1/9) || (((x - (9/4)))^2) + (((y - (9/4)))^2) < 
            (1/9) || ((x^2) + (y^2) < 25 && (((x - 2))^2) + (((y - 2))^2) > 
                1 && (((x + 2))^2) + (((y - 2))^2) > 
                1 && ((x^2) + (((y + 1))^2) > 5 || 
                  (x^2) + (((y - 1))^2) < 9))], {x, -5, 
        5}, {y, -5, 5}, BoundaryStyle -> None, 
      ColorFunction -> (Hue[#(3/2) + (1/4)] &), 
      Method -> {"SymbolicPiecewiseSubdivision" -> 
            False}]

Methods

For polynomial systems the CAD algorithm can be used to decompose regions into finitely many parts, each bounded by explicit algebraic functions. In fact, since measure zero sets do not matter here, having an explicit description of a set which differs from the solution set by lower dimensional components is quite sufficient, hence a simplified version of CAD can be used. For symbolic integration over regions described using transcendental inequalities, full description of the solution set is computed,  using the methods described before, and if the result contains infinitely many parts, symbolic summation is used to combine the partial integrals.

Exact Global Optimization

Mathematica uses real equation and inequality solving algorithms for solving exact global optimization problems.

This finds the point on the cubic curve (x^3) + (y^2) - x == (1/4) which is closest to the origin. This computation uses CAD.

Minimize[{(x^2) + (y^2), (x^3) + (y^2) - x == (1/4)}, {x, 
        y}]

This uses the Lagrange multiplier method and real equation and inequality solving to find the stationary and singular points.

Minimize[{y + Sin[5 x], 
        (y^3) == Cos[5 x] && -5 <= x <= 5 && -5 <= y <= 5}, {x, 
        y}]

N[%,  20]

Here Mathematica uses exact global optimization to prove that a composition of piecewise functions Mod, Ceiling and Floor, with a potentially infinite number of cases has in fact just five piecewise cases. (The case with value zero is suppressed in the output.)

PiecewiseExpand[
      Mod[(x^2), Ceiling[x] + 1] Boole[2  x - Floor[x]^2 > 0]]

Methods

For linear systems the simplex algorithm is used.

For univariate problems equation and inequality solving methods are used to find the constraint solution set and discontinuity points and zeros of the derivative of the objective function within the set.

When the inequality constraints define a bounded box on which the objective function is (C^1), equation and inequality solving methods are used to find singular points of the set described by the equational constraints and stationary points of the objective function on the set described by the equational constraints. This computation is done within the bounding box, and on its faces of all dimensions. If the number of singular and stationary points is finite, the extremum must be attained at one of them.

General polynomial optimization problems are solved using an algorithm based on CAD. The problem of minimizing an arbitrary polynomial (or rational function) f is reduced to minimizing the value of a new variable y on the semialgebraic set described by the constraints and by y = f((Subscript[x, 1]),  ... , (Subscript[x, n])). The new variable y is set to be the last variable in the projection order. The cell construction phase of CAD starts with the lowest values of y and is stopped as soon as a cell satisfying the constraints is found. If all constraints are strict inequalities a simplified projection operator is used.

Quantifier Elimination and Piecewise Calculus

Mathematica provides functions for eliminating quantifiers from systems of equations and inequalities. For instance, we can use the definition of limit to compute the limit of a rational function.

Resolve[ (Subscript[∀ , ϵ, ϵ > 0])(Subscript[∃ , M])(Subscript[∀ , x, x > M])-ϵ < ((3 (x^2) - 5 x + 7)/(7 (x^2) - 2 x + 11)) - λ < ϵ, Reals]

Whenever possible Resolve eliminates quantifiers without solving the resulting system. Here it eliminates the quantifier using the Weispfenning virtual substitution method.

The quadratic case of Weispfenning's virtual substitution algorithm is implemented in the development version of Mathematica which I used at the conference. Mathematica 5.1 would use CAD here, just like it does in the next example. The answer and timing below have been generated with the development version.

Resolve[(Subscript[∃ , x])(a (x^2) + b x + c) (c (x^2) + b x + 
                  a) ≤ 0, Reals] // Timing

Quantified systems can be also given to Reduce which eliminates quantifiers and solves the resulting system.

Reduce[(Subscript[∃ , x])(a (x^2) + b x + c) (c (x^2) + b x + 
                  a) ≤ 0, {a, b, c}, Reals] // Timing

Mathematica uses quantifier elimination internally in piecewise function calculus. For instance to compute a limit of a piecewise function with \(TraditionalForm\`x -> \(0\^+\)\) one needs to find a piecewise condition c for which values of parameters \(TraditionalForm\`∃ \_\(ϵ > 0\)\(∀ \_\(0 < x < ϵ\)\ c(x)\)\).

Here quantifier elimination is used to to prove that there is an ϵ>0 for which the interval (i+ϵ, i) is outside the unit circle.

Limit[(Piecewise[{{1, Abs[z] ≤ 1}, {2, True}}]),z -> I]

There is an ϵ>0 for which the interval (i - ϵ (e^(((i π)/6))), i) is inside the unit circle.

Limit[(Piecewise[{{1, Abs[z] ≤ 1}, {2, True}}]),z -> I,Direction -> (E^(((I π)/6)))]

Here quantifier elimination is used to find the value of Max[a (x^2) + (x/2) + (a^2), ((x^2)/2) + a x + 1] in a neigbourhood of 0. The answer depends on the parameter a.

Series[(E^Max)[a (x^2) + (x/2) + (a^2), ((x^2)/2) + a x + 1], {x, 0, 2}, Assumptions -> (a | x) ∈ Reals]

This experimental functionality uses quantifier elimination over (Subscript[ℤ, 7]) to find a counterexample to the cnjecture that the axiom \(TraditionalForm\`∀ \_{a, b}f(f(f(b, b), a), b) == 
          a\)   implies the commutativity of f.

The following example requires the development version of Mathematica.

\(EquationalLogic`FindCounterexample[
      f[a, b] == 
        f[b, a], ∀ \_{a, b}f[f[f[b, b], a], b] == a]\)

FullSimplify[f[f[f[b, b], a], b] == a /. %[[1]],  
      Element[a | b,  Integers] && 0 ≤ a < 7]

f[a, b] == 
          f[b, a] /. %%[[1]] /. %%[[2]]

Methods

To eliminate quantifiers from complex polynomial systems Groebner basis methods are used. For real algebraic systems and mixed real-complex polynomial systems methods based on the CAD algorithm are used. The Loos-Weispfenning linear quantifier elimination algorithm is used to eliminate linear variables from the innermost quantifiers, as long as there are any such variables. If there are quantifiers left, the resulting nonlinear system is handled by CAD. The quadratic case of Weispfenning's virtual substitution algorithm is used when we do not need result in a solved form or there are at least three variables and the size of the system does not grow too much after the elimination. This choice of default methods was made based on experiments, as a most often, but not always, fastest method combination. Other method combinations can be specified with options.

Conclusions

Algorithms for solving systems of equations and inequalities have applications in multiple areas of symbolic and numeric computation. Having them available within a computer algebra system allows for development of many new features. The algorithms provide computational engine used for instance in such seemingly unrelated areas as numeric integration, special function simplification, or piecewise function calculus and visualization.

Bibliography

P. Aubry, F. Rouillier,and M. Safey El Din. Real solving for positive dimensional systems. J. Symb. Comput., 34:543--560, 2002.

B. Buchberger. An algorithmic criterion for the solvability of algebraic systems of equations. Aequationes Mathematicae, 4:374--383, 1970.

T. Becker and V. Weispfenning. Groebner bases-A computational approach to commutative algebra. Springer Verlag, New York, 1993.

B. Caviness and J. Johnson,editors. Quantifier Elimination and Cylindrical Algebraic Decomposition, New York, 1998. Springer Verlag.

G. E. Collins. Quantifier elimination for the elementary theory of real closed fields by cylindrical algebraic decomposition. Lect. Notes Comput. Sci., 33:134--183, 1975.

G. B. Dantzig. Linear Programming and Extensions. Princeton University Press, 1963.

R. Loos and V. Weispfenning. Applying linear quantifier elimination. The Computer Journal, 36(5):450--462, 1993.

Yu. V. Matiyasevich. The diophantineness of enumerable sets. Soviet Math.Dokl., 11:354--358, 1970.

A. Strzebonski. Solving systems of strict polynomial inequalities. J. Symbolic Comp., 29:471--480, 2000.

A. Strzebonski. Advanced documentation: Equations and inequalities. Mathematica 5.1 online help browser, Wolfram Research, Inc., 2004.

V. Weispfenning. Quantifier elimination for real algebra-the cubic case. In ISSAC 94, 258--263, 1994.

V. Weispfenning. Quantifier elimination for real algebra-the quadratic case and beyond. Appl. Algebra Eng. Commun. Comput., 8(2):85--101, 1997.

S. Wolfram. The Mathematica Book. Wolfram Media, 5th edition, 2003.

G. X. Zeng and X. N. Zeng. An effective decision method for semidefinite polynomials.J. Symb. Comput., 37:83--99, 2004.