A nontrivial example

Notice from Figure 5 that has a root between 1 and 2 (of course, this root is , but we feign ignorance for a moment). A general algorithm for nonlinear root-finding is the method of bisection, which takes a function and an interval on which function changes sign, and repeatedly bisects the interval until the root is trapped in a very small interval.

A function implementing the method of bisection illustrates many of the important techniques of programming in Matlab. The first important technique, without which a useful bisection routine cannot be written, is the ability to pass the name of one function to another function. In this case, bisect needs to know the name of the function whose root it is to find. This name can be passed as a string (the alternative is to ``hard-code'' the name in bisect.m, which means that each time one wants to use bisect with a different function, the file bisect.m must be modified. This style of programming is to be avoided.).

The built-in function feval is needed to evaluate a function whose name is known (as a string). Thus, interactively

>> fcn(2)
ans =
   -0.7568

and

>> feval('fcn',2)
ans =
   -0.7568

are equivalent (notice that single quotes are used to delimit a string). A variable can also be assigned the value of a string:

>> str = 'fcn'
str =
f
>> feval(str,2)
ans =
   -0.7568

See help strings for information on how Matlab handles strings.

The following Matlab program uses the string facility to pass the name of a function to bisect. A % sign indicates that the rest of the line is a comment.

function c = bisect(fn,a,b,tol)

% c = bisect('fn',a,b,tol)
%
%   This function locates a root of the function fn on the interval
%   [a,b] to within a tolerance of tol.  It is assumed that the function
%   has opposite signs at a and b.

% Evaluate the function at the endpoints and check to see if it
% changes sign.

fa = feval(fn,a);
fb = feval(fn,b);

if fa*fb >= 0
   error('The function must have opposite signs at a and b')
end

% The flag done is used to flag the unlikely event that we find
% the root exactly before the interval has been sufficiently reduced.

done = 0;

% Bisect the interval

c = (a+b)/2;

% Main loop

while abs(a-b) > 2*tol & ~done

   % Evaluate the function at the midpoint

   fc = feval(fn,c);

   if fa*fc < 0       % The root is to the left of c
      b = c;
      fb = fc;
      c = (a+b)/2;
   elseif fc*fb < 0   % The root is to the right of c
      a = c;
      fa = fc;
      c = (a+b)/2;
   else               % We landed on the root
      done = 1;
   end

end

Assuming that this file is named bisect.m, it can be run as follows:

>> x = bisect('fcn',1,2,1e-6)
x =
    1.7725
>> sqrt(pi)-x
ans =
  -4.1087e-07

Not only can new Matlab commands be created with m-files, but the help system can be automatically extended. The help command will print the first comment block from an m-file:

>> help bisect

  c = bisect('fn',a,b,tol)
 
    This function locates a root of the function fn on the interval
    [a,b] to within a tolerance of tol.  It is assumed that the function
    has opposite signs at a and b.

(Something that may be confusing is the use of both fn and 'fn' in bisect.m. I put quotes around fn in the comment block to remind the user that a string must be passed. However, the variable fn is a string variable and does not need quotes in any command line.)

Notice the use of the error function near the beginning of the program. This function displays the string passed to it and exits the m-file.

At the risk of repeating myself, I want to re-emphasize a potentially troublesome point. In order to execute an m-file, Matlab must be able to find it, which means that it must be found in a directory in Matlab's path. The current working directory is always on the path; to display or change the path, use the path command. To display or change the working directory, use the cd command. As usual, help will provide more information.