MATLAB Examples

Root sharpening after a homotopy

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

Our goal is to find the intersection of a hyperbola and parabola:

$$ x^2+y^2=1, \quad x=y^2+1. $$

The solutions are two nonsingular points $\left(0,\pm\sqrt{-1}\right)$ and a point $(1,0)$ of multiplicity 2.

polysyms x y
hyperbola = x^2 - y^2 - 1;
parabola = x - y^2 - 1;

poly_system = BertiniLab('function_def',[hyperbola; parabola],'variable_group',[x y]);
poly_system = solve(poly_system);
sols = poly_system.match_solutions('nonsingular_solutions');

The solutions approximating $(0,\pm i)$ are shown below:

xsols = double(sols.x);
ysols = double(sols.y);
fprintf('%17s %35s\n','x','y')
fprintf('%15.11f + %15.11fi  %15.11f + %15.11fi\n', ...
    [real(xsols) imag(xsols) real(ysols) imag(ysols)].')
                x                                   y
 -0.00000000000 +  -0.00000000000i    0.00000000000 +  -1.00000000000i
  0.00000000000 +   0.00000000000i    0.00000000000 +   1.00000000000i

Now sharpen to 20 digits (which we expect to fail for the singular solution):

poly_system.config = struct('SharpenOnly',1,'SharpenDigits',20);
poly_system.interactive_choices = {'1'};
poly_system = poly_system.solve;

For comparison, here are the revised values of x:

sols = poly_system.match_solutions('nonsingular_solutions');
disp([x y; sols.x sols.y])
                                                                             x                                                                            y
0.383098700203159704858806206684e-78+(-0.450700145517885898187670732226e-78)*I 0.631088724176809444329382852226e-29+(-0.100000000000000000000000000000e1)*I
-0.237097864677219098477898889898e-77+(0.311207928631001598647304138103e-78)*I   -0.631088724176809444329382852226e-29+0.100000000000000000000000000003e1*I

The result is much closer to the correct value of zero.

To sharpen the singular solution, the endgame must be chosen instead of Newton's method:

poly_system.interactive_choices = {'6','2','1'};
poly_system = poly_system.solve;

% Display the summary of the run
result = poly_system.solve_summary;
istart = strfind(result,'Endpoints Considered');
Endpoints Considered for Sharpening: 4
 Finite Sharpening Successes: 4