CPSC/Math 402 Lab

Numerical Differentiation

The program derivative.cc uses three different formulas for approximating the derivative. Three functions (sin(x), log(x), and 1 + x + x3) are built into the program (note: the third function wasn't in the version handed out in class). The three formulas are:

The program lets you choose the function to use, the value of x, the starting value for h and the number of iterations. At each iteration h is divided by some number for the next iteration - you also enter that number (typically you would want to divide by 2 or 4 or 10). The program prints tables of results for each approximation forumla. The tables contain the iteration number, the estimate of f' given by the formula and the actual error in the estimate (which in this case can be computed because we actually know what f' is and don't really need to do numerical differentiation!).

Save the program to your directory either by clicking on the link or using the Unix copy command to copy to your current directory: 

      cp ~cpsc/public_html/Spring2006/CPSC402A/derivative.cc  .
The command to compile is: 
      g++ derivative.cc
or
      g++ derivative.cc  -o diff
The first command puts the executable in a file named a.out while the second puts it in a file named diff (you can substitute any name here).

To run the program type: 

        ./a.out
or
        ./diff
depending on what you named the executable.

Run the program as directed below and answer the associated questions. You should get a printout (use script) of the output that you are examining for each question so you can explicitly show what you are talking about.

  1. Use the function sin(x) and estimate the derivative at x = 0.5. Run the program for 30 iterations (at least). Start h at 1 and divide by 4 at each iteration.
    1. Since the derivative is a limit as h approaches 0, then in theory as h gets smaller with each iteration the estimate of the derivative should get closer to the actual derivative. In fact, the forward difference formula is an O(h) estimate so since h is being divided by 4 at each iterate one would expect the error to be reduced by about a factor of 4 (taking not quite two iterates to gain a digit of accuracy). So what actually happens? Examine the sequence of iterates for forward difference and determine whether or not the theoretical results hold at any point in the sequence (it should be obvious that they do not hold throughout the sequence).
    2. Both of the other formulas are not O(h), they are O(h2). So, at each iterate as h is divided by 4 what would you expect the error to be divided by?
    3. Is there a part of the sequence of iterations where the error in the second order algorithms matches the theoretical results? Explain.
    4. Clearly things go awry if one gets too carried away with iteration. The theoretical error mentioned above is only the truncation error; that is, the error in using an approximation formula (the difference between the actual value and the value the approximation formula would get if the approximation could be accurately computed). Roundoff error rears its ugly head during the calculation of the approximation so the error is compounded. For each of the three formulas, when does roundoff error seriously kick in and start making matters worse? Which iteration and what is h at that iteration (note h = 4-n for some value of n at each iteration)? (NOTE: A theoretical analysis of truncation versus roundoff for the central difference formula is derived on pages 443 and 444 of the text. The graph show the truncation error decreasing as h decreases - goes right to left in the graph - while the round off error increases. There is an optimal point where the overall error is minimized. That should be evident from looking at the output of the program.)
    5. Explain why the roundoff error gets so bad. What is happening in the calculation of the approximation formulas? Why does the approximation eventually become 0 if you make h small enough? (By the way, if you did a lot of iterations - say 540 - the approximation would eventually become NaN, Not a Number - why?)
    6. In theory, the error in using the central difference approximation is about half that in using the other second order approximation. Is this evident in practice? Explain!

  2. Do the same as above with the function log(x) and estimate the derivative at x = 1. Answer the same questions for parts (a) through (d) above.

  3. Use the program to do #13 on page 446. To answer part (c) think about what the error formula is for the forward difference approximation and relate that to your observed error for the two different values of x0.

  4. Add code to the program to approximate the derivative using the formula in problem 4 on page 445. Run the program to estimate the order of the approximation. For homework do parts (a) and (b) - that is, derive the formula and the error term. Be sure to comment on whether or not your theoretical estimate of order was verified numerically by the program - explain!!! (Don't just say yes it was or no it wasn't!)

  5. Suggest at least one stopping criterion so the program would stop when the iterates are close to the correct answer. NOTE: Don't base your stopping criteria on knowing the exact error!! Assume you don't know the actual derivative so you don't know the exact error. All you know is the function, the iterates, and h - some combination of those is what you need to base the stopping criterion on.

Completed work is due Wednesday, April 12.