MATLAB Examples

Numerical accuracy and function residuals

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

Solve the polynomial

$$p(z) = z^{10} - 30z^9 + 2, $$

which has large function residuals for the default tolerance.

This example illustrates the use of a function in calls to BertiniLab.

f = @(x) x.^10-30*x.^9+2; % The .^ allows vector inputs
polysyms z
poly_system = BertiniLab('function_def',f(z),'variable_group',z);
poly_system = solve(poly_system);
sols = poly_system.match_solutions('raw_solutions');

A polysym object, z, was input to the function to create the Bertini call. Now we can use the same function on the numerical results to see how close to zero it is.

zsols = double(sols.z);
fprintf('%17s %36s\n','z','f(z)')
fprintf('%15.11f + %15.11fi %15.11f + %15.11fi\n',[real(zsols) imag(zsols) real(f(zsols)) imag(f(zsols))].')
                z                                 f(z)
  0.74221900043 +   0.00000000000i   0.00000000000 +   0.00000000000i
  0.56732724477 +   0.47779035001i   0.00000000000 +   0.00000000000i
  0.12660396467 +   0.72957585154i   0.00000000000 +  -0.00000000000i
 -0.37105944022 +   0.63923724056i   0.00000000000 +  -0.00000000000i
 -0.69398126943 +   0.25187220550i   0.00000000000 +   0.00000000000i
 -0.69398126943 +  -0.25187220550i  -0.00000000000 +   0.00000000000i
 -0.37105944022 +  -0.63923724056i   0.00000000000 +  -0.00000000000i
  0.12660396467 +  -0.72957585154i   0.00000000000 +  -0.00000000000i
  0.56732724477 +  -0.47779035001i   0.00000000000 +   0.00000000000i
 30.00000000000 +   0.00000000000i   0.00000000000 +   0.00000000000i

The book reports that f(z) is nowhere near zero for the solution that is near 30. Repeated runs of this example indicate that this seems to be the typical behavior, though the value of f(z) does vary wildly between runs, perhaps due to the random choices made by Bertini during the run.

A greater accuracy may be obtained for the point using a different configuration setting:

poly_system.config = struct('TrackType',-4,'MPType',1,'Precision',100);

The point must be assigned to the property starting_points as class polysym. So that BertiniLab knows which values are assigned to which variables, it is paired with the variable z.

[~,idx] = min(abs(zsols-30));
starting_points = poly_system.unpack_solutions(sols);
poly_system.starting_points = starting_points(:,idx); %[z polysym(zsols)];
poly_system = poly_system.solve;
fvals = poly_system.read_solutions('function');
disp(fvals)
0.31700000000059347371461626607925e-1

We redo the calculations at a higher precision and compare the results:

poly_system.config = struct('FinalTol',1e-19);
poly_system = solve(poly_system);
sols = poly_system.match_solutions('raw_solutions');

zsols = double(sols.z);
[~,idx] = min(abs(zsols-30));

starting_points = poly_system.unpack_solutions(sols);
poly_system.starting_points = starting_points(:,idx); %[z polysym(zsols)];
poly_system.config = struct('TrackType',-4,'MPType',1,'Precision',100);
poly_system = poly_system.solve;

fvals = poly_system.read_solutions('function');
disp([z sols.z(idx)])
disp([polysym('fval') fvals])
z 0.299999999999998983894731614066e2+(-0.954488878703767090457086436084e-23)*I
fval 0.23602719778637037961743772029877e-9+(-0.18787204599525102254812522960387e-9)*I