Creative Commons License Copyright © Michael Richmond. This work is licensed under a Creative Commons License.

Debugging: How to find the errors ....

As some of you may have discovered, there are times when -- despite your carefully crafted pseudocode -- you make a mistake in translating your algorithm into actual Scilab code. We call these mistakes bugs. The process of removing them from your program is called debugging.

You can find many, many references and guides to debugging all over the Net; going to Google and typing debugging techniques returns about 1,340,000 items. There are good chapters on debugging in some of the references for the course.

Let me show you just one, very simple technique that I use frequently. It isn't particular to Scilab, but can be used in just about any programming language.


A buggy program

Consider the following program, which tries to find the root via the bisection technique, but fails.
function root = root_bisect(lower_bound, upper_bound, epsilon, func)
//
// Find the root of the given function within the given interval,
// using the bisection method.
//
// Arguments:
//
//      lower_bound     (input)    lower bound of interval
//                                     on which to search for root
//
//      upper_bound     (input)    upper bound of interval
//
//      epsilon         (input)    when fractional change in root
//                                     reaches this limit, stop
//
//      func            (input)    name of function whose root
//                                     we seek -- should take
//                                     a single argument
//
//      root           (output)    argument at which the given
//                                     function has value zero
//
// MWR 3/25/2002
//

// sanity checks
if (lower_bound >= upper_bound)
   error('lower_bound must be less than upper_bound');
end
if (epsilon <= 0)
   error('epsilon must be greater than zero');
end

// if there is no guaranteed root between the given bounds, 
// then quit with an error message
y_low = feval(lower_bound, func);
y_high = feval(upper_bound, func);
if (sign(y_low) == sign(y_high))
   error('function has same sign at both bounds -- quitting');
end

done = 0;
iterations = 0;
old_mid = lower_bound;
while (done == 0) 
   y_low = feval(lower_bound, func);
   y_high = feval(upper_bound, func);
   
   // Pick the half of the interval in which the root must exist
   // 'mid' is going to be next root candidate
   mid = (upper_bound - lower_bound)/2.0;
   y_mid = feval(mid, func);
   if (sign(y_mid) == sign(y_low))
      // root is in upper half
      lower_bound = mid;
   else
      // root is in lower half
      upper_bound = mid;
   end
   
   // now, check to see if we are close enough
   if (mid == 0)
      // we can't calculate fractional change here, so use criterion
      // that absolute value of the change in candidates is less than epsilon
      if (abs(mid - old_mid) < epsilon)
         done = 1;
      end
   else
      frac_change = abs((mid - old_mid)/mid);
      if (frac_change < epsilon)
         done = 1;
      end
   end
   
   // we need to save a copy of the current root candidate, so we
   // can calculate the fractional change since the previous one
   old_mid = mid;
   
   iterations = iterations + 1;
end

// when we get here, the latest root candidate is in 'mid'
root = mid;

endfunction

Now, if we run the program with a very simple function funcz.sci


function y = funcz(x)
y = x*x - 5;

on the interval [0, 5], we ought to find the root sqrt(5) = 2.236. But when we try, what happens?

Obviously, something is wrong. But what? And where in the source code? Can you find the error? It's pretty hard when all you have to use is the source code itself.


Print and pause to find the bug

Here's a second version of exactly the same program -- it still contains the error. However, this version includes some extra code which is designed to track it down.

function root = root_bisect(lower_bound, upper_bound, epsilon, func)
//
// Find the root of the given function within the given interval,
// using the bisection method.
//
// Arguments:
//
//      lower_bound     (input)    lower bound of interval
//                                     on which to search for root
//
//      upper_bound     (input)    upper bound of interval
//
//      epsilon         (input)    when fractional change in root
//                                     reaches this limit, stop
//
//      func            (input)    name of function whose root
//                                     we seek -- should take
//                                     a single argument
//
//      root           (output)    argument at which the given
//                                     function has value zero
//
// MWR 3/25/2002
//

// set this to 1 to watch diagnostic messages as the function runs
verbose = 1;

// sanity checks
if (lower_bound >= upper_bound)
   error('lower_bound must be less than upper_bound');
end
if (epsilon <= 0)
   error('epsilon must be greater than zero');
end

// if there is no guaranteed root between the given bounds, 
// then quit with an error message
y_low = feval(lower_bound, func);
y_high = feval(upper_bound, func);
if (sign(y_low) == sign(y_high))
   error('function has same sign at both bounds -- quitting');
end

done = 0;
iterations = 0;
old_mid = lower_bound;
while (done == 0) 
   y_low = feval(lower_bound, func);
   y_high = feval(upper_bound, func);
   
   if (verbose > 0)
      mprintf('iter %5d   %11.7e %9.4e   %11.7e %9.4e ', ...
         iterations, lower_bound, y_low, upper_bound, y_high);
      // pause;
   end
   
   // Pick the half of the interval in which the root must exist
   // 'mid' is going to be next root candidate
   mid = (upper_bound - lower_bound)/2.0;
   y_mid = feval(mid, func);
   if (sign(y_mid) == sign(y_low))
      // root is in upper half
      lower_bound = mid;
   else
      // root is in lower half
      upper_bound = mid;
   end
   
   // now, check to see if we are close enough
   if (mid == 0)
      if (verbose > 0) 
         mprintf('  abs change %9.4e ... ', abs(mid - old_mid));
      end
      // we can't calculate fractional change here, so use criterion
      // that absolute value of the change in candidates is less than epsilon
      if (abs(mid - old_mid) < epsilon)
         if (verbose > 0)
            mprintf(' stop \n');
         end
         done = 1;
      else
         if (verbose > 0)
            mprintf(' keep going \n');
         end
      end
   else
      frac_change = abs((mid - old_mid)/mid);
      if (verbose > 0) 
         mprintf('  %9.4e \n', frac_change);
      end
      if (frac_change < epsilon)
         if (verbose > 0)
            mprintf(' stop \n');
         end
         done = 1;
      else
         if (verbose > 0)
            mprintf(' keep going \n');
         end
      end
   end
   
   // we need to save a copy of the current root candidate, so we
   // can calculate the fractional change since the previous one
   old_mid = mid;
   
   iterations = iterations + 1;
end

// when we get here, the latest root candidate is in 'mid'
root = mid;


endfunction

Look near the top of the main loop. There's a new bit of code that looks like this:

   
   if (verbose > 0)
      mprintf('iter %5d   %11.7e %9.4e   %11.7e %9.4e ', ...
         iterations, lower_bound, y_low, upper_bound, y_high);
      pause;
   end

What does it do?

Please copy this program to your local disk, rename it to root_bisect.sci, and run it. Watch what happens.

Now that you've seen the output of the program, can you explain why it isn't finding the root properly? Can you now find the place in the program which is most likely to contain the error?


Leaving debugging mode

It is often useful to include extra mprintf statements in your program as you are developing it. Sometimes, you may add a LOT of these debugging "print" lines. What happens when it's time to turn the program in?

Try editing your version of the program: change the verbose variable to zero, and run it again. Now what happens?