Newton-Maehly method: polynomials with real zeroes only

Directly to the input form

• In the following it is always assumed that the polynomial has real and simple zeroes only (for example a polynomial from the set of polynomials orthogonal with respect to some positive weight function on an interval). We also assume the degree n being >= 2 since otherwise the discussion here is senseless. In principle one can use the classical Newton's method for finding the algebraically largest of these zeroes, starting at the right, since one has monotonic and quadratic convergence in this case: since the zeroes of all derivatives of order one to n intersect the zeroes of the derivative of order zero to n-1 respectively right of the largest zero the polynomial and all of its derivatives up to order n have the same sign. Let z₁ be this largest zero and assume w.l.o.g. that p(x₀) > 0. Then in the interval [z₁,x₀] the polynomial p will be strictly convex and hence strictly above its tangent. Clearly the zero of the tangent will be larger than z₁. From
  x_k+1 = x_k - p(x_k)/p'(x_k)
  it follows that the sequence will strictly decrease and x_k+1 > z₁. But since
  p(x_k+1) = p(x_k) + p'(x_k)*(x_k+1-x_k) + (1/2)*p''(xi(k))(x_k+1-x_k)²
  = (1/2)*p''(xi(k))*(p(x_k)/p'(x_k))²
  and, by our assumptions, the denominator here will be bounded from below by some strictly positive value convergence of the function value to zero and hence that of x_k to z₁ follows. Also the quadratic convergence can be seen from this.
• Having found one zero, we can deflate p and since also the zero just found is a legal first approximation of the next one we can simply repeat this as long as the degree is larger than 1.
• Unfortunately this deflation process of zeroes is subject to roundoff effects which will disturb the coefficients of the next polynomial and these perturbations can have disastrous effects on the zeroes themselves. A famous example is Wilkinsons 'perfidious' polynomial with the zeroes 1,2,...,20. In this code you can see this process since we compute for comparison reasons the zeroes this way too.
• It was Maehly's idea to write down Newton's method for the polynomial
  q(x) = p(x)/((x-z₁)*...*(x-z_k))
  with k zeroes of p already found, which after some manipulation results in the simple formula
  x_j+1 = x_j - p(x_j)/(p'(x_j) - p(x_j)*sum_{i=1 to k} 1/(x_j-z_i))
  with the obvious effect that only original values of p enter here, so that the deflation cannot produce the effect mentioned above.
• The situation here is especially satisfying, since violation of strictly monotonic convergence indicates that we are in the region of limiting accuracy and can terminate automatically and safely. Should in this case the polynomial value be clearly above the roundoff level we terminate with an error message. However we also terminate if the polynomial value is below the computed upper bound for its roundoff error using the standard model of floating point arithmetic.
  1.06*machine_eps*sum_{i=0 to n} (2*(n-i)+1)|x|^n-i|a_i|
  
• In order to find a first zero we need a reasonable upper bound for all zeoes. For this we use here Mardens bound:
  2 max{ |a_i/a₀|^1/i } .
  Here a₀ is the highest coefficient. For the zero z_k we cannot use z_k-1 as initial value. Since the zeroes of q' separate those of q we take a Newton step for q' with x₀=z_k-1 and the result of this as the initial value for z_k.

Input

• You specify the degree n first.
• You have now the choice either to specify n real zeroes and to compute the polynomial coefficients from these or to specify the n+1 coefficients a₀, ... ,a_n, in this order. a₀ is the highest coefficient (the one for xⁿ). The first option is useful if you want to experiment with the method regarding its sensitivity w.r.t. the distribution of zeroes. Important: a_n*a₀ < > 0.
• You can require the printed table of all Newton-Maehly steps. This shows the current value of x, p(x), p'(x) and the upper bound for the rounding error in p(x) named ''errbound'' there.

Output

• You get the coefficients of p.
• The list of zeroes as computed by Newton-Maehly, with the correponding polynomial value.
• For comparison reasons the same for Newton's method with explicit deflation.
• A function plot showing p and the sequence of iterates.
• If required, a table of all intermediate steps for Newton-Maehly.
• The possible error messages ''assumptions not satisfied'' mean:
  

  nm1	Newton-Maehly: current x >= zero already found
  nm2	Newton-Maehly: correction leads outside Marden's bound
  nm3	Newton-Maehly: x_new >= x_cur and |p| > 10*errorbound
  nm4	Newton-Maehly: new zero >= zero already found
  nm5	Newton-Maehly: correction for q' leads outside Marden's bound
  nm6	Newton-Maehly: correction for q' >= zero found-1.0e-8*Mardens bound
  d1	Newton deflation: correction leads outside Marden's bound
  d2	Newton deflation: x_new >= x_cur and |p| > 10*errorbound

  All these messages imply that probably there is either a complex conjugate zero or a multiple zero.

Questions?!

• Compare the two variants w.r.t. computing effort and precision!
• Identify situations where Newton-Maehly clearly outperforms deflation.
• What if the polynomial doesn't satisfy the assumptions?
• Multiple zeroes: are they a problem practically?

To the input form