MATLAB Examples

Performing a Newton iteration

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

The option TrackType = -2 will tell Bertini to perform a single Newton iteraction on a given set of points. An option of -1 will perform the iteration and also approximate the condition number.

This example uses the circle-parabola system again. Again, a start file is needed. Starting points must be provided in a matrix with one column per point.

circle_parabola_intersection.config  = struct('TrackType',-1,'MPType',1,'Precision',64);
circle_parabola_intersection.starting_points = [.75 2; 1e-15 .25]';
circle_parabola_intersection = solve(circle_parabola_intersection);
results = circle_parabola_intersection.solve_summary;
istart = strfind(results,'Evaluating point');
disp(results(istart:end))
Evaluating point 0 of 2

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

successFlag:    Success or failure of Newton iterations.
newPoints:      The point resulting from the Newton iteration.
newtonResidual: The Newton residual for each point.
CN:             Approximation of the condition number at each point.
------------------------------------------------------------------------------------------

All the solutions are classified as real and finite, so we will ignore the imaginary part.

sols = circle_parabola_intersection.match_solutions('newPoints');
condnum = circle_parabola_intersection.read_solutions('CN').';
successFlag = circle_parabola_intersection.read_solutions('successFlag').';
disp([x y polysym('condition number') polysym('flag'); ...
    sols.x sols.y polysym([condnum successFlag])])
                          x                         y      condition number flag
                                                                                
  0.930431071276141698754e0 0.154417413572343149809e1 4.407326555258761e+00    0
                                                                                
0.999999999999999999994e-15 0.250000000000000000000e0                    -1   -1

To make the Newton iteration work, we redo this example with a higher precision:

circle_parabola_intersection.config.Precision  = 128;
circle_parabola_intersection = solve(circle_parabola_intersection);
sols = circle_parabola_intersection.match_solutions('newPoints');
newtonResidual = circle_parabola_intersection.read_solutions('newtonResidual').';
successFlag = circle_parabola_intersection.read_solutions('successFlag').';
disp([x y polysym({'Newton residual','flag'});...
    sols.x sols.y polysym([newtonResidual successFlag])])
                                             x                                              y                              Newton residual flag
                                                                                                                                               
  0.9304310712761416986769099445155783183919e0   0.1544174135723431498079385403329065300897e1 4.902372793100671314134911519201955119926e-1    0
                                                                                                                                               
-0.1713709677419354838709677419310405827262e14 -0.9274193548387096774193548387111342351739e-1 1.713709677419354838709677419753147762744e13    0

The large Newton residual for the second solution indicates that a large step must be taken, which is why it did not work at the lower precision.