AMTH247 Lecture 15

Nonlinear Equations II

Reading: Heath §5.5.1, 5.5.3, 5.5.4, 5.5.7

Contents

• Algorithms
  • Bisection Algorithm
  • Newton's Method
  • Secant Method
  • Practical Algorithms
  
  
• About this document ...

Notes: Again, there will be no formal practical this week.

Algorithms

Bisection Algorithm

This algorithm is based on the simple idea that if a continuous function $f(x)$ has opposite signs at two points $a$ and $b$, then it must have a root on the interval $[a, b]$.

Take the midpoint $m$ of the interval $[a, b]$. Depending on the sign of $f(m)$ we can determine in which of the subintervals $[a, m]$ and $[m, b]$ the root lies. Repeating this process we get the bisection algorithm.

The bisection algorithm is easy to implement in Scilab. The following function computes the midpoint of the current interval [a,b], tests the sign of f(m), and adjusts the endpoints until the length of the interval is less than tol.

function [a,b] = bisect (f, a0, b0, tol)
  if (sign(f(a0)) == sign(f(b0)))
    error ("function has same sign at both endpoints")
  end
  a = a0
  b = b0;
  while (abs(a-b) > tol)
    m = a + (b-a)/2
    if (sign(f(m)) == sign(f(a)))
      a = m
    else
      b = m
    end
  end
endfunction

We will use the example of Kepler's equation from the previous lecture

$$f(E) = E - 1 - 0.5 \sin(E) = 0 .$$

-->function y = kepler (e)
-->   y = e - 1 - 0.5*sin(e)
-->endfunction

-->[a,b] = bisect(kepler, 1, 2, 1e-5)
 b   =

    1.498703
 a   =

    1.4986954

-->[a,b] = bisect(kepler, 1, 2, 1e-10)
 b   =

    1.4987011
 a   =

    1.4987011

It is easy to see that the bisection algorithm is linearly convergent with rate of convergence $1/2$, since at each iteration the size of the interval containing the root is decreased by a factor of 2. Further, it is easy to calculate the number of iterations needed to achieve a given accuracy. Starting with an interval $[a, b]$, after $k$ iterations the length of the interval containing the root is $(b - a)/2^k$, so to obtain an interval of length $t$ requires

$$\log_2 \left( \frac{b-a}{t} \right)$$

iterations.

One important property of the bisection algorithm, not shared by most other algorithms, is that it always converges provided the function is continuous and we start from interval on which the function has opposite signs at the two endpoints.

Newton's Method

The idea behind Newton's method is equally simple. Given an approximation $x_k$ to a solution of $f(x) = 0$, approximate $f(x)$ by its tangent at $x_k$ and take the point at which the tangent intersects the $x$-axis as the next approximation $x_{k+1}$

A simple geometric calculation then gives

$$x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)}$$

which is the iteration formula for Newton's method.

Newton's method is easy to program in Scilab. The following function performs n iterations of Newton's method starting with the initial approximation x0. Note that we need the derivative of the function $f(x)$ as well as the function itself.

function x = newton (x0, f, df, n)
  x = x0
  for k = 1:n
    x = x - f(x)/df(x)
  end
endfunction

Applying this to Kepler's equation gives

-->function dy = dkepler (e)
-->   dy = 1 - 0.5*cos(e)
-->endfunction

-->x = newton(1.5, kepler, dkepler, 5)
 x   =

    1.4987011

Newton's method can be viewed as a fixed point iteration

$$x_{k+1} = g(x_k)$$

with

$$g(x) = x - \frac{f(x)}{f'(x)}$$

The derivative of $g(x)$ is

$$g'(x) = \frac{f(x)f''(x)}{f'(x)^2}$$

A solution of $f(x) = 0$ for which $f'(x) = 0$ is called a multiple root of $f(x)$. When $f(x) = 0$, and provided $f'(x) \neq 0$, we also have $g'(x) = 0$. So, by the results of the previous lecture, Newton's method is quadratically convergent except for multiple roots. For a multiple root it turns out that Newton's method is only linearly convergent.

To examine the convergence of Newton's method, we can modify our function newton so that it returns the sequence of iterations:

function x = newtons (x0, f, df, n)
  x = zeros(1,n+1)
  x(1) = x0
  for k = 1:n
    x(k+1) = x(k) - f(x(k))/df(x(k))
  end
endfunction

For our example:

-->x = newtons(1.5, kepler, dkepler, 5)
 x   =

          column 1 to 5

!    1.5     1.4987016     1.4987011     1.4987011     1.4987011 !

          column 6

!    1.4987011 !

Comparing the iterations with the final value

-->x-x(6)
 ans   =

!    0.0012989     4.361E-07     4.907E-14     0.     0.     0. !

shows that it has converged to full accuracy in 3 iterations! The quadratic convergence, that is the error at each iteration is proportional to the square of the error at the previous iteration, is also apparent.

Beware that Newton's method, like other fixed point methods, does not necessarily converge to a solution. To ensure convergence we need to start with an initial approximation sufficiently close to the solution.

Secant Method

The secant method is similar to Newton's method except that we approximate the function by a secant rather than by the tangent. Given two approximations $x_{k-1}$ and $x_k$ to a solution of $f(x) = 0$, approximate $f(x)$ by its secant at $x_{k-1}$ and $x_k$ and take the point at which the secant intersects the $x$-axis as the next approximation $x_{k+1}$.

Again a simple geometric calculation gives

$$x_{k+1} = x_k - \frac{f(x_k)(x_k - x_{k-1})}{f(x_k) - f(x_{k-1})} .$$

which is the iteration formula for the secant method. Another way to get this formula is to replace $f'(x)$ in Newton's method by the approximation

$$f'(x) \approx \frac{f(x_k) - f(x_{k-1})}{x_k - x_{k-1}} .$$

Below is a Scilab function for the secant method. It stores the iterations and takes two initial approximations x0 and x1 to get started.

function x = secant (x0, x1, f, n)
  x = zeros(1, n+2)
  x(1) = x0
  x(2) = x1
  for k = 2:n+1
    x(k+1) = x(k) - f(x(k))*(x(k) - x(k-1))/(f(x(k)) - f(x(k-1)))
  end
endfunction

For Kepler's equation it gives

-->x = secant(2.0, 1.5,   kepler, 5)
 x   =

          column 1 to 6

!    2.     1.5     1.498849     1.4987012     1.4987011     1.4987011 !

          column 7

!    1.4987011 !

We have to a bit careful using the function secant since it will result in a division by zero when two successive iterates are the same, a situation which will occur whenever it has converged to full accuracy.

The secant algorithm has order of convergence $r \approx 1.62$ (see Heath, p223) except at multiple roots when, like Newton's method, it is only linearly convergent. Again, like Newton's method, it requires good initial approximations to ensure convergence.

Practical Algorithms

The essential properties of the algorithms we have studied can be summarized as follows

Method	Always	Order of	Requires
	Converges	Convergence	Derivatives
Bisection	Yes	Linear	No
Secant	No	1.62	No
Newton	No	Quadratic	Yes

It is clear that each of these algorithms has its advantages and disadvantages. The methods used in practice, for example Scilab's fsolve, use a combination of methods to obtain the advantages of each. These algorithms tend to be quite sophisticated and rely on years of accumulated experience.

A simple way to combine the bisection method with other more rapidly convergent algorithms is to start with an interval known to contain the solution and apply the faster method as long as it produces results in the interval. When that fails, bisection can be applied to reduce the width of the interval and then the faster method tried again.

Scilab's fsolve is easy to use, it simply requires an initial approximation

-->x = fsolve(1.5, kepler)
 x  =
 
    1.4987011

About this document ...

AMTH247 Lecture 15

Nonlinear Equations II

This document was generated using the LaTeX2HTML translator Version 2K.1beta (1.61)

Copyright © 1993, 1994, 1995, 1996, Nikos Drakos, Computer Based Learning Unit, University of Leeds.
Copyright © 1997, 1998, 1999, Ross Moore, Mathematics Department, Macquarie University, Sydney.

The command line arguments were:
latex2html -split +0 lecture.tex

The translation was initiated by amth247 on 2003-05-13