Real zeroes of arbitrary real polynomials

Directly to the input form

• Sometimes it is necessary to check whether there are real zeroes of a polynomial or to compute only these. If one knows that there are only real zeroes, then the Newton-Maehly method would be the choice.
• In order to find only the real zeroes if there could be none or in addition also complex ones a different approach must be taken. The following one is rather simple and robust.
• We use underestimating (if the polynomial is positive at some point) or overestimating (in case it is negative) by a parabola. This because we can obtain quadratic convergence locally like in Newton's method, whereas a linear under/over-estimator would yield slow linear convergence only. This under-/overestimating process is done in not too large subintervals in order to be efficient. Since an interval which contains all the zeroes, if any, is easily computed, we might begin at a point x₀ ''on the left'' and provide a stepsize h > 0. We now compute the complete Taylorexpansion of p at x₀. In this program we use the Shaw-Traub algorithm for this: Let
  p(x) = ∑_{i=0 to n}a_ix^n-i and a(0)=1.
  Since this is already the Taylor expansion for x=0 we assume in the following x₀ < > 0. Then we fill the upper triangular part of a matrix of dimension n+2 , in this case numbered with indices ranging from -1 to n as follows:
  t(i,i) = x₀ⁿ , i=0,...,n
  t(-1,i) = a_i+1x₀^n-(i+1), i=0,...,n-1 .
  This makes 2n multiplications. Now the rows 0 to n-1 are filled, including the up to now uninitialized column n in the following way: each element is obtained as the sum of the element left of it and the one above the element left of it:
  t(j,i) = t(j,i-1)+t(j-1,i-1), i=j+1,...,n ; j=0,...,n-1
  This makes n(n+1)/2 additions. Then there holds
  p^(j)(x₀)/(j!) = b(j) = t(j,n)/x₀^j, j=0,..,n.
  which requires n divisions. Then
  p(z) = ∑_{i=0 to n} b_i(z-x₀)ⁱ .
  For not too small n this is obviously more efficient than Horners complete scheme even if the division is much slower than the multiplication. The under-/overestimation is now obtained by modifying the second order term in p:
  For x in [x₀,x₀+h] there holds:
  p(x) is in
  p(x₀)+p'(x₀)*(x-x₀)+(x-x₀)²*( p''(x₀)/2+(-)∑ _{k=3 to n} |p^(k)(x₀)|*h^k-2/(k!))
  
  For p > 0 we take the lower and for p < 0 the upper bound. This is a parabola whose zeroes are computed easily in a stable fashion using Vieta's rule. Should there be no real zero in [x₀,x₀+h] we replace x₀ by x₀+h and if there is a real zero z < x₀+h we replace h by z-x₀. Should there be a real zero z in [x₀,x₀+h] we replace x₀ by z. Since only the second order term is modified this maintains local quadratic convergence. If a real zero is found we deflate p by Horners scheme, move a little bit back and restart the method: this even allows us to identify multiple zeroes with their limiting accuracy.
• This process is continued until all real zeroes are found resp. the predefined search interval is exhausted.

Input

• The degree n of the polynomial.
• You now have the choice either to specify the zeroes of p as n complex numbers (since we take the real part of the coefficients only, be careful: if you forget to specify the zeores either with imaginary part zero or in complex conjugate pairs nonsense will be computed) or to specify the n+1 coefficients a₀, ... ,a_n of p. Important: The leading coefficient is a₀ (it is multiplied with xⁿ) and must be nonzero. You don't need to have it equal to one!
• An interval in which the zeroes should be sought or otherwise the information that all real zeroes must be found. In this case we use Marden's bound
  2 max_{i=1 to n}{ |a_i/a₀|^1/i }
  for defining the search interval (symmetric to zero).

Output

• A list of the zeroes found.
• A graph of the normalized polynomial p(x)/a₀ with the testpoints used. For the smallest zero found also the values of the bounding parabola are shown. Sometimes we use some zooming of the plot in order to improve visibility (parts of the interval with no zeroes but large function values are omitted).

Questions?

• How can you show the second order convergence in case of a simple zero?
• What about the speed of convergence in case of multiple zeroes?
• Why do we use explicit deflation here? Imagine you wanted to work with the original polynomial throughout: which problem must you solve then?

To the input form