Below is an outline for a possible solution to the project. The computations were done on a mac, Powerbook G4, with Octave.
1. Solution to Assignment One
We apply run_mueller three times on polynomials with random roots of degrees 3 to 10 and record the number of iterations to reach the first root in the matrix A. The (i,j)-th entry of A lists the number of iterations to reach the first root of the i-th polynomial in the j-th trial.
d trials --------------- 3 9 8 8 4 12 13 15 5 15 15 9 6 18 12 10 7 17 16 12 8 15 17 17 9 21 25 15 10 24 38 26The averages over the three columns are
16.375 18.000 14.000We see that the global convergence is fairly stable.
To examine the local convergence we generate a polynomial with known roots and start mueller with x0, x1, and x2 at a distance of 1.0e-2 from the first root. Also here we conduct 3 trials for degrees ranging from 3 to 10.
d trials ------------- 3 3 3 3 4 3 3 5 5 3 3 3 6 3 4 2 7 3 3 4 8 3 3 3 9 3 4 4 10 3 1 2The average over the three columns are
3.0000 3.0000 3.2500Compared to the global convergence, the number of steps is more uniform, and of course faster.
Starting at 1.0e-2 accuracy of the root, it takes about three steps to reach a residual less than 1.0-12. This corresponds to a superlinear or almost quadratic convergence.
2. Solution to Assignment Two
We take a polynomial with 5 roots:
1.3079 2.6067 3.2021 1.7204 3.9319In both runs, the accuracy of run_mueller and ref_mueller was 1.0e-8 to see the effect of the accumulation of errors on the roots better.
The reported residuals of run_mueller are
3.8294e-11 4.0297e-10 3.4722e-11 1.2826e-14 0.0000e+00But the actual residuals are much higher:
3.8294e-11 2.0308e-10 4.8764e-10 2.9000e-10 7.0822e-10We observe that the actual residuals are higher for the roots which are computed last.
And here is the difference between the roots
2.4950e-11
-2.7882e-12
-1.0685e-10
1.1276e-10
-2.8072e-11
Since the roots are sorted, we no longer distinguish between roots
computed first and last.
ref_mueller reports the residuals
1.9895e-13 6.8378e-21 1.0042e-22 1.5632e-13 0.0000e+00wich correspond to the actual residuals
1.9895e-13 6.8378e-21 1.0042e-22 1.5632e-13 8.5265e-13and all roots are computed with greater accuracy
6.6613e-15
-1.6209e-14
9.9476e-14
4.4409e-15
-9.5035e-14
3. Solution to Assignment Three
If we apply ref_mueller on (x-1)^d, for d = 3:10, we see that the distance between the computed roots to the exact root is about the same for all roots. The distance increases as d grows, and is bounded by 10^(-1/m). A plot of the computed roots shows they are clustered around the exact root, on a circle with radius of size about 10^(-1/m). Increasing the precision does not help much and only leads to division by zero.
The application of a modified Newton's method using the known multiplicity leads to more accurate results.
The data are below:
A = Columns 1 through 6: 4.9666e-05 3.6954e-05 5.4574e-05 0.0000e+00 0.0000e+00 0.0000e+00 9.4057e-04 9.0701e-04 9.2884e-04 8.9160e-04 0.0000e+00 0.0000e+00 3.4605e-03 3.4112e-03 3.4048e-03 3.3935e-03 3.4317e-03 0.0000e+00 9.0691e-03 9.0431e-03 9.0354e-03 9.0636e-03 8.9250e-03 8.9951e-03 1.7770e-02 1.7769e-02 1.7770e-02 1.7710e-02 1.7695e-02 1.7859e-02 2.9953e-02 2.9837e-02 3.0026e-02 3.0070e-02 2.9828e-02 2.9881e-02 4.4007e-02 4.3798e-02 4.3939e-02 4.3666e-02 4.3615e-02 4.3641e-02 6.0912e-02 6.0991e-02 6.1068e-02 6.0941e-02 6.0919e-02 6.1207e-02 Columns 7 through 10: 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.7790e-02 0.0000e+00 0.0000e+00 0.0000e+00 2.9911e-02 3.0038e-02 0.0000e+00 0.0000e+00 4.3768e-02 4.3876e-02 4.3955e-02 0.0000e+00 6.1093e-02 6.1147e-02 6.1091e-02 6.0913e-02Calculation of 1/log10(A(i,1))
ans = -2.3235e-01 ans = -3.3040e-01 ans = -4.0636e-01 ans = -4.8961e-01 ans = -5.7133e-01 ans = -6.5636e-01 ans = -7.3720e-01 ans = -8.2284e-01After application of a modified Newton's method:
C = Columns 1 through 6: 1.1548e-08 1.5135e-08 3.0269e-08 0.0000e+00 0.0000e+00 0.0000e+00 6.1573e-08 2.2204e-16 1.1102e-16 1.0000e+00 0.0000e+00 0.0000e+00 1.4780e-06 2.0015e-16 4.4604e-16 4.1012e-04 2.0005e+00 0.0000e+00 1.6638e-05 3.1402e-16 1.0116e-12 7.8129e-14 5.5511e-17 3.0000e+00 7.4504e-05 1.2676e-14 8.9509e-16 9.9301e-16 1.5073e-14 4.0972e-13 2.5161e-04 5.7524e-11 6.2696e-15 9.9762e-12 4.6875e-14 2.3379e-08 1.1814e-04 4.9921e-14 8.8630e-15 1.6514e-12 8.0398e-14 5.3448e-05 1.0401e-03 5.0952e-16 4.1581e-09 6.3101e-11 1.3771e-08 1.8706e-06 Columns 7 through 10: 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 4.0000e+00 0.0000e+00 0.0000e+00 0.0000e+00 1.2524e-13 5.0000e+00 0.0000e+00 0.0000e+00 6.2563e-11 1.3527e-14 6.0000e+00 0.0000e+00 3.2236e-07 2.5358e-10 2.6220e-12 7.0000e+00
4. Solution to Assignment Four
We run ref_mueller on p = poly(1:d), for d from 3 to 22. The maximal error at each degree d is
d max(error) --------------- 3 2.2204e-16 4 7.1055e-15 5 2.6645e-14 6 1.2790e-13 7 1.1680e-12 8 5.4614e-12 9 4.0462e-11 10 7.2676e-10 11 4.4514e-09 12 5.5132e-09 13 3.0954e-08 14 3.4696e-07 15 1.1700e-06 16 8.4886e-06 17 8.0682e-05 18 7.1884e-04 19 1.9306e-03 20 1.0748e-02 21 1.0282e-01 22 6.6215e-01We see that the maximal error of the last polynomial is already too large, the roots are, sorted in increasing order
1.00000 + 0.00000i 2.00000 + 0.00000i 3.00000 - 0.00000i 4.00000 + 0.00000i 5.00000 - 0.00000i 6.00000 + 0.00000i 6.99995 - 0.00000i 8.00023 + 0.00000i 8.99937 - 0.00001i 10.00036 + 0.00002i 11.00646 - 0.00006i 11.96470 + 0.00002i 13.14654 - 0.00006i 13.73587 + 0.00070i 15.45578 - 0.38765i 15.46511 + 0.39030i 17.54919 - 0.33963i 17.56899 + 0.33733i 19.14212 + 0.00316i 19.96595 - 0.00708i 20.99796 + 0.00367i 22.00143 - 0.00072i
As the maximal error is already larger than .6, we can no longer distinguish between 15 and 16 and between 17 and 18. So d = 21 is the largest degree possible.
This polynomial is so hard to solve because the coefficients vary so widely in magnitude. The constant term is d!, which is 2.4e+18, for d = 20. A tiny change in one of the coefficients, caused by a representation error introduces imaginary parts in the roots. The roots of this polynomial are ill-conditioned as small changes in the coefficients cause large changes in the roots.