MATLAB Examples

Contents

Printing a Witness Set

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

We begin by generating an irreducible decomposition of the union of a line and a circle:

polysyms x y
f = (x-y)*(x^2+y^2-1);
poly_system = BertiniLab('variable_group',[x y],'function_def',f, ...
    'config',struct('TrackType',1));
poly_system = poly_system.solve('LineCircle.input');

Bertini is able to store the witness set and linear slicing equations in a form that can be easily reused by Bertini. To do this, set TrackType to 4. The configuration MPType determines the precision to use (0 for double; 1 for fixed higher precision set by Precision; 2 for adaptive precision).

poly_system.config = struct('TrackType',4,'MPType',0);
dim = 1;
component = -2;
file_name = 'witness_points';
linear_system_file = 'linear_system';
poly_system.interactive_choices = {dim,component,file_name,linear_system_file};
poly_system = poly_system.solve('LineCircle.print');

The interactive dialog looks like this:

results = poly_system.solve_summary;
istart = strfind(results,'*************** Components');
disp(results(istart:end))
*************** Components to Print ****************

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


Please select a dimension to print (-1 to quit): 
Dimension 1: 2 classified components
-----------------------------------------------------
   component 0 has degree 1
   component 1 has degree 2


Please select a component to print (-1 to quit, -2 to print all): Enter the name of the file where to write the witness points (max of 255 characters): Enter the name of the file where to write the linear system (max of 255 characters): 
Writing the witness points to 'witness_points' and the linear system to 'linear_system'.


This makes the file witness_points with three points in it and linear_system:

type('linear_system')
variable_group x,y;
function linear1;
constant const1_1,const1_2,const1_3;

const1_1 = 9.467142488389329e-03-2.082095711650707e-01*I;
const1_2 = 2.459452749846511e-01-1.489497140967446e+00*I;
const1_3 = 1.143368969518281e+00+1.007590722015548e+00*I;

linear1 = const1_1 + const1_2*x + const1_3*y;

END;

Computing specific codimensions

The polynomial system

$$ x(x^2-y-z^2) $$

$$ x(x+y+z^2-2)(y-5) $$

$$ x(x^2+x-2)(z-2) $$

has the following irreducible components:

  • The $x=0$ plane
  • Two quadratic curves: $(1,y,z)$, where $y+z^2=1$; and $(-2,y,z)$, where $y+z^2=4$.
  • Two points $(3,5,2)$ and $(-3,5,2)$.

Suppose we are only interested in irreducible components of dimension at least 1. Since the curves are embedded in a three-dimensional space, the maximum codimension desired is 2. The following code finds those components.

config = struct;
config.TrackType = 1;      % compute a numerical irred. decomp.
config.MaxCodimension = 2; % for components of codim <= 2
polysyms x y z
f1 = x*(x^2-y-z^2);
f2 = x*(x+y+z^2-2)*(y-5);
f3 = x*(x^2+x-2)*(z-2);
poly_system = BertiniLab('variable_group',[x y z],'function_def',[f1; f2; f3], ...
    'config',config);

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

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

Dimension 1: 2 classified components
-----------------------------------------------------
   degree 2: 2 components

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


Now we will find just the one-dimensional components (codimension 2).

config = struct;
config.TrackType = 1;           % compute a numerical irred. decomp.
config.SpecificCodimension = 2; % for components of codim = 2
config.WitnessGenType = 1;      % use dimension-by-dimension
poly_system.config = config;

poly_system = solve(poly_system);

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

Dimension 1: 2 classified components
-----------------------------------------------------
   degree 2: 2 components

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