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.
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
on the interval [0, 5], we ought to find the root sqrt(5) = 2.236. But when we try, what happens?function y = funcz(x) y = x*x - 5;
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.
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?
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?
you can switch the program from its "debugging" state to its "production" state, or back.// set this to 1 to watch diagnostic messages as the function runs verbose = 1;
Try editing your version of the program: change the verbose variable to zero, and run it again. Now what happens?