MATH/APPM 4650 Notes: February 8, 2008

MATH/APPM 4650
Class notes, February 8

Root-Finding Review

We have studied two types of methods: bounding methods (bisection and false position) and iterative methods (Newton and secant).

The basic idea of a bounding method is to start with an interval in which a solution is known. (Using the Intermediate Value Theorem.)
Then move one or both endpoints closer to the other, while still bounding the root.
Advantages:

Always have both an upper and lower estimate for the root, so error is known precisely.
Always converges given valid initial endpoints. (For bisection, both endpoints converge to the root. For false position, at least one endpoint converges.)

Disadvantages:

Usually slower to converge than an iterative method.
Requires two initial conditions, which must bound a root.

Rate of convergence:

Bisection will always converge linearly with α = 1 and λ = 1/2. Both endpoints converge at the same rate.
False position will usually have one endpoint converging and one endpoint not moving. The moving endpoint will usually converge linearly with α = 1 and λ = 1 - f'(p)/m, where m is the slope of the secant between the unmoving endpoint and the actual root.

The basic idea of an iterative method is to transform the root problem f(p) = 0 into the fixed-point problem g(p) = p.
If g(x) = x + h( f(x) ) for some one-to-one function h with h(0) = 0, then this is equivalent.
If h is well-chosen, the recursive sequence p_n+1 = g(p_n) converges to p as long as p₀ is sufficiently close to p.
Advantages:

Lots of freedom to choose g.
If it converges, it will often converge faster than a bounding method.
Requires only one starting approximation.

Disadvantages:

No guarantee it will converge to the desired root for a given function and initial condition.
No estimate of the error at the end.

An iterative method may be one-step (like Newton's method) or multi-step (like the secant method).

Error analysis:
For a one-step method, it is essential that g be chosen so that at the fixed point, |g'(p)|<1. There are usually many ways to do this.

If g'(p) ≠ 0, then a one-step iterative method is order α = 1 with asymptotic error constant λ = |g'(p)|.
If g'(p) = 0 but g''(p) ≠ 0, then the method is order α = 2. Since α > 1, we don't care what λ is. (This is the case for Newton's method, when the original function f has f'(p) ≠ 0 and f''(p) ≠ 0.)
If g'(p) = 0 and g''(p) = 0, the method will be at least order 3. Such algorithms are rare. (This is the case for Newton's method when f has an inflection point at the root.)

For a two-step method, the analysis is substantially more involved. The most popular two-step method is the secant method, which has order α = 1.618... as long as f'(p) ≠ 0 and f''(p) ≠ 0. This can be derived from the estimate p_n+1 - p = [f''(p)/(2f'(p))] (p_n - p) (p_n-1 - p).

The textbook mentions a three-step method called Muller's method, which works the same way as the secant method. Instead of finding the secant line joining the two previous (x, y) values, we find a parabola joining the three previous (x, y) values. Then we find where this parabola crosses the axis to get the new x value. This method involves more computation at each step, but the method is order α = 1.839...

Guidelines for using these methods:

If f and f' can both be computed to good accuracy, and if you don't need an error formula, Newton's method is best. This is true for all polynomials, and for simple functions if you have Maple or Mathematica.
If f is known exactly but f' is not known (for example, you have a formula for f but you're using MATLAB), and if you don't need an error formula, the secant method is a good compromise.
If f is not known exactly, you may want to use the bisection method. This is because bisection depends only on knowing the sign of f at any point, not on the explicit values of f.
If you need not just an estimate of the root, but also an interval containing the root (so that you know the error), you should use bisection.

The best methods (and the ones implemented in Maple, Mathematica, and MATLAB by default) use a combination of these methods. For example, if the function is known to be concave up at the root, you can use the secant method to get the left bound and Newton's method to get the right bound. (See example on the board.)

Although all mathematical software has a combination of root-finding algorithms, they can be unpredictable, and there may be good reasons to program your own.

For example, all standard algorithms will have trouble with the following function. Bisection because the function hardly changes sign at all; Newton's method because the derivative is close to zero near the roots; and secant method because there are equal function values. The latter two are likely to suffer from divide-by-zero errors or roundoff errors.

Acceleration of convergence:

All of the algorithms above will give a sequence (hopefully) converging to the root.
If we are stuck with slow (linear) convergence, we may want to use Aitken's acceleration trick.
This says to replace three terms of a sequence (p_n) which is known to converge linearly with the new sequence p^_n = p_n - [p_n+1 - p_n]² / [p_n+2 - 2 p_n+1 + p_n]. If p_n converges linearly with some asymptotic error constant λ, then we know p^_n will either converge linearly with a smaller asymptotic error constant, or will converge with higher order.

If p_n is known to come from iteration of a function, then the accelerated sequence will converge quadratically.

We can derive other acceleration formulas for sequences known to converge quadratically or with some other order.

However, it is important to realize that acceleration can create roundoff errors.
If p_n+2, p_n+1, and p_n are all close to the limit p, then the numerator and denominators of the fraction will both be close to zero. Thus we expect roundoff or divide-by-zero errors if we get too close to the answer.

Acceleration should only be used on the bad terms of a sequence, while they are still somewhat far from the limit. Make sure you have a stopping criterion to avoid dividing by very small numbers.

Steffensen's method is a general name for a broad class of methods, involving the combination of acceleration with a standard root-finding method.
In general, the technique is:

Start with an approximation.
Use your favorite linearly-convergent method to generate two more steps from the initial approximation.
Use Aitken's acceleration to find a new approximation from these three terms.
Go back to the beginning.

When used with a linearly convergent method, Steffensen's method gives a quadratically convergent sequence. (But it requires more computations per step.)

Finally, we briefly discuss Horner's method. This actually has nothing to do with root-finding.
It's just a way to compute polynomials efficiently.

Recall from Section 1.2 that when evaluating a polynomial f(x) = a_n xⁿ + a_n-1 x^n-1 + ... + a₁ x + a₀, there are roughly n²/2 multiplications and n additions if you do it directly.

This is a waste, because once you've computed x², you might as well take advantage of this to compute x³. (For example.)

Instead one should factor terms out of the polynomial.

f(x) = [[...[a_n x + a_n-1] x + ... + a₁] x + a₀. This is algebraically equivalent, but computationally much faster, because it only involves n multiplications and n additions. It's also less likely to result in roundoff error.

Horner's method is basically just this idea, along with the observation that you can use the same method to compute the derivative. It gives a slightly faster method to compute both f(x₀) and f'(x₀) at the same time. This is useful for Newton's method, which requires both to be evaluated at the same time.

See board for details.