MATLAB Examples

Witness sets and the numerical irreducible decomposition

Example 8.11 from Numerically solving polynomial systems with Bertini, by Daniel J. Bates, Jonathan D. Haunstein, Andrew J. Sommese and Charles W. Wampler (SIAM 2013).

We wish to solve the system of polynomial equations

$$ (y-x^2)(x^2+y^2+z^2-1)(x-2) = 0 $$

$$ (z-x^3)(x^2+y^2+z^2-1)(y-2) = 0 $$

$$ (z-x^3)(y-x^2)(x^2+y^2+z^2-1)(z-2) = 0. $$

The complex solution of this system consists of a surface ($x^2+y^2+z^2=1$), the tangent bundle to a sphere in complex space; three lines ($x=2$ and $z=8$ with $y$ free; $y=2$ and $x=\pm\sqrt{2}$ with $z$ free); one nonlinear curve ($y=x^2$ and $z=x^3$) and an isolated point $(2,2,2)$. To find these solutions (called an irreducible decomposition), Bertini uses witness sets (see the book for details).

config = struct('TrackType',1); % This invokes the decomposition
polysyms x y z
S = x^2+y^2+z^2-1;
T = y-x^2;
U = z-x^3;
poly_system = BertiniLab('function_def',[T*S*(x-2); U*S*(y-2); T*U*S*(z-2)], ...
    'variable_group',[x y z],'config',config);
poly_system = poly_system.solve;

results = poly_system.solve_summary;
istart = strfind(results,'************** Decomposition by Degree');
disp(results(istart:end))
************** Decomposition by Degree **************

Dimension 2: 1 classified component
-----------------------------------------------------
   degree 2: 1 component

Dimension 1: 4 classified components
-----------------------------------------------------
   degree 1: 3 components
   degree 3: 1 component

Dimension 0: 1 classified component
-----------------------------------------------------
   degree 1: 1 component

*****************************************************


The surface has degree 2 and dimension 2; the nonlinear curve has degree 1 and dimension 3.

Components can be sampled by setting TrackType to 2:

poly_system.config = struct('TrackType',2);
poly_system.interactive_choices = {'2','0','100','0','sphere_sample'};
poly_system = poly_system.solve;
points_on_sphere = poly_system.match_solutions('sphere_sample');

results = poly_system.solve_summary;
istart = strfind(results,'Generating samples');
disp(results(istart:end))
Generating samples points on component 0 of dimension 2: 100 points to generate.
Generating 0 of 100
Generating 20 of 100
Generating 40 of 100
Generating 60 of 100
Generating 80 of 100

Writing the sample points to 'sphere_sample'.