Solving Equations and Inequalities

Adam Strzebonski

Wolfram Research, Inc.

Extended functionality of Solve and NSolve

In Mathematica 8 Solve and NSolve have been extended to allow the same input syntax as Reduce, including inequalities, domain specifications and quantifiers.

Solve[expr,vars] attempts to solve the system expr of equations or inequalities for the variables vars. 
Solve[expr,vars,dom] solves over the domain dom. Common choices of dom are Reals, Integers and Complexes.
NSolve[expr,vars] attempts to find numerical approximations to the solutions of the system expr of equations or inequalities for the variables vars.
NSolve[expr,vars,Reals] finds solutions over the domain of real numbers.

Solve examples

Solve a system of equations over the reals.

Solve[(x^3)+(y^3)+(z^3)==1&&3x-2y==-z&&x+ y+ z==3,{x,y,z},Reals]

Find a condition for a quartic to have all roots equal.

Solve[(Subscript[∃, x,(x^4)+a (x^3)+b (x^2)+c x+d==0])((Subscript[∀, y,(y^4)+a (y^3)+b (y^2)+c y+d==0])x==y),{b,c,d}]

Solve a Diophantine equation.

Solve[(x^3)+3x (y^2)+(y^3)==811,{x,y},Integers]

NSolve examples

Find approximations of real roots of a polynomial.

NSolve[(x^1000)+1.23 (x^777)-4.56 (x^345)+7.89 x-1==0,Reals]

Find approximate solutions of univariate transcendental equations.

NSolve[(E^x)-x==7,x,Reals]

NSolve[2 Sin[(E^(2x))]-Cos[2π x]==(3/2)&&-1<x<1,x,Reals]

NSolve[Gamma[x]-Log[x]==(I/2)&&Abs[x-2]<(3/2),x]

Systems of real polynomial equations and inequalities.

NSolve[x y==z^2-x&&x y z==2&&x^2+y^2+z^2==7&&x<1,{x,y,z},Reals]

NSolve[Exists[x, x^2+a x+b==0&&2 x+a==0&&a^2+x b^2==1],{a,b},Reals]

Systems of transcendental equations and inequalities.

NSolve[Sin[x+y]==1/2&&E^x-y==1&&x^2≤10,{x,y},Reals]

NSolve[Sin[x-Cos[x]] y^3-x==1&&x^2+y^2==1,{x,y},Reals]

Representation of solutions as replacement rules

Representing solutions of equations may require non-equational conditions on parameters.

e=(x^2)+x y+(y^2);
ContourPlot[e==1,{x,-2,2},{y,-2,2}]

Non-equational conditions on parameters are given using new ConditionalExpression constructs.

sol=Solve[e==1,y,Reals]

Conditions on parameters are preserved when solutions are used for replacement.

e/.sol

Expand[%]

cexpr=y/.sol[[1]]

ConditionalExpression is defined when the condition is satisfied.

cexpr/.x->1/2

Plot[cexpr,{x,-2,2}]

cexpr/.x->2

Generic solutions

Solve gives generic solutions only. Solutions that require continuous parameters to satisfy inequations are given without explicit conditions and solutions that require continuous parameters to satisfy equations are not given.

This implicitly assumes that a≠0.

Solve[a (x^5)+b (x^2)+c x+1==0,x]

With MaxExtraConditions->k, Solve gives solutions that require at most k equational conditions on continuous parameters. Inequation conditions required for validity of solutions are not explicitly given.

Solve[a (x^5)+b (x^2)+c x+1==0,x,MaxExtraConditions->1]

Solve[a (x^5)+b (x^2)+c x+1==0,x,MaxExtraConditions->2]

With MaxExtraConditions->All, Solve gives all solutions and includes all conditions required for their validity.

Solve[a (x^5)+b (x^2)+c x+1==0,x,MaxExtraConditions->All]

Equational conditions on discrete parameters are allowed by default.

Solve[x^2+y^2≤2,y,Integers]

“Syntactic assumptions” feature removed

In previous versions of Mathematica Solve treated equations not containing Solve variables as “assumptions”.

This system has a solution only if a==b, hence it has no generic solutions.

Solve[x==a&&x==b,x]

In Mathematica 7 and earlier, Solve would return a solution if an extra equation  a==b was added.

Solve[x==a&&x==b&&a==b,x]

The feature was based on syntactic properties of the input equations. With a mathematically equivalent equation  a+x==b+x added, no solution was returned.

Solve[x==a&&x==b&&a+x==b+x,x]

This feature was not well-defined mathematically and will not be present in Mathematica 8. MaxExtraConditions option can be used instead.

Solve[x==a&&x==b,x,MaxExtraConditions->1]

Normal can be used to remove ConditionalExpression.

Normal[%]

Multiple solutions

The new Solve finds solution multiplicity only for univariate equations.

Solve[(((x-1))^2) (((x-2))^3)==0,x]

Solution multiplicities are given for non-polynomial analytic equations as well.

Solve[Cos[x]-1==0 && -1<x<1,x]

Solve[(E^((((E^x)-Log[x]-3))^3))==1,x,Reals]

Extended functionality of InverseFunction

In  Mathematica 8 InverseFunction will try to compute closed from inverses of pure functions.

InverseFunction[((a #1+b)/(c #1+d))&]

If the function is not one-to-one, InverseFunction gives an error message and picks one of the solutions.

InverseFunction[((#1)^2)+2&]

Over the positive reals, the function is one-to-one.

InverseFunction[ConditionalExpression[((#1)^2)+2,#1>0]&]

InverseFunction of a function with restricted domain computes the domain of the inverse.

InverseFunction[ConditionalExpression[Sin[((#1)^2)-(1/2)],0<#1<1]&]

InverseFunction wrt. the second argument.

InverseFunction[((#1)^2)+((#2)^2)-1&,2,2]

There is no closed form formula for this InverseFunction.

f=InverseFunction[Sin[#]-#&]

Exact point values can be computed.

f[2]

However, the inverse may not be unique.

Reduce[Sin[x]-x==2&&Abs[x]<5,x]

InverseFunction wrt. the first argument of a two-argument function.

g=InverseFunction[Function[{x,y},Beta[x,y]-x+y],1,2]

g[3,4]

Here we cannot find a closed form formula for the InverseFunction or an exact point value.

h=InverseFunction[(∑_k^∞)(k^(-#1 k))&]

h[3]

But we can compute an arbitrarily accurate approximation of the point value.

N[h[3],50]

Cylindrical algebraic formulas

A semialgebraic set is the solution set of a system of polynomial equations and inequalities.

Each semialgebraic set can be represented as a finite union of cylindrically arranged cells bounded by graphs of algebraic functions.

ineqs=-5 (((y-3))^3)+7 (x^2) (y-3)<1&&(x^2)+2 (((y-4))^2)<21;

pl=RegionPlot[ineqs,{x,-4.5,4.5},{y,-0.5,7.75},AspectRatio->Automatic,PlotPoints->100];
cad=CylindricalDecomposition[ineqs,{x,y}];
bds=Select[cad/.{Or->List,And->List,HoldPattern[Inequality[a_,__,b_]]:>{a,b}},ListQ[#[[1]]]&];
pts=Sort[Union@@(First/@bds),Less];
cbs=Table[{Graphics[{c=Hue[(i-1)/7],Thickness[0.005],Line[{{pts[[i]],0},{pts[[i+1]],0}}]}],Plot[bds[[i,2]],{x,bds[[i,1,1]],bds[[i,1,2]]},ColorFunction->(c&),PlotStyle->Thickness[0.005],PlotRange->{{-4.5,4.5},{-0.5,7.75}}]},{i,7}];
Show[{pl,cbs,Graphics[{Thickness[0.002],Gray,Line[{{#,-0.5},{#,7.75}}]&/@pts}]},Axes->True,Frame->False]

Solution sets of polynomial equations and inequalities computed by Reduce (or CylindricalDecomposition) are given as finite unions of cylindrically arranged cells, represented by cylindrical algebraic formulas.

CAFView[CylindricalDecomposition[(x^2)+(y^2)+(z^2)≤1,{x,y,z}],{x,y,z}]

The cylindrical algebraic formula representation of a semialgebraic set A can be used to:

• decide whether A is nonempty,

• find the minimal and maximal values of the first coordinate of elements of A,

• generate an arbitrary element of A,

• find a graphical representation of A,

• compute the volume of A,

• compute multidimensional integrals over A,

• ...

CAFView code

cview[cad_,vars_,rul_,dvars_,ind_]:=
Module[{k=Length[ind],cc,ee,a,b,ans=False,ha,hb,x,na,rr},
ha=Hue[k/Length[vars],1/3];
hb=Hue[(k+1)/Length[vars],1/3];
If[rul=!={},
a=cad/.rul;
If[a===True||a===False,Return[Row[{a},Background->ha,Frame->True]]]];
x=dvars[[k+1]];
cc=If[Head[cad]===Or,List@@cad,{cad}];
If[k==Length[vars]-1,
Return[Or@@(Row[{Tooltip[#,N[#/.rul]]},Background->ha,Frame->True]&/@cc)]];
ee=cc[[1]];
If[Head[ee]===And,ee=ee[[1]]];
If[!FreeQ[ee,Alternatives@@Drop[vars,k+1]],
Return[Row[{True},Background->ha,Frame->True] &&Row[{Opener[Dynamic[Evaluate[x]]],PaneSelector[{False->Row[{Subscript["b",Sequence@@ind,1]},Frame->True,Background->hb],True->Row[{cview[cad,vars,rul,dvars,Append[ind,1]]},Frame->True]},Dynamic[Evaluate[x]],ImageSize->Automatic]}]]];
Do[ee=cc[[i]];
{a,b}=If[Head[ee]===And, {ee[[1]],Drop[ee,1]},{ee,True}];
na=N[a/.rul];
rr=If[na===True,rul,{}];
If[na=!=False,ans=ans||Row[{Tooltip[a,na]},Background->ha,Frame->True] &&Row[{Opener[Dynamic[Evaluate[x]]],PaneSelector[{False->Row[{Subscript["b",Sequence@@ind,i]},Frame->True,Background->hb],True->Row[{cview[b,vars,rr,dvars,Append[ind,i]]},Frame->True]},Dynamic[Evaluate[x]],ImageSize->Automatic]}]],
{i,Length[cc]}];
ans]

CAFView[cad_,vars_,rul_:{}]:=With[{dvars=Table[Unique["dynamicx"],{Length[vars]}]},cview[cad,vars,rul,dvars,{}]]

Combining cylindrical formulas

The following new function (in Mathematica 8 it is in Reduce` context, since the interface design is not ready) allows to compute set-theoretic operations with semialgebraic sets represented by cylindrical algebraic formulas, that is, to compute Boolean operations on the formulas.

Reduce`CADCombine[bf,cads,vars] gives a cylindrical algebraic formula in vars representing Apply[bf, cads].

cads is a list of cylindrical algebraic formulas in vars and bf is a Boolean pure function with the number of variables equal to the length of cads.

bf may contain quantifiers. In this case cads should be cylindrical algebraic formulas in vars and quantifier variables, ordered starting with the outermost quantifiers.

Combining sets with simple intersections

Compute Boolean combinations of the following sets.

A=(((x+1))^4)+(y^4)+(z^4)<4&&(((x+2))^2)+(y^2)+(z^2)<5;
B=(((x-1))^4)+(y^4)+(z^4)<4&&(((x-2))^2)+(y^2)+(z^2)<5;

Show[{RegionPlot3D[ineq1,{x,-3,3},{y,-2,2},{z,-2,2},PlotStyle->Directive[Red,Opacity[0.5]],Mesh->None,PlotPoints->50],RegionPlot3D[ineq2,{x,-3,3},{y,-2,2},{z,-2,2},PlotStyle->Directive[Blue,Opacity[0.5]],Mesh->None,PlotPoints->50]}]

Direct computation of A⋂B and A⋃B.

(i1=CylindricalDecomposition[A∧B,{x,y,z}]);//Timing

(u1=CylindricalDecomposition[A∨B,{x,y,z}]);//Timing

Computation of A⋂B and A⋃B using CADCombine.

(cadA=CylindricalDecomposition[A,{x,y,z}]);//Timing

(cadB=CylindricalDecomposition[B,{x,y,z}]);//Timing

(i2=Reduce`CADCombine[#1&&#2&,{cadA,cadB},{x,y,z}]);//Timing

(u2=Reduce`CADCombine[#1||#2&,{cadA,cadB},{x,y,z}]);//Timing

i1===i2

Reduce`CADCombine[#1⧦#2&,{u1,u2},{x,y,z}]//Timing

Length[BooleanConvert[#]]&/@{u1,u2}

Independent quantifier elimination subproblems

Problem: Suppose a cubic polynomial function f=a x+b+(x^3) maps [-1,1] into [-1,1]. What is the largest possible distance between its real roots?

Direct approach - single quantifier elimination problem.

Reduce[(Subscript[∃, {a,b}])(Subscript[∀, z,-1≤z≤1])(Subscript[∃, {x,y}])((x^3)+a x+b==0&&(y^3)+a y+b==0&&x-y==d&&-1≤(z^3)+a z+b≤1),d,Reals]//Timing

The problem can be solved faster by solving two independent quantifier elimination subproblems and combining the results in a final quantifier elimination step.

For what values of a and b, f maps [-1,1] into [-1,1]?

(cad1=CylindricalDecomposition[ForAll[z,-1≤z≤1,-1≤z^3+a z+b≤1],{a,b}])//Timing

For what values of a and b, the distance between roots of f can be equal to d?

(cad2=CylindricalDecomposition[Exists[{x,y},x^3+a x+b==0&&y^3+a y+b==0&&x-y==d],{d,a,b}])//Timing

For what values of d there exist a and b satisfying the computed conditions.

Reduce`CADCombine[Exists[{a,b},#1&&#2]&,{cad1,cad2},d]//Timing