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.
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? Do not spend time trying to make this function fancy.)
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 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," then 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.
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. There are two reasons for writing it as a function instead of a script m-file (without the function header) or by typing everything at the command line. The first reason is a practical one: we will be executing the algorithm several times, with differing end points, and functions allow us to use temporary variables without worrying about whether or not they have already been used before. The second reason is pedantic: you should learn how to do this now because it is an important tool in your toolbox.
Atkinson describes the bisection algorithm (Section 2.1) in the following manner.
Bisect(f,a,b,root,epsilon)
1. Define c:=(a+b)/2.
2. If b-c < epsilon, then accept root:=c, and exit.
3. If sign(f(b)) * sign(f(c))£0, then a:=c; otherwise, b:=c.
4. Return to step 1.
This language of this description is based on a very old computer language called Algol. One objective of this lab is to show how to write it in the Matlab language. As a starting point, let's fix cosy as the function and leave it out of the calling sequence. Let's also fix epsilon=0.000001. Note also that, in Matlab, returned values are separated from arguments, so the variable root, which I will call approx_root to make a point, is written on the left of an equal sign. The first line of the file bisect.m should be:
function approx_root = bisect( a, b)
The second line should be a comment explaining what the function does.
Now think about the iteration. Atkinson's algorithm employs the "return to Step 1" device to indicate uncounted repetition of Steps 1-4. Matlab provides a more elegant construct for uncounted repetition: the while loop. In Atkinson's Step 2, the algorithm is terminated. Matlab while loops usually terminate at their end, so we continue iterating for as long as (b-c) is larger than epsilon, rather than stopping when (b-c) is smaller. This reversal of the test is sometimes confusing, but it is how it is done.
Here is 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.
function approx_root = bisect ( a, b )
% bisect finds an approximate root of the function cosy using bisection
fa = cosy(a);
fb = cosy(b);
while ( abs ( b - a ) > 0.000001 )
c = ( a + b ) / 2;
approx_root = c;
fc = cosy(c);
%
% What follows is just a nice way to print out a little table. It does not add to the
% algorithm itself, it only makes it easier to see what is going on at runtime.
%
[ a, c, b;
fa, fc, fb ]
if ( sign(fb) * sign(fc) <= 0 )
a = c;
fa = fc;
else
b = c;
fb = fc;
end
end
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?
Review:
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)?
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 )
is acceptable to MATLAB and produces the same root as before.
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 |
INTERVAL |
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.
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. Complete the table sketched above.
Back to the MATH2070 page.
Last revised on 20 September 1999.