LAB #3: Roots of Equations

TABLE OF CONTENTS

• Introduction
• A Sample Problem
• The Bisection Idea
• Programming Bisection
• A Bisection Program
• Variable Function Names
• Some Nonlinear Functions
• ASSIGNMENT

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

abs ( f(x) )

A Sample Problem

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

cos ( x ) = x

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

        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 ] )

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

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 )

cos ( x ) = x?

Discussion:

• What does the sign(x) function do? What if x is 0?
• What does the statement [ a, c, b; fa, fc, fb ] do? What other ways could this be done?
• Under what circumstances could this program return an incorrect result? That is, APPROX_ROOT is not close to a root of the function.
• Under what circumstances could this program run "forever?" How do you stop MATLAB when it gets into an infinite loop?
• If I compute x = bisect ( a, b), then what can I claim is true about x and about f(x)?

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:

• In the "declaration line" for bisect, include a dummy function name. We might as well call it f, so the line becomes

  function approx_root = bisect ( f, a, b )

• Where we evaluated the function by writing f(x), we now have to use the MATLAB "helper" function called feval. For instance, the line

  fa = cosy(a);

  must be rewritten as:

  fa = feval ( f, a );

• When we use bisect, we must include the function name, in quotes:

  x = bisect ( 'cosy', 0, 3 )

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 )

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.