Examples for Workshop: Symbolic Computation with Mathematica

Adam Strzeboński

Wolfram Research

Formulation of problems for algebraic solvers

Main classes of problems that are solvable algorithmically

• Systems of linear equations (linear algebra)
• Systems of polynomial equations and inequations in complex variables (Groebner bases)
• Systems of polynomial equations and inequalities in real variables (cylindrical algebraic decomposition)
• Real univariate equations and inequalities involving tame elementary functions (real root isolation)

Mathematica solvers also implement several algorithms for solving certain classes of Diophantine problems (linear, binary quadratic equations, Thue equations, etc.). For other systems Mathematica tries various heuristics in order to reduce the problem to a combination of solving polynomial systems and using known solutions to simple transcendental problems.

Example from StackExchange

f=y E^(x-z);
g=9 x^2+4 y^2+36 z^2-36;
h=x y+y z-1;

Find the stationary points of f over the constraints g==0 ∧h==0.

Direct solution (unsuccessful)

gradf=Grad[f,{x,y,z}];
gradg=Grad[g,{x,y,z}];
gradh=Grad[h,{x,y,z}];

TimeConstrained[sol=Solve[gradf==L gradg+M gradh && g==0 && h==0,{x,y,z,L,M},Reals],300]

Two-step solution (successful)

eq1=Eliminate[gradf==L gradg+M gradh && g==0 && h==0,{L,M}]

sol1=Solve[eq1,{x,y,z},Reals];

y E^(x-z)/.sol1//N

Polynomial formulation

gradf

Since stationary point conditions depend only on directions of gradients and not on their magnitudes, one can divide gradients by nonzero functions.

sol2=Solve[gradf/E^(x-z)==L gradg+M gradh && g==0 && h==0,{x,y,z,L,M},Reals];

y E^(x-z)/.sol2//N

Formulation of problems for exact optimization

Main classes of optimization problems that are solvable algorithmically

• Linear objective functions and constraints (linear programming)
• Polynomial objective functions and constraints (cylindrical algebraic decomposition)
• Real univariate objective functions and constraints involving tame or periodic elementary functions (differentiation/limits/real root isolation)
• Some analytic objective functions and analytic equation constraints in a bounded box (Lagrange multipliers/equation solving)

For other systems Mathematica tries various heuristics in order to reduce the problem to a solvable form.

Example

Find the distance from the point p={3,3,3} to the torus given by the parametrization

p={3,3,3};
par={Cos[a](2+Cos[b]),Sin[a](2+Cos[b]),Sin[b]};

ParametricPlot3D[Evaluate[par],{a,0,2Pi},{b,0,2Pi}]

Direct solution (unsuccessful)

MinValue[{Norm[par-p],0≤a<2Pi&&0≤b<2Pi},{a,b}]

Reduce to analytic objective function in a bounded box

The “analytic objective function in a bounded box” method is not used for the last input, because the use of square root in the objective function adds an inequality constraint. The constraint is always true, but the method is used only if there are no inequalities other than box bounds (checked syntactically).

MinValue[{(par-p).(par-p),0≤a<2Pi&&0≤b<2Pi},{a,b}]//Timing

Specifying the bounds on a and b is necessary.

MinValue[(par-p).(par-p),{a,b}]

Reduce to a polynomial problem

eqs=Join[Thread[{x,y,z}-par],{Sin[a]^2+Cos[a]^2-1,Sin[b]^2+Cos[b]^2-1}];
GroebnerBasis[eqs,{x,y,z},{Sin[a],Cos[a],Sin[b],Cos[b]}]//Timing

poly=%[[2,1]]

MinValue[{Norm[{x,y,z}-p],poly==0},{x,y,z}]//Timing

Region formulations

Trigonometric parametric description with no parameter bounds

The heursitcs successfully reduce the problem, but it takes a long time.

R1=ParametricRegion[par,{a,b}]

RegionDistance[R1,{3,3,3}]//Timing

N[%]

Trigonometric parametric description with parameter bounds

R2=ParametricRegion[{par,0≤a<2Pi&&0≤b<2Pi},{a,b}]

RegionDistance[R2,{3,3,3}]//Timing

N[%]

Polynomial implicit description

R3=ImplicitRegion[poly==0,{x,y,z}]

RegionDistance[R3,{3,3,3}]//Timing

N[%]

Solving for real variables vs. solving over the reals

Solve[(Sqrt[x]) (Sqrt[y])==(Sqrt[z])&&(Sqrt[x-y]) (Sqrt[x+y])==2 (Sqrt[z-1])&&(x^3)+(y^3)+(z^3)==7&&(x|y|z)∈Reals,{x,y,z}]//Timing

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