Parameter homotopy: Sextic Ab Initio
Example 6.1 from Numerically solving polynomial systems with Bertini, by Daniel J. Bates, Jonathan D. Haunstein, Andrew J. Sommese and Charles W. Wampler (SIAM 2013).
As a simple example of using a parameter homotopy, we solve a sextic polynomial in one variable. The first step is to do an ab initio run to get the starting parameters.
This example departs a little from the one in the book to illustrate how MATLAB handles polynomial coefficients. In the book, the parameters are
polysyms a0 a1 a2 a3 a4 a5 a6 x % and the equation is f = a0+x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*a6))))); disp(f)
a0+x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*a6)))))
This can be set up in vector form using the overloaded method POLYVAL, but the conventions are different in MATLAB. Indices start at 1 and the lowest coefficient is multiplied by the highest power:
p = polysym('a',[1 7]);
f = polyval(p,x);
disp(f)
x*(x*(x*(x*(x*(x*a1+a2)+a3)+a4)+a5)+a6)+a7
We proceed with the ab initio run.
config = struct('ParameterHomotopy',1); poly_system = BertiniLab('function_def',f,'variable_group',x, ... 'parameter',p,'config',config); poly_system = solve(poly_system,'sextic.input');
Now for the parameter homotopy run:
poly_system.config = struct('ParameterHomotopy',2); % Make nonsingular solutions the starting points poly_system.starting_points = poly_system.read_solutions('nonsingular_solutions'); % Define final parameter values pvals = [1.1 2.4 0.8 3.6 -0.52 -1.8 4.4]; poly_system.final_parameters = pvals; poly_system = solve(poly_system,'sextic.input'); % Display summary of results summary = poly_system.solve_summary; iStart = strfind(summary,'NOTE: nonsingular'); disp(summary(iStart:end))
NOTE: nonsingular vs singular is based on condition number and identical endpoints | Number of real solns | Number of non-real solns | Total ------------------------------------------------------------------------------------------ Non-singular | 2 | 4 | 6 Singular | 0 | 0 | 0 ------------------------------------------------------------------------------------------ Total | 2 | 4 | 6 Finite Multiplicity Summary Multiplicity | Number of real solns | Number of non-real solns ------------------------------------------------------------------------------------------ 1 | 2 | 4 ------------------------------------------------------------------------------------------ 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
The solutions can be compared with the output of the MATLAB function ROOTS, which solves single-variable polynomials.
sols = poly_system.read_solutions('finite_solutions');
r = roots(pvals);
disp(num2str(sort([double(sols).' r])))
0.67389-0.56532i 0.67389-0.56532i 0.67389+0.56532i 0.67389+0.56532i -1.114+4.1633e-17i -1.114+0i -0.015442-1.3949i -0.015442-1.3949i -0.015442+1.3949i -0.015442+1.3949i -2.3847-1.1102e-16i -2.3847+0i
It would be easy to automate this process for a large number of parameter values, in one or more loops. The package Paramotopy is designed to do this efficiently. We hope to extend the BertiniLab interface to Paramotopy in a future release.