Remark: Polynomial Rootfinding is not a stable problem.
 

• Ex: (x-2)² =0 has solution x=2. But (x-2)² = 0.000001 has solutions x=2.001 and x=1.997 in exact arithmetic. So, a change in the data of 0.000001 leads to a change in the final answer of 0.001, a relatively large change.
• Ex. (Wilkinson) f(x)=(x-20)(x-19)(x-18)...(x-2)(x-1) = x²⁰ + 210x¹⁹ + ... obviously has roots 1,2,3,...,19,20. However, it can be shown that if the coefficient 210 is changes by a small change on the order of machine roundoff that the resulting function has in exact arithmetic only 10 real roots and 10 complex roots and these roots are completely different from 1,2,3,...,19,20. Therefore, we need to be careful about how we interpret the output from the methods below and would like to have additional information gathered (mathematically) to support the final output of each method.

Problem: Given a continuous function f(x), determine a value p such that f(p)=0 (approximately). Use IVT...

BISECTION METHOD

Assume you have real numbers a and b, with a<b such that f(a) * f(b) < 0. Then, IVT imples a root of f(x) exists in the interval [a,b]. Guess this root to be the midpoint of the interval c = (a+b)/2. Notice, this "bisects" the original interval into two intervals, each half the size of the first. Reapply this procedure on the piece which is guaranteed to contain a root (by IVT again).

In general, consider the sequences a_n, p_n, and b_n given by a₀=a, b₀=b, p_n+1=(a_n+b_n)/2 and a_n+1 and b_n+1 chosen to be the endpoints of the half of the interval that is kept. (ie, either a_n+1= a_n and b_n+1= p_n+1, OR a_n+1=p_n+1 and b_n+1=b_n.)

Stopping Criterion for Bisection:  How do you determine when this recursive process should end? How can you determine if f(p_n+1)=0 in the presence of roundoff?

• Error bound after n steps is | x^*-p_n+1 | < (b-a) 2^-n. Hence, if you can choose the number of steps n such that (b-a) 2^-n < e, for some small e, then p_n+1 must be at least e-close to x^*.
• Notice, one also could perhaps stop the iteration when | f(p_n+1) | < tolerance.
• Should check sign( f(a) ) sign( f(p_n) ), instead of computing f(a) f(p_n).

HOMEWORK: page 29, #8

Review: Taylor's Theorem

NEWTON'S METHOD - pg 35

Suppose p is the unknown actual root of f(x)=0 and p₀ is an approximation to p.  From Taylor's Theorem, f(p) = f(p₀) + (p-p₀)f'(p₀) + (p-p₀)²f''(c)/2!, where c is unknown but between p and p₀.
If p₀ is close to p, p-p₀ is "small" and (p-p₀)² is even "smaller".
Then 0 = f(p) @ f(p₀) + (p-p₀)f'(p₀), or by solving for p, we obtain p @ p₀ - f(p₀)/f'(p₀).

Comments on Newton's method:

• Derived graphically using a triangle with hypotenuse given by the tangent line to (a,f(a)).
• Newton's method is not as robust as Bisection but converges much faster. See page 37, 38.
• Online Java applet describing Newton's method is located at http://cs.jsu.edu/mcis/faculty/leathrum/JavaMath/topframe.html

From Taylor's theorem above, we have 0 = f(p) = f(p_k) + (p-p_k)f'(p_k) + (p-p_k)²f''(c)/2!, or
(p-p_k) + f(p_k) / f'(p_k) = -(p-p_k)²f''(c) / {2 * f'(p_k) }, or
p - p_k+1 = -(p-p_k)² [ f''(c) / {2 * f'(p_k) } ]. Hence, if [ f''(c) / {2 * f'(p_k) } ] < M is bounded near where f(p)=0, then we have | p - p_k+1| < M |p-p_k|² , which says that the iteration must converge as it progresses, so long as the bound exists.
Therefore, for a stopping criteria, be sure to continue the iteration until successive guesses p_k and p_k are very close to each other.

SECANT METHOD - pg 30

• Use Newton's by replacing the explicit evaluation of f'(a) with a difference quotient approximation.
• Requires no derivative but is not as fast as Newton's and needs two starting values instead of one.

Mixes the ideas of the Secant Method and the Bisection Method.

Convergence of this method can be hard to determine and even slower than Bisection.
 

RATES OF CONVERGENCE: If p is the exact value being successively approximated by the values p_n, then the error at step n is denoted by e_n = | p - p_n |.  The rate of convergence is the constant r such that e_n+1 < M e_n^r, for some constant M, or for which the ratio e_n+1 / e_n^r  is approximately a constant.
Special Cases:

• Linear Convergence | p - p_n+1 | < M | p - p_n |, for some constant M.   (Note the error formula for Bisection method.)
• Quadratic Convergence | p - p_n+1 | < M | p - p_n |², for some constant M.  (Note the error formula for Newton'smethod.)

HORNER'S ALGORITHM: Each calculation introduces additional roundoff and computer time. To minimize these for a given x value, rearrange the computations using a judicious choice of nested parenthesis. (Equivalent to Synthetic Substitution)
 

MUELLER'S METHOD: Takes 3 initial guesses p₀, p₁ and p₂, passes a quadratic through them and determines the zero closest to p₂. The quadratic: p(x) = a(p- p₂)² + b(p- p₂) + c. Plug in p₀, p₁ and p₂ to get a, b, and c using the formulas on page 46. Use the new improved roundoff quadratic-formula formulas and choose sign to make bottom largest possible.

• requires a square root at each step => possibly more expensive to implement
• only requires original function, not derivative
• can approximate complex roots easily by using complex arithmetic and complex square roots.
• is third order when original function is a polynomial