Previous Contents Next

Chapter 5  Root Finding and Simple Optimization

5.1  Root Finding:

Let's discuss root finding using both the bisection and Newton method.

5.1.1  The Bisection Method:

We need a simple function to find the root of a nice function f of the real variable x using what is called bisection. The method is actually quite simple. We know that if f is a continuous function on the finite interval [a,b] then f must have a zero inside the interval [a,b] if f has a different algebraic sign at the endpoints a and b. This means the product f(a) f(b) is not zero. So we assume we can find an interval [a,b] on which this change in sign satisfies f(a) f(b) £ 0 (which we can do by switching to -f if we have to!) and then if we divide the interval [a,b] into two equal pieces [a,m] and [m,b], f(m) can't have the same sign as both f(a) and f(b) because of the assumed sign difference. So at least one of the two halves has a sign change.

Note that if f(a) and f(b) was zero then we still have f(a) f(b) £ 0 and either a or b could be our chosen root and either half interval works fine. If only one of the endpoint function values is zero, then the bisection of [a,b] into the two halves still finds the one half interval that has the root.

So our prototyping Matlab code should use tests like f(x) f(y) £ 0 rather than f(x) f(y) < 0 to make sure we catch the root.

The Bisection Matlab Code:

Here is a simple Matlab function to perform the Bisection routine.

Bisection.m 
 1 function root = Bisection(fname,a,b,delta)
 2 %
 3 % fname is a string that gives the name of a 
 4 %       continuous function f(x) 
 5 % a,b   this is the interval [a,b] on
 6 %       which f is defined and for which
 7 %       we assume that the product
 8 %       f(a)*f(b) <= 0
 9 % delta this is a non negative real number
10 %
11 % root  this is the midpoint of the interval
12 %       [alpha,beta] having the property that
13 %       f(alpha)*f(beta) <= 0 and
14 %       |beta-alpha| <= delta + eps * max(|alpha|,|beta|)
15 %       and eps is machine zero
16 %
17   disp(' ')
18   disp('      k  |       a(k)    |      b(k)     |       b(k) - a(k) ')
19   k = 1;
20   disp(sprintf(' %6d  | %12.7f  | %12.7f  | %12.7f',k,a,b,b-a));           
21   fa = feval(fname,a);
22   fb = feval(fname,b);
23   while abs(a-b) > delta + eps*max(abs(a),abs(b))
24     mid = (a+b)/2;
25     fmid = feval(fname,mid);
26     if fa*fmid <= 0
27       % there is a root in [a,mid]
28       b = mid;
29       fb = fmid;
30     else
31       % there is a root in [mid,b]
32       a = mid;
33       fa = fmid;
34     end
35     k = k+1;
36     disp(sprintf(' %6d  | %12.7f  | %12.7f  | %12.7f',k,a,b,b-a));
37   end
38   root = (a+b)/2;
syntax highlighted by Code2HTML, v. 0.8.8b

We should look at some of these lines more closely. First, to use this routine, we need to write a function definition for the function we want to apply bisection to. We will do this in a file called func.m (Inspired Name, eh?) An example would be the one we wrote for the function
f(x) =
tan(
x
4
) - 1;

which is coded in Matlab by

func.m 
1 function y = func(x)
2 %
3 % x  real input
4 % y  real output
5 %
6   y = tan(x/4)-1;
syntax highlighted by Code2HTML, v. 0.8.8b

So to apply bisection to this function on the interval [2,4] with a stopping tolerance of say 10-4, in Matlab, we would type the command 
 root = Bisection('func',2,4,10^-4)}
Note that the name of our supplied function, the uninspired choice func is passed in as the first argument in single quotes as it is a string.

Also, in the Bisection routine, we have added the code to print out what is happening at each iteration of the while loop. Matlab handles prints to the screen a little funny, so do set up a table of printed values we use this syntax: 
  % this prints a blank line and then a table heading.
  % note disp prints a string only
  disp(' ')
  disp('      k  |       a(k)    |      b(k)     |       b(k) - a(k) ')
  % now to print the k, a, b and b-a, we must first
  % put their values into a string using the c like function
  % sprintf and then use disp to disply that string.
  % so we do this 
  % disp( sprintf(' output specifications here ',variables here))
  % so inside the while loop we use
  disp(sprintf(' %6d  | %12.7f  | %12.7f  | %12.7f',k,a,b,b-a));

Running the Code:

As mentioned above, we will test this code on the function
f(x) =
tan(
x
4
) - 1;

on the interval [2,4] with a stopping tolerance of d = 10-6. Our function has been written as the Matlab function func supplied in the file func.m.

The Matlab run time looks like this:

Bisection.prompts 
 1 >> eps
 2 
 3 ans =
 4 
 5    2.2204e-16
 6 >> root = Bisection('func',2,4,10^-6)
 7  
 8       k  |       a(k)    |      b(k)     | b(k) - a(k) 
 9       1  |    2.0000000  |    4.0000000  |    2.0000000
10       2  |    3.0000000  |    4.0000000  |    1.0000000
11       3  |    3.0000000  |    3.5000000  |    0.5000000
12       4  |    3.0000000  |    3.2500000  |    0.2500000
13       5  |    3.1250000  |    3.2500000  |    0.1250000
14       6  |    3.1250000  |    3.1875000  |    0.0625000
15       7  |    3.1250000  |    3.1562500  |    0.0312500
16       8  |    3.1406250  |    3.1562500  |    0.0156250
17       9  |    3.1406250  |    3.1484375  |    0.0078125
18      10  |    3.1406250  |    3.1445312  |    0.0039062
19      11  |    3.1406250  |    3.1425781  |    0.0019531
20      12  |    3.1406250  |    3.1416016  |    0.0009766
21      13  |    3.1411133  |    3.1416016  |    0.0004883
22      14  |    3.1413574  |    3.1416016  |    0.0002441
23      15  |    3.1414795  |    3.1416016  |    0.0001221
24      16  |    3.1415405  |    3.1416016  |    0.0000610
25      17  |    3.1415710  |    3.1416016  |    0.0000305
26      18  |    3.1415863  |    3.1416016  |    0.0000153
27      19  |    3.1415863  |    3.1415939  |    0.0000076
28      20  |    3.1415901  |    3.1415939  |    0.0000038
29      21  |    3.1415920  |    3.1415939  |    0.0000019
30      22  |    3.1415920  |    3.1415930  |    0.0000010
31 
32 root =
33 
34     3.1416
syntax highlighted by Code2HTML, v. 0.8.8b

Exercises:

Well, you have to practice this stuff to see what is going on. So here are two problems to sink your teeth into!
  1. Use bisection to find the first five positive solutions of the equation x = tan(x). You can see where this is roughly by graphing tan(x) and x simultaneously. Do this for tolerances {10-1, 10-2, 10-3, 10-4, 10-5, 10-6, 10-7}. For each root, choose a reasonable bracketing interval [a,b], explain why you chose it, provide a table of the number of iterations to achieve the accuracy and a graph of this number vs. accuracy.
  2. Use the Bisection Method to find the largest real root of the function f(x) = x6 - x - 1. Do this for tolerances {10-1, 10-2, 10-3, 10-4, 10-5, 10-6, 10-7}. Choose a reasonable bracketing interval [a,b], explain why you chose it, provide a table of the number of iterations to achieve the accuracy and a graph of this number vs. accuracy.

5.1.2  Newton's Method:

When Do We Do A Newton Step?

The following code uses a simple test to see if we should use a bisection or a newton method step in our zero finding routine.

StepIsIn.m 
 1 function ok = StepIsIn(x,fx,fpx,a,b)
 2 %
 3 % x       current approximate root
 4 % fx      value of f at the approximate root
 5 % fpx     value of derivative of f at the approximate
 6 %         root
 7 % a,b     the interval the root is in.
 8 %
 9 % ok      1 if the Newton Step x - fx/fpx is in [a,b]
10 %         0 if not
11 %
12 if fpx > 0
13   ok = ((a-x)*fpx <= -fx) & (-fx <= (b-x)*fpx);
14 elseif fpx < 0
15   ok = ((a-x)*fpx >= -fx) & (-fx >= (b-x)*fpx);
16 else
17   ok = 0;
18 end
syntax highlighted by Code2HTML, v. 0.8.8b

A Global Newton Method:

GlobalNewton.m 
 1 function [x,fx,nEvals,aF,bF] = ...
 2          GlobalNewton(fName,fpName,a,b,tolx,tolf,nEvalsMax)
 3 %
 4 % fName        a string that is the name of the function f(x)
 5 % fpName       a string that is the name of the functions derivative
 6 %              f'(x)
 7 % a,b          we look for the root in the interval [a,b]
 8 % tolx         tolerance on the size of the interval
 9 % tolf         tolerance of f(current approximation to root)
10 % nEvalsMax    Maximum Number of derivative Evaluations
11 %
12 % x            Approximate zero of f
13 % fx           The value of f at the approximate zero
14 % nEvals       The Number of Derivative Evaluations Needed
15 % aF, bF       the final interval the approximate root lies in,
16 %              [aF,bF]
17 %
18 % Termination  Interval [a,b] has size < tolx
19 %              |f(approximate root)| < tolf
20 %              Have exceeded nEvalsMax derivative Evaluations
21 %
22 fa = feval(fName,a);
23 fb = feval(fName,b);
24 x = a;
25 fx = feval(fName,x);
26 fpx = feval(fpName,x);
27 
28 nEvals = 1;
29 k = 1;
30 disp(' ')
31 disp(' Step     |       k  |       a(k)    |      x(k)     |       b(k) ')
32 disp(sprintf('Start     |  %6d  | %12.7f  | %12.7f  | %12.7f',k,a,x,b));
33 while   (abs(a-b)>tolx) & (abs(fx)> tolf) & (nEvals<nEvalsMax)| (nEvals==1)
34   %[a,b] brackets a root and x=a or x=b
35   check = StepIsIn(x,fx,fpx,a,b);
36   if check
37     %Take Newton Step
38     x = x - fx/fpx;
39   else
40     %Take a Bisection Step:
41     x = (a+b)/2;
42   end
43   fx = feval(fName,x);
44   fpx = feval(fpName,x);
45   nEvals = nEvals+1;
46   if fa*fx<=0
47     %there is a root in [a,x].  Use right endpoint.
48     b = x;
49     fb = fx;
50   else
51     %there is a root in [x,b].  Bring in left endpoint.
52     a = x;
53     fa = fx;
54   end
55   k = k+1;
56   if(check)
57     disp(sprintf('Newton    |  %6d  | %12.7f  | %12.7f  | %12.7f',k,a,x,b));
58   else
59     disp(sprintf('Bisection |  %6d  | %12.7f  | %12.7f  | %12.7f',k,a,x,b));
60   end
61 end
62 aF = a;
63 bF = b;
syntax highlighted by Code2HTML, v. 0.8.8b

A Run Time Example:

We will apply our global newton method root finding code to a simple example: find a root for f(x) = sin(x) in the interval [-7 p/2, 15 p + 0.1]. We code the function and its derivative in two simple Matlab files; f1.m and f1p.m. These are

f1.m 
1 function y = f1(x)
2 y = sin(x);
syntax highlighted by Code2HTML, v. 0.8.8b

and

f1p.m 
1 function y = f1p(x)
2 y = cos(x);
syntax highlighted by Code2HTML, v. 0.8.8b

To run this code on this example, we would then type a phrase like the one below: 
[x,fx,nEvals,aLast,bLast] = GlobalNewton('f1','f1p',-7*pi/2,15*pi+.1,...
                                         10^-6,10^-8,200)


Here is the runtime output:

GlobalMethod.prompts 
 1 >> [x,fx,nEvals,aLast,bLast] = GlobalNewton('f1','f1p',-7*pi/2,15*pi+.1,...
 2                                             10^-6,10^-8,200)
 3  
 4  Step     |       k  |       a(k)    |      x(k)     |       b(k) 
 5 Start     |       1  |  -10.9955743  |  -10.9955743  |   47.2238898
 6 Bisection |       2  |  -10.9955743  |   18.1141578  |   18.1141578
 7 Bisection |       3  |  -10.9955743  |    3.5592917  |    3.5592917
 8 Newton    |       4  |    3.1154761  |    3.1154761  |    3.5592917
 9 Newton    |       5  |    3.1154761  |    3.1415986  |    3.1415986
10 Newton    |       6  |    3.1415927  |    3.1415927  |    3.1415986
11 
12 x =
13 
14     3.1416
15 
16 
17 fx =
18 
19    1.2246e-16
20 
21 
22 nEvals =
23 
24      6
25 
26 
27 aLast =
28 
29     3.1416
30 
31 
32 bLast =
33 
34     3.1416
syntax highlighted by Code2HTML, v. 0.8.8b

Some Exercises:

  1. Use the Global Newton Method to find the first five positive solutions of the equation x = tan(x). You can see where this is roughly by graphing tan(x) and x simultaneously. Do this for tolerances {10-1, 10-2, 10-3, 10-4, 10-5, 10-6, 10-7}. For each root, choose a reasonable bracketing interval [a,b], explain why you chose it, provide a table of the number of iterations to achieve the accuracy and a graph of this number vs. accuracy.
  2. Use the Global Newton Method to find the largest real root of the function f(x) = x6 - x - 1. Do this for tolerances {10-1, 10-2, 10-3, 10-4, 10-5, 10-6, 10-7}. Choose a reasonable bracketing interval [a,b], explain why you chose it, provide a table of the number of iterations to achieve the accuracy and a graph of this number vs. accuracy.

5.1.3  Adding Finite Difference Approximations to the Derivative:

We can also choose to replace the derivative function for f with a finite difference approximation. We will use
f' (x) approx
f(xc + dc) - f(xc)
dc

to approximate the value of the derivative at the point xc. As we have discussed earlier, some care is required to pick a size for dc so that round-off errors do not destroy the accuracy of our finite difference approximation to f' .

The simple Matlab code to implement this is given below: 
fval = feval(fname,x);
fpval = (feval(fname,x+delta) - fval)/delta;
We can also use a secant approximation as follows:
f' (x) approx
f(xc) - f(x-
xc - x-

where x- is the previous iterate from our routine. The Matlab fragment we need is then: 
fpc = (fc - f_)/(xc - x_);

A Finite Difference Global Newton Method:

We add the finite difference routines into our Global Newton's Method as follows:

GlobalNewtonFD.m 
 1 function [x,fx,nEvals,aF,bF] = ...
 2          GlobalNewton(fName,a,b,tolx,tolf,nEvalsMax)
 3 %
 4 % fName        a string that is the name of the function f(x)
 5 % a,b          we look for the root in the interval [a,b]
 6 % tolx         tolerance on the size of the interval
 7 % tolf         tolerance of f(current approximation to root)
 8 % nEvalsMax    Maximum Number of derivative Evaluations
 9 %
10 % x            Approximate zero of f
11 % fx           The value of f at the approximate zero
12 % nEvals       The Number of Derivative Evaluations Needed
13 % aF, bF       the final interval the approximate root lies in,
14 %              [aF,bF]
15 %
16 % Termination  Interval [a,b] has size < tolx
17 %              |f(approximate root)| < tolf
18 %              Have exceeded nEvalsMax derivative Evaluations
19 %
20 fa = feval(fName,a);
21 fb = feval(fName,b);
22 x = a;
23 fx = feval(fName,x);
24 delta = sqrt(eps)* abs(x);
25 fpval = feval(fName,x+delta);
26 fpx = (fpval-fx)/delta;
27 
28 nEvals = 1;
29 k = 1;
30 disp(' ')
31 disp(' Step     |       k  |       a(k)    |      x(k)     |       b(k) ')
32 disp(sprintf('Start     |  %6d  | %12.7f  | %12.7f  | %12.7f',k,a,x,b));
33 while   (abs(a-b)>tolx) & (abs(fx)> tolf) & (nEvals<nEvalsMax)| (nEvals==1)
34   %[a,b] brackets a root and x=a or x=b
35   check = StepIsIn(x,fx,fpx,a,b);
36   if check
37     %Take Newton Step
38     x = x - fx/fpx;
39   else
40     %Take a Bisection Step:
41     x = (a+b)/2;
42   end
43   fx = feval(fName,x);
44   fpval = feval(fName,x+delta);
45   fpx = (fpval-fx)/delta;  
46   nEvals = nEvals+1;
47   if fa*fx<=0
48     %there is a root in [a,x].  Use right endpoint.
49     b = x;
50     fb = fx;
51   else
52     %there is a root in [x,b].  Bring in left endpoint.
53     a = x;
54     fa = fx;
55   end
56   k = k+1;
57   if(check)
58     disp(sprintf('Newton    |  %6d  | %12.7f  | %12.7f  | %12.7f',k,a,x,b));
59   else
60     disp(sprintf('Bisection |  %6d  | %12.7f  | %12.7f  | %12.7f',k,a,x,b));
61   end
62 end
63 aF = a;
64 bF = b;
syntax highlighted by Code2HTML, v. 0.8.8b

Note, we use for our finite difference stepsize sqrt(emachine) |x|.

A Run Time Example:

We will apply our finite difference global newton method root finding code to the same simple example: find a root for f(x) = sin(x) in the interval [-7 p/2, 15 p + 0.1]. We only need the code for the function now which is as usual in the file f1.m.

To run this code on this example, we would then type a phrase like the one below: 
[x,fx,nEvals,aLast,bLast] = GlobalNewtonFD('f1',-7*pi/2,15*pi+.1,...
                                         10^-6,10^-8,200)


Here is the runtime output:

GlobalNewtonFD.prompts 
 1 >> [x,fx,nEvals,aLast,bLast] = GlobalNewtonFD('f1',-7*pi/2,15*pi+.1,...
 2                                             10^-6,10^-8,200)
 3  
 4  Step     |       k  |       a(k)    |      x(k)     |       b(k) 
 5 Start     |       1  |  -10.9955743  |  -10.9955743  |   47.2238898
 6 Bisection |       2  |  -10.9955743  |   18.1141578  |   18.1141578
 7 Bisection |       3  |  -10.9955743  |    3.5592917  |    3.5592917
 8 Newton    |       4  |    3.1154761  |    3.1154761  |    3.5592917
 9 Newton    |       5  |    3.1154761  |    3.1415986  |    3.1415986
10 Newton    |       6  |    3.1154761  |    3.1415927  |    3.1415927
11 
12 x =
13 
14     3.1416
15 
16 
17 fx =
18 
19   -4.3184e-15
20 
21 
22 nEvals =
23 
24      6
25 
26 
27 aLast =
28 
29     3.1155
30 
31 
32 bLast =
33 
34     3.1416
syntax highlighted by Code2HTML, v. 0.8.8b

Some Exercises:

  1. Use the Finite Difference Global Newton Method to find the second positive solution of the equation x = tan(x). Do this for tolerances 10-8. This time alter the GlobalNewtonFD code to allow the finite difference step size delta to be a parameter and do a parametric study on the effects of delta. Note that the code now uses the reasonable choice of sqrt(emachine) |x| but you need to use the additional d choices {10-4, 10-6, 10-8, 10-10}. This will give you five d choices. Provide a table and a graph of d vs. accuracy of the root approximation.
  2. Use the Finite Difference Global Newton Method to find the largest real root of the function f(x) = x6 - x - 1. Do this for tolerances 10-8. Again use altered GlobalNewtonFD code with the finite difference step size delta as a parameter and do a parametric study on the effects of delta. Note that the code now uses the reasonable choice of sqrt(emachine) |x| but you need to use the additional d choices {10-4, 10-6, 10-8, 10-10}. This will give you five d choices. Provide a table and a graph of d vs. accuracy of the root approximation.
  3. Do the same thing for the problems above, but replace the Finite Difference Global Newton Code with a Secant Global Newton Code. This will only require a few lines of code to change really, so don't freak out!
Previous Contents Next