Root Finding

1-D Root Finding

• Given some function, find location where f(x)=0
• Need:
• Starting position x₀, hopefully close to solution
• Ideally, points that bracket the root

f(x₊) > 0

f(x_–) < 0

1-D Root Finding

• Given some function, find location where f(x)=0
• Need:
• Starting position x₀, hopefully close to solution
• Ideally, points that bracket the root
• Well-behaved function

What Goes Wrong?

Tangent point:

very difficult�to find

Singularity:

brackets don’t�surround root

Pathological case:

infinite number of�roots – e.g. sin(1/x)

Example: Press et al., Numerical Recipes in C�Equation (3.0.1), p. 105

2

2 ln((Pi - x) )

f := x -> 3 x + ------------- + 1

4

Pi

> evalf(f(Pi-1e-10));

30.1360472472915692...

​

> evalf(f(Pi-1e-100));

25.8811536623525653...

​

> evalf(f(Pi-1e-600));

2.24285595777501258...

​

> evalf(f(Pi-1e-700));

-2.4848035831404979...

​

> evalf(f(Pi+1e-700));

-2.4848035831404979...

​

> evalf(f(Pi+1e-600));

2.24285595777501258...

​

> evalf(f(Pi+1e-100));

25.8811536623525653...

​

> evalf(f(Pi+1e-10));

30.1360472510614803...

​

Bisection Method

• Given points x₊ and x_– that bracket a root, find� x_half = ½ (x₊+ x_–)�and evaluate f(x_half)
• If positive, x₊ ← x_half else x_– ← x_half
• Stop when x₊ and x_– close enough
• If function is continuous, this will succeed�in finding some root

Bisection

• Very robust method
• Convergence rate:
• Error bounded by size of [x₊… x_–] interval
• Interval shrinks in half at each iteration
• Therefore, error cut in half at each iteration:� |ε_n+1| = ½ |ε_n|
• This is called “linear convergence”
• One extra bit of accuracy in x at each iteration

Faster Root-Finding

• Fancier methods get super-linear convergence
• Typical approach: model function locally by something whose root you can find exactly
• Model didn’t match function exactly, so iterate
• In many cases, these are less safe than bisection

Secant Method

• Simple extension to bisection: interpolate or extrapolate through two most recent points

1

2

3

4

Secant Method

• Faster than bisection:� |ε_n+1| = const. |ε_n|^1.6
• Faster than linear: number of correct bits multiplied by 1.6
• Drawback: the above only true if sufficiently close to a root of a sufficiently smooth function
• Does not guarantee that root remains bracketed

False Position Method

• Similar to secant, but guarantee bracketing

​

​

​

​

​

• Stable, but linear in bad cases

1

2

3

4

Other Interpolation Strategies

• Ridders’s method: fit exponential to�f(x₊), f(x_–), and f(x_half)
• Van Wijngaarden-Dekker-Brent method:�inverse quadratic fit to 3 most recent points�if within bracket, else bisection
• Both of these safe if function is nasty, but�fast (super-linear) if function is nice

Newton-Raphson

• Best-known algorithm for getting quadratic�convergence when derivative is easy to evaluate
• Another variant on the extrapolation theme

1

2

3

4

Slope = derivative at 1

Newton-Raphson

• Begin with Taylor series

​

• Divide by derivative (can’t be zero!)

Newton-Raphson

• Method fragile: can easily get confused

​

​

​

​

• Good starting point critical
• Newton popular for “polishing off” a root found approximately using a more robust method

Newton-Raphson Convergence

• Can talk about “basin of convergence”:�range of x₀ for which method finds a root
• Can be extremely complex:�here’s an example�in 2-D with 4 roots

Yale site

D.W. Hyatt

Popular Example of Newton: Square Root

• Let f(x) = x² – a: zero of this is square root of a
• f’(x) = 2x, so Newton iteration is

​

​

• “Divide and average” method

Reciprocal via Newton

• Division is slowest of basic operations
• On some computers, hardware divide not available (!): simulate in software

​

​

​

​

• Need only subtract and multiply

Rootfinding in >1D

• Behavior can be complex: e.g. in 2D

Rootfinding in >1D

• Can’t bracket and bisect

​

​

​

​

​

​

• Result: few general methods

Newton in Higher Dimensions

• Start with

​

​

• Write as vector-valued function

Newton in Higher Dimensions

• Expand in terms of Taylor series

​

​

• f’ is a Jacobian

Newton in Higher Dimensions

• Solve for δ

​

​

• Requires matrix inversion (we’ll see this later)
• Often fragile, must be careful
• Keep track of whether error decreases
• If not, try a smaller step in direction δ