MATLAB Examples

Code used in manuscript

Contents

Section 5.1: Class polysym

x = polysym('x'); y = polysym('y');
disp(x^2+cos(y))
x^2+cos(y)

Equivalence of different polysym definitions:

clear
x = polysym('x'); y = polysym('y'); %#ok<NASGU>
whos
  Name      Size            Bytes  Class      Attributes

  x         1x1               106  polysym              
  y         1x1               106  polysym              

clear
polysyms('x','y')
whos
  Name      Size            Bytes  Class      Attributes

  x         1x1               106  polysym              
  y         1x1               106  polysym              

clear
polysyms x y
whos
  Name      Size            Bytes  Class      Attributes

  x         1x1               106  polysym              
  y         1x1               106  polysym              

Definition of polysym vector:

v = polysym('v',[1 2]); disp(v)
v1 v2
M = polysym('M',2); disp(M)
M1_1 M1_2
         
M2_1 M2_2

Matrix operations have the same syntax as for numeric calculations:

disp(M*v.')
M1_1*v1+M1_2*v2
               
M2_1*v1+M2_2*v2

Numbers can be entered directly or combined with a polysym object:

x = polysym(1+1i);
disp(x*(1-1i))
(1+1*I)*(1-1*I)

polysym variables can be input to a MATLAB function:

f = @(y,z) z^3 + y^2*z + 1;
polysyms y z
disp(f(cos(y),z))
z^3+cos(y)^2*z+1

Importance of the order of evaluation:

polysyms x
disp(9/16*x)
disp(x*9/16)
0.5625*x
x*9/16

Error results from an attempt to assign a polysym object to a component of a numeric array:

y = polysym('y',2);
x = zeros(2);
try
    x(1) = y(1); %#ok<NASGU>
catch error
    disp(error.message)
end
In an assignment  A(I) = B, the number of elements in B and I must be the same.

Ways of initializing an array so it is compatible with polysym:

x = polysym(0,size(y));
disp(x)
0 0
   
0 0
x = y*0;
disp(x)
0 0
   
0 0

Section 5.2: Class BertiniLab

% Solve the system $x^2-1=0$, $y^2-4=0$:
poly_system = BertiniLab('variable_group',{'x','y'}, ...
                'function_def',{'x^2-1';'y^2-4'});
poly_system = poly_system.solve;

Read solutions and match to variables

sols = poly_system.match_solutions('real_finite_solutions');

Convert output to double precision

x = real(double(sols.x));
y = real(double(sols.y));
disp([x(:) y(:)])
    1.0000    2.0000
    1.0000   -2.0000
   -1.0000    2.0000
   -1.0000   -2.0000

Solve $(29/16)z_1^2 -2z_1z_2 = 0$, $z_2-z_1^2=0$.

polysyms z1 z2
f = @(x,y) [(29/16)*x^3 - 2*x*y; y - x^2];
poly_system = BertiniLab('function_def',f(z1,z2), ...
'variable_group',[z1 z2]);
poly_system = poly_system.solve;
disp(poly_system.solve_summary)
    Bertini(TM) v1.4
   (October 23, 2013)

 D.J. Bates, J.D. Hauenstein,
 A.J. Sommese, C.W. Wampler

(using GMP v5.1.3, MPFR v3.1.2)



NOTE: You have requested to use adaptive path tracking.  Please make sure that you have
setup the following tolerances appropriately:
CoeffBound: 3.905043000000e+00, DegreeBound: 3.000000000000e+00
AMPSafetyDigits1: 1, AMPSafetyDigits2: 1, AMPMaxPrec: 1024

Tracking path 0 of 6

Finite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	0		|		0		|   0
Singular	|	1		|		0		|   1
------------------------------------------------------------------------------------------
Total		|	1		|		0		|   1


Finite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	3	|	1		|	0


------------------------------------------------------------------------------------------
The following files may be of interest to you:

main_data:             A human-readable version of the solutions - main output file.
raw_solutions:         A list of the solutions with the corresponding path numbers.
raw_data:              Similar to the previous, but with the points in Bertini's homogeneous
                         coordinates along with more information about the solutions.
real_finite_solutions: A list of all real finite solutions.
finite_solutions:      A list of all finite solutions.
nonsingular_solutions: A list of all nonsingular solutions.
singular_solutions:    A list of all singular solutions.
------------------------------------------------------------------------------------------

Paths Tracked: 6
 Truncated infinite paths: 3 - try adjusting SecurityMaxNorm or set SecurityLevel to 1 in the input file
   Please see 'failed_paths' for more information about these paths.


Solve again with security level reset:

poly_system.config = struct('SecurityLevel',1);
poly_system = poly_system.solve;
disp(poly_system.solve_summary)
    Bertini(TM) v1.4
   (October 23, 2013)

 D.J. Bates, J.D. Hauenstein,
 A.J. Sommese, C.W. Wampler

(using GMP v5.1.3, MPFR v3.1.2)



NOTE: You have requested to use adaptive path tracking.  Please make sure that you have
setup the following tolerances appropriately:
CoeffBound: 6.725360000000e+00, DegreeBound: 3.000000000000e+00
AMPSafetyDigits1: 1, AMPSafetyDigits2: 1, AMPMaxPrec: 1024

Tracking path 0 of 6

Finite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	0		|		0		|   0
Singular	|	1		|		0		|   1
------------------------------------------------------------------------------------------
Total		|	1		|		0		|   1


Finite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	3	|	1		|	0


Infinite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	0		|		0		|   0
Singular	|	0		|		1		|   1
------------------------------------------------------------------------------------------
Total		|	0		|		1		|   1


Infinite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	3	|	0		|	1


------------------------------------------------------------------------------------------
The following files may be of interest to you:

main_data:             A human-readable version of the solutions - main output file.
raw_solutions:         A list of the solutions with the corresponding path numbers.
raw_data:              Similar to the previous, but with the points in Bertini's homogeneous
                         coordinates along with more information about the solutions.
real_finite_solutions: A list of all real finite solutions.
finite_solutions:      A list of all finite solutions.
nonsingular_solutions: A list of all nonsingular solutions.
singular_solutions:    A list of all singular solutions.
------------------------------------------------------------------------------------------

Paths Tracked: 6

Solve using multihomogenous homotopy (needs 5 paths instead of 6):

polysyms z1 z2
f = @(x,y) [(29/16)*x^3 - 2*x*y; y - x^2];
poly_system = BertiniLab('function_def',f(z1,z2), ...
'variable_group',{z1,z2},'config',struct('SecurityLevel',1));
poly_system = poly_system.solve;
disp(poly_system.solve_summary)
    Bertini(TM) v1.4
   (October 23, 2013)

 D.J. Bates, J.D. Hauenstein,
 A.J. Sommese, C.W. Wampler

(using GMP v5.1.3, MPFR v3.1.2)



NOTE: You have requested to use adaptive path tracking.  Please make sure that you have
setup the following tolerances appropriately:
CoeffBound: 3.856657000000e+00, DegreeBound: 4.000000000000e+00
AMPSafetyDigits1: 1, AMPSafetyDigits2: 1, AMPMaxPrec: 1024

Tracking path 0 of 5

Finite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	0		|		0		|   0
Singular	|	1		|		0		|   1
------------------------------------------------------------------------------------------
Total		|	1		|		0		|   1


Finite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	3	|	1		|	0


Infinite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	0		|		0		|   0
Singular	|	0		|		1		|   1
------------------------------------------------------------------------------------------
Total		|	0		|		1		|   1


Infinite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	2	|	0		|	1


------------------------------------------------------------------------------------------
The following files may be of interest to you:

main_data:             A human-readable version of the solutions - main output file.
raw_solutions:         A list of the solutions with the corresponding path numbers.
raw_data:              Similar to the previous, but with the points in Bertini's homogeneous
                         coordinates along with more information about the solutions.
real_finite_solutions: A list of all real finite solutions.
finite_solutions:      A list of all finite solutions.
nonsingular_solutions: A list of all nonsingular solutions.
singular_solutions:    A list of all singular solutions.
------------------------------------------------------------------------------------------

Paths Tracked: 5

Section 6: Customizing BertiniLab

% A function:
type('findRealRoots')
function [sols,varargout] = findRealRoots(polyFun)
%findRealRoots Find the real, finite, isolated roots of a polynomial
%system.
%
% sols = findRealRoots(polyFun) inputs a set of equations (a cell array of
% strings, polysym or sym object) and returns a structure with the
% solutions for each variable in a field with the name of that variable.
% The variable names are automatically obtained from polyFun.
%
% [sols,poly_system] = ... also returns a BertiniLab object, POLY_SYSTEM,
% with the definition of the problem.
%
% Example:
%   sols = findRealRoots({'x^2-1'; 'y^2-1'});
%
% See also: polysolve, BertiniLab

% Find variables automatically
polyFun = polysym(polyFun);
vars = symvar(polyFun);

% Solve and retrieve finite solutions.
poly_system = BertiniLab('function_def',polyFun,'variable_group',vars);
poly_system = solve(poly_system);
sols = poly_system.match_solutions('real_finite_solutions');

% Convert solutions to double precision
sols = structfun(@double,sols,'UniformOutput',false);
sols = structfun(@real,sols,'UniformOutput',false);

% Warn user if there is a problem.
if poly_system.has_failed_paths || poly_system.has_path_crossings
    results = poly_system.solve_summary;
    istart = strfind(results,'Paths Tracked');
    warning(results(istart:end));
end

% Return BertiniLab object, if desired.
if nargout > 1
    varargout{1} = poly_system;
end

A call to the function:

sols = findRealRoots({'x^2-1'; 'y^2-4'});

Section 7: Subclassing and parameter homotopy

anisotropyFcn = @(x,K) 2*[x(1)*(x(2)^2+x(3)^2); x(2)*(x(1)^2+x(3)^2); x(3)*(x(1)^2+x(2)^2) - K*x(3)];
constrFcn = @(x) deal(sum(x.^2)-1,2*x(:));
%We solve this system for a trial value of K:
Kval = 0.2; N = 3;
gradFcn = @(x) anisotropyFcn(x,Kval);
poly_system = solveWithConstraints(gradFcn,N,constrFcn,'config',struct('SecurityLevel',1));
poly_system = solve(poly_system);
disp(poly_system.solve_summary)
    Bertini(TM) v1.4
   (October 23, 2013)

 D.J. Bates, J.D. Hauenstein,
 A.J. Sommese, C.W. Wampler

(using GMP v5.1.3, MPFR v3.1.2)



NOTE: You have requested to use adaptive path tracking.  Please make sure that you have
setup the following tolerances appropriately:
CoeffBound: 3.406454000000e+00, DegreeBound: 4.000000000000e+00
AMPSafetyDigits1: 1, AMPSafetyDigits2: 1, AMPMaxPrec: 1024

Tracking path 0 of 54
Tracking path 20 of 54
Tracking path 40 of 54

Finite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	26		|		0		|   26
Singular	|	0		|		0		|   0
------------------------------------------------------------------------------------------
Total		|	26		|		0		|   26


Finite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	1	|	26		|	0


Infinite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	0		|		0		|   0
Singular	|	0		|		14		|   14
------------------------------------------------------------------------------------------
Total		|	0		|		14		|   14


Infinite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	2	|	0		|	14


------------------------------------------------------------------------------------------
The following files may be of interest to you:

main_data:             A human-readable version of the solutions - main output file.
raw_solutions:         A list of the solutions with the corresponding path numbers.
raw_data:              Similar to the previous, but with the points in Bertini's homogeneous
                         coordinates along with more information about the solutions.
real_finite_solutions: A list of all real finite solutions.
finite_solutions:      A list of all finite solutions.
nonsingular_solutions: A list of all nonsingular solutions.
singular_solutions:    A list of all singular solutions.
------------------------------------------------------------------------------------------

Paths Tracked: 54

%Now we do the ab initio solve for parameter homotopy:
polysyms K
gradFcn = @(x) anisotropyFcn(x,K);
poly_system = solveWithConstraints(gradFcn,N,constrFcn,'config', ...
    struct('ParameterHomotopy',1,'SecurityLevel',1),'parameter',K);
% poly_system = solve(poly_system,'paramotopy.input');
ia = [];
% Sometimes this method obtains the wrong number of starting points.
while (numel(ia) ~= 5)
    poly_system = poly_system.solve('paramotopy.input');
    vars = poly_system.read_solutions('nonsingular_solutions');
    v = round(double(vars)*1e10);
    v = [sort(v(1:2,:)); v(3:4,:)];
    [~,ia] = unique(abs(v.'),'rows');
end
vars_reduced = vars(:,ia(:));

An example of a run using parameter homotopy follows:

poly_system.starting_points = vars_reduced;
poly_system.config.ParameterHomotopy = 2;
poly_system.final_parameters = 0.2;
poly_system = solve(poly_system,'paramotopy.input');
disp(poly_system.solve_summary)
    Bertini(TM) v1.4
   (October 23, 2013)

 D.J. Bates, J.D. Hauenstein,
 A.J. Sommese, C.W. Wampler

(using GMP v5.1.3, MPFR v3.1.2)


NOTE: Please make sure that you have also setup the proper start points file for your parameter homotopy
      with the starting parameters in 'start_parameters' and the target parameters in 'final_parameters'.



NOTE: You have requested to use adaptive path tracking.  Please make sure that you have
setup the following tolerances appropriately:
CoeffBound: 1.901255000000e+00, DegreeBound: 4.000000000000e+00
AMPSafetyDigits1: 1, AMPSafetyDigits2: 1, AMPMaxPrec: 1024

Tracking path 0 of 5

Finite Solution Summary

NOTE: nonsingular vs singular is based on condition number and identical endpoints

		| Number of real solns	|  Number of non-real solns	|  Total
------------------------------------------------------------------------------------------
Non-singular	|	5		|		0		|   5
Singular	|	0		|		0		|   0
------------------------------------------------------------------------------------------
Total		|	5		|		0		|   5


Finite Multiplicity Summary

  Multiplicity	|  Number of real solns	|  Number of non-real solns
------------------------------------------------------------------------------------------
	1	|	5		|	0


------------------------------------------------------------------------------------------
The following files may be of interest to you:

main_data:             A human-readable version of the solutions - main output file.
raw_solutions:         A list of the solutions with the corresponding path numbers.
raw_data:              Similar to the previous, but with the points in Bertini's homogeneous
                         coordinates along with more information about the solutions.
real_finite_solutions: A list of all real finite solutions.
finite_solutions:      A list of all finite solutions.
nonsingular_solutions: A list of all nonsingular solutions.
singular_solutions:    A list of all singular solutions.
------------------------------------------------------------------------------------------

Paths Tracked: 5

Create the solution curves.

N = 200;
Kval = linspace(-1,1,N);
mag = zeros(3,5,N);
lam = zeros(1,5,N);
problemCases = struct.empty;
nCases = 0;
for ii=1:N
    poly_system.final_parameters = Kval(ii);

    poly_system = solve(poly_system,'paramotopy.input');
    sols = poly_system.read_solutions('raw_solutions');

    if poly_system.has_failed_paths || poly_system.has_path_crossings
        nCases = nCases+1;
        problemCases(nCases).number = ii;
        problemCases(nCases).diagnosis = poly_system.solve_summary;
    end

    mag(:,:,ii) = sols(1:3,:);
    lam(:,:,ii) = sols(4,:);
end
if ~isempty(problemCases)
    warning('There were some failed paths or path crossings. See problemCases for more information.')
end

Plot magnetization

map = pmkmp(N,'IsoL');
%map = gray(N);
for ii=1:N
    I = all(abs(imag(mag(:,:,ii)))<1e-12);
    mx = real(mag(1,I,ii)); my = real(mag(2,I,ii)); mz = real(mag(3,I,ii));
    m1 = [mx -mx mx -mx mx -mx mx -mx my -my my -my my -my my -my];
    m2 = [my my -my -my my my -my -my mx mx -mx -mx mx mx -mx -mx];
    m3 = [mz mz mz mz -mz -mz -mz -mz mz mz mz mz -mz -mz -mz -mz];
    line(m1,m2,m3,'Marker','.','LineStyle','none','Color',map(ii,:))
end
colorbar
caxis([min(Kval) max(Kval)])
colormap(map)
view(148,26)
axis off equal

% Add axes and an equator
line([0 0],[0 0],[0 1],'Color','k')
line([0 0],[0 1],[0 0],'Color','k')
line([0 1],[0 0],[0 0],'Color','k')
phi = linspace(0,2*pi,200);
x = cos(phi); y = sin(phi); z = zeros(size(phi));
line(x,y,z,'Color','k')

% Make the colorbar a bit smaller.
h_bar = findobj(gcf,'Tag','Colorbar');
initpos = get(h_bar,'Position');
set(h_bar, ...
   'Position',[initpos(1)+initpos(3)*0.15 initpos(2)+initpos(4)*0.15 ...
      initpos(3)*0.7 initpos(4)*0.7])
%In real space there are 16 distinct curves. How many are there in complex
%space? We can find out by running a positive-dimensional solve for the
%system. As written, |solveWithConstraints| does not handle this, so we
%will call BertiniLab directly:
m = sym('m',[3 1]);
syms K lambda
E = K*(1-m(3)^2) + m(1)^2*m(2)^2+m(2)^2*m(3)^2+m(3)^2*m(1)^2  ...
+ lambda*(m(1)^2+m(2)^2+m(3)^2-1);
f = [diff(E,m(1)); diff(E,m(2)); diff(E,m(3)); diff(E,lambda)];
poly_system = BertiniLab('variable_group',[m; lambda; K], ...
'function_def',f,'config',struct('TrackType',1));
poly_system = solve(poly_system);
disp(poly_system.solve_summary)
    Bertini(TM) v1.4
   (October 23, 2013)

 D.J. Bates, J.D. Hauenstein,
 A.J. Sommese, C.W. Wampler

(using GMP v5.1.3, MPFR v3.1.2)



NOTE: You have requested to use adaptive path tracking.  Please make sure that you have
setup the following tolerances appropriately:
CoeffBound: 5.535781000000e+00, DegreeBound: 3.000000000000e+00
AMPSafetyDigits1: 1, AMPSafetyDigits2: 1, AMPMaxPrec: 1024


Tracking regeneration codim 1 of 4: 3 paths to track.
Tracking path 0 of 3

Sorting codimension 1 of 4: 3 paths to sort.
Sorting 0 of 3

Preparing regeneration codim 2 of 4: 6 witness points to move.
Moving 0 of 6

Tracking regeneration codim 2 of 4: 9 paths to track.
Tracking path 0 of 9

Sorting codimension 2 of 4: 9 paths to sort.
Sorting 0 of 9

Preparing regeneration codim 3 of 4: 18 witness points to move.
Moving 0 of 18

Tracking regeneration codim 3 of 4: 27 paths to track.
Tracking path 0 of 27
Tracking path 20 of 27

Sorting codimension 3 of 4: 27 paths to sort.
Sorting 0 of 27
Sorting 20 of 27

Preparing regeneration codim 4 of 4: 24 witness points to move.
Moving 0 of 24
Moving 20 of 24

Tracking regeneration codim 4 of 4: 48 paths to track.
Tracking path 0 of 48
Tracking path 20 of 48
Tracking path 40 of 48

Sorting codimension 4 of 4: 48 paths to sort.
Sorting 0 of 48
Sorting 20 of 48
Sorting 40 of 48


************ Regenerative Cascade Summary ************

NOTE: nonsingular vs singular is based on rank deficiency and identical endpoints

|codim|   paths   |witness superset| nonsingular | singular |nonsolutions| inf endpoints | other bad endpoints
----------------------------------------------------------------------------------------------------------------
| 1   |   3       |   0            |  0          |  0       |  3         |   0           |  0
| 2   |   9       |   0            |  0          |  0       |  9         |   0           |  0
| 3   |   27      |   0            |  0          |  0       |  24        |   3           |  0
| 4   |   48      |   26           |  26         |  0       |  0         |   22          |  0
----------------------------------------------------------------------------------------------------------------
|total|   87

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



*************** Witness Set Summary ****************

NOTE: nonsingular vs singular is based on rank deficiency and identical endpoints

|codim| witness points | nonsingular | singular 
-------------------------------------------------
| 4   |   26           |  26         |  0       
-------------------------------------------------

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


Calculating traces for codimension 4.
Calculating 0 of 26
Calculating 20 of 26

Using monodromy to decompose codimension 4.
Performing monodromy loops: 16 points left to classify
Performing monodromy loops: 16 points left to classify
Performing monodromy loops: 16 points left to classify
Performing monodromy loops: 16 points left to classify
Performing monodromy loops: 16 points left to classify
Performing monodromy loops: 16 points left to classify
Performing monodromy loops: 16 points left to classify

Using combinatorial trace test to decompose codimension 4.


************* Witness Set Decomposition *************

| dimension | components | classified | unclassified
-----------------------------------------------------
|   1       |   14       |   26       |  0
-----------------------------------------------------

************** Decomposition by Degree **************

Dimension 1: 14 classified components
-----------------------------------------------------
   degree 1: 10 components
   degree 4: 4 components

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


BertiniClean % Delete all the files created by Bertini