Assignment One
Looking at the result of 'i' following [x,r,i] = graeffe(p,1.0e-4,10), we see a nice uniform pattern. The separation of roots works very well. One could say we gain about one decimal place per iteration. There is not much difference in the number of iterations, for small (e.g., 2 or 3) and large degrees (e.g., 7 or 8). The table below summarizes the experiments.
+-----+------------------------+
| | Number of Iterations, |
| d | for three trials |
+-----+------------------------+
| 2 | 3 3 4 |
| 3 | 5 4 5 |
| 4 | 4 4 3 |
| 5 | 4 4 4 |
| 6 | 4 3 4 |
| 7 | 3 4 3 |
| 8 | 3 3 4 |
+-----+------------------------+
Assignment Two
We observe that, although the residuals get very small, the distance to the exact root 1 gets larger as the degree increases. In particular, the error appears to be (1.0e-10)^(1/d), for d = 2,3,..,8. The 1.0e-10 is an estimate for the residual. This error can be explained by the flat shape of the polynomial around 1: residuals can get very small at relatively large distance from the exact root. The triangular table of errors for increasing degrees d = 2,3,..,8 illustrates the point.
Errors with tolerance = 1.0e-08 8.5e-05 8.5e-05 2.1e-03 0.0e+00 2.1e-03 5.4e-03 1.6e-03 1.6e-03 5.4e-03 1.3e-02 5.4e-03 0.0e+00 5.4e-03 1.2e-02 2.8e-02 1.4e-02 4.5e-03 4.5e-03 1.4e-02 2.8e-02 6.3e-02 3.5e-02 1.6e-02 0.0e+00 1.6e-02 3.4e-02 5.9e-02 6.7e-02 4.0e-02 2.2e-02 7.0e-03 6.9e-03 2.1e-02 3.8e-02 6.3e-02 7.1e-02 4.4e-02 2.7e-02 1.3e-02 0.0e+00 1.3e-02 2.6e-02 4.2e-02 6.6e-02 1.5e-01 9.9e-02 6.3e-02 3.6e-02 1.1e-02 1.1e-02 3.4e-02 5.9e-02 9.0e-02 1.3e-01
Note the zero error with odd degrees, which may be explained by b1 -- see formula (6) -- being the sum of the roots. The pattern of approximations around a multiple root is such that errors cancel each other out. The center of a cluster around a multiple root is a much sharper approximation for the root.
To show the relationship with multiplicity, for the error = (1.0e-10)^(1/d), consider the table below:
+-----+-------------------+
| d | (1.0e-10)^(1/d) |
+-----+-------------------+
| 2 | 1.0e-05 |
| 3 | 4.6e-04 |
| 4 | 3.2e-03 |
| 5 | 1.0e-02 |
| 6 | 2.2e-02 |
| 7 | 3.7e-02 |
| 8 | 5.6e-02 |
| 9 | 7.7e-02 |
| 10 | 1.0e-01 |
+-----+-------------------+
When we let the method run with zero tolerance, the errors become zero,
and we get the exact multiple root. When the polynomial has only one
multiple root, division by zero or overflow does not occur and the
method keeps converging to the exact root. The errors are below:
Errors with tolerance = 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00 0.0e+00Assignment Three
The highest degree we could still compute the roots 1:d with two correct decimal places was 12. Looking at the coefficients of p = poly(1:12), we observe huge numbers, in particular, the constant term of p equals 12! = 4.79E+8. This explains the huge residuals and the rapid overflow. The table below summarizes the experiment.
+----+------------+
| d | max|error| |
+----+------------+
| 2 | 1.4552e-11 |
| 3 | 2.1730e-07 |
| 4 | 1.2556e-05 |
| 5 | 1.2376e-04 |
| 6 | 5.4819e-04 |
| 7 | 1.5755e-03 |
| 8 | 3.4865e-03 |
| 9 | 6.5072e-03 |
| 10 | 1.0799e-02 |
| 11 | 1.6465e-02 |
| 12 | 2.3567e-02 |
| 13 | 2.0321e+00 |
+----+------------+
Appendix: scripts for the three assignments.
The data in the tables above is produced by running three scripts:
% Script for Assignment One
nit = zeros(7,3);
for d=2:8
for t=1:3
p = poly(rand(1,d));
[x,r,i] = graeffe(p,1.0e-4,10);
nit(d-1,t) = i;
end;
end;
printf("Number of Iterations :\n");
nit
% Script for Assignment Two
tol = 1.0e-8;
for k=1:2
err = zeros(9,10);
for d=2:10
p = poly(ones(1,d));
[x,r,i] = graeffe(p,tol,100);
for j=1:d
err(d-1,j) = abs(x(j)-1);
end;
end;
printf("Errors with tolerance = %.1e\n",tol);
for d=2:10
for j=1:d
printf(" %.1e", err(d-1,j));
end;
printf("\n");
end;
tol = 0;
end;
printf("Relationship with multiplicity, error = (1.0e-10)^(1/d)\n");
for d=2:10
printf(" %.1e\n", (1.0e-10)^(1/d));
end;
% Script for Assignment Three err = zeros(12,1); for d=2:13 p = poly(1:d); [x,r,i] = graeffe(p,1.0e-8,5); err(d-1,1) = max(sort(x)-[1:d]); end; err(:,1)