LAB #3: Roots of Equations


TABLE OF CONTENTS

Introduction

The root finding problem is the task of finding one or more values for a variable x so that an equation

f(x) = 0
is satisfied. As a computational problem, we are interested in effective and efficient ways of finding one value x so that the residual error, that is,
abs ( f(x) )
is small. We hope this means that the approximation error is small. Today we will look at some nonlinear equations and a very simple method of finding an approximate solution. At the end of this session, there is an assignment.

A Sample Problem

Suppose we want to know if there is a solution to:

cos ( x ) = x
if the solution is unique, and what that solution is. This is not an algebraic equation, and there is little hope of forming an explicit expression for the solution.

Since we can't solve this equation exactly, it's worth knowing whether there is anything we can say about it. In fact, if there is a solution x* to the equation, we can probably find a computer representable number x which approximately satisfies the equation (has low residual error) and which is close to x* (has a low approximation error).

Exercise: Plot the functions cos(x) and x together in the same graph. What does your graph suggest about the solubility of our sample problem?

Exercise: write an M file cosy.m that defines the function

f(x) = cos ( x ) - x
(Remember how, in lab #2, we talked about how to set up a function?)
        function result = name ( input )
        result = ...
      

Note that MATLAB has a special command for plotting functions that's easier than what we normally do. For instance, we could plot cosy from -4.0 to 4.0 by typing:

fplot ( 'cosy', [ -4.0, 4.0 ] )
Note the quotes around the function name, and the square brackets that make the x limits a vector. Now would be a good time to try the fplot command out, and see what the cosy function looks like. Can you add an x axis using hold on and the appropriate plot command?

The Bisection Idea

The idea of the bisection method is very simple. We assume that we are given two values A and B, and that the function F(X) is positive at one of these values (say at A) and negative at the other. (When this occurs, we say that [A,B] is a change-of-sign interval for the function. So we know that (assuming F is continuous) that there must be one (at least one) root in the interval.

Consider the point

C = ( A + B ) / 2
If F(C)=0 we are done. (This is pretty unlikely). Otherwise, depending on the sign of F(C), we know that the root lies in [A,C] or [C,B]. In any case, our change-of-sign interval is now half as large as before. Repeat this process with the new change of sign interval, until the interval is "sufficiently small", say no larger than TOL. Declare victory.

We are guaranteed to converge. We can even compute the maximum number of steps this will take. We know in advance how well we will approximate the root X*. These are very powerful facts, which make bisection a robust algorithm - that is, it is very hard to defeat it.

Discussion: If we know the start points A and B and the interval size tolerance TOL, we can predict beforehand the number of steps the bisection code will (probably) take. What is this formula?

Discussion: Name a function which has a root in the interval [-1, 1], but for which bisection could not be used.

Programming Bisection

Before we look at a sample bisection program, let's discuss some things. If you haven't done much programming before, this is a good time to try to understand the logic behind how we choose variables to set up, what names to give them, and how to control the logic.

First of all, we can write the bisection algorithm as a MATLAB function. The input can be the two endpoints A and B, and we expect the output to be the approximation to the root. Inside of the bisection code, we should probably evaluate the function at both points, and save these values in FA and FB.

Discussion: Now think about the iteration. We control the iteration with a while loop, because we are going to keep doing it until something is true (and what is that?) I'm sure you can see how to compute C and FC, but what do we do to check the signs and update the interval. (The MATLAB sign(x) function is part of the answer).

A Bisection Program

Here is a sketch of a MATLAB function that carries out the bisection algorithm for our cosy function. The endpoints of the initial interval are variables that you specify when you invoke the function. Please note that this "sketch" is not yet correct! We'll make some corrections and improvements to it.

        function approx_root = bisect ( a, b )

        fa = cosy(a); 
        fb = cosy(b);

        while ( abs ( b - a ) > 0.000001 )

          c = ( a + b ) / 2;
          fc = cosy(c);

          %
          %  What follows is just a nice way to print out a little table:
          %
          [  a,  c,  b;
            fa, fc, fb ]

          if ( fc == 0 ) 
            approx_root = c;
            return
          elseif ( sign ( fc ) == sign ( fa ) )
            a = c;
            fa = fc;
          elseif ( sign ( fc ) == sign ( fb ) )
            b = c;
            fb = fc;
          end

        end
      

Discussion: there is a very serious thing wrong with this program. It almost never will give you the right answer! A single line is missing. What is it?

Exercise: Type in this function as the file bisect.m. Then try the command:

x = bisect ( 0, 3 )
Is the value x a root of the equation
cos ( x ) = x?

Discussion:

Variable Function Names

Suppose we have written an M file called bisect.m and we want to use it on some function. Then we have to create an M file called cosy.m, because that is what the bisect function expects. Suppose that we have 6 files, called f1.m through f6.m, each of which evaluates a function that we want to use with bisect. We could rename f1.m to cosy.m, and edit its contents, (because the name f1 occurs in the text of the file as well), and repeat this for the other files. But who wants to do that? Besides, cosy is a stupid name!

There is a much more convenient way to do this. The bisect routine treats the endpoints of the interval as variables, and gives them "dummy names" that we fill in when we use the routine. We can do the same thing with functions...almost. Here are the steps required:

Exercise: Follow the directions, and modify the file bisect.m so that the function name is a variable. Once you are done, see if the command

x = bisect ( 'cosy', 0, 3 )
is acceptable to MATLAB and produces the same root as before.

Some Nonlinear Functions

Here is a small collection of nonlinear functions. For each function f(x), we give a formula and the endpoints of a search interval.

FORMULA

f1(x) = x^2 - 9
f2(x) = x^6 - x - 1
f3(x) = 2 * x - exp ( - x )
f4(x) = ( x + 3 ) * ( x - 1 )^2
f5(x) = 20 * x / ( 100 * x^2 + 1 )
f6(x) = ( 16 - x^4 ) / ( 16 * x^4 )

INTERVAL

[ 0, 10]
[ 1, 2]
[ -10, 100]
[-1000, 1000]
[ -10, 10]
[ -1, 1000]

COMPUTED ROOT

____________
____________
____________
____________
____________
____________

NUMBER OF STEPS

______________
______________
______________
______________
______________
______________

Exercise: make 6 little function files, called f1.m and so on, to evaluate each of these functions. Take a look at the plots of these functions over the suggested intervals.

Assignment

Run your version of bisect on these functions. Write down the values of the computed root, and the number of steps the algorithm took. One of the functions may cause bisection to fail; if you notice this, please explain what happened. When you are done, prepare a short table with the function number, the computed root, and the number of steps, and mail it to me. You should finish this assignment, and mail me the results, before our next class on 22 September.


Last revised on 20 September 1999.