C Course Homework

Homework 4 - Integration

Things You Will Need to Know/Learn:

• Functions (Chapter 4)
• Using C's Random Number Generator (see Hw4 hints )
• Integration
  • Mini-refresher on integrals.
  • Using the Simpson Rule to estimate an integral.
  • Using the Monte Carlo Method to estimate an integral.

Program Description

Simpson Rule

The first technique will use the Simpson Rule to evaluate an integral. The Simpson rule uses rectangles to approximate the area under a curve. Use the following prototype for this function.

     double integrate_rect(double a, double b, double F(double), int n);
      /* The interval (a,b) is partitioned into n parts.  The function    */
      /*  F(x) is computed on the midpoint of each partition to calculate */
      /*  the area of each rectangle.  The areas of all the rectangles    */
      /*  are then added to approximate the entire area under the curve.  */

Monte Carlo Method

The second technique will use a random number generator to do a Monte Carlo estimation of the area under a curve. The theory of this method will be explained in class. Use the following prototype for this function.

     double integrate_mc(double a, double b, double max, double F(double), int n);
      /* The interval of integration is (a,b) and max is a float larger   */
      /*  than any function value in this interval.  The dimensions       */
      /*  (b - a) * max define a rectangle within which F(x) is drawn.    */
      /*  n points are sampled at random using a random number generator  */
      /*  to produce (x,y).  If y <= F(x), then we count this as a     */
      /*  "yes", but if y > F(x) then we count this as a "no".   	  */
      /* The number of "yesses"/n * (b - a) * max is the area estimate.   */

Testing Both Methods

In addition to implementing the two routines above, you will have to write the function F(x) to test both methods. This function will be passed as an argument to the two routines you write to implement each method. An easy function to use for initial testing is F(x) = x in the interval (0,1). Try not to exceed n = 1000, as this will take a long time to run. Here's the function definition for F(x) = x:

     double F(double x)
     {
        return x;
     }

Main

Your main routine should set the interval you wish to check and the maximum value which is larger than F(x) in the interval (for the Monte Carlo method). You can hard-code (set) the interval and the maximum to values--you don't need to do any input from the user. You should then run your functions for five different values of n (4, 16, 64, 256, 1024). A for loop works nicely here. You should print the results from each iteration as well as a difference between the two methods.

Take a look at the sample output .

Grading

Style -- 2 points total

• proper comments and variable names 1 point
  • Name, class, date, instructor, program description, and useful in-line comments.
    Variable names which hint towards the use of the variable.
  • proper indentation, use of whitespace, and overall appearance 1 point

Correctness -- 6 points total

• correct main() 1 point
  • For 5 different n values (4, 16, 64, 256, and 1024), the main program should make a call to the two integration methods (integrate_rect and integrate_mc), print out each area estimate, and print out the difference between the two estimates. The main program should not accept any user input.
• correct implementation of integrate_rect 2 points
  
  • Part 1 - overall (1 point): The function should have the prototype:
    double integrate_rect(double a, double b, double F(double), int n);
    The student should not have changed the prototype.
    The function should contain an "area-calculated-so-far" variable initialized to 0.
    The function should contain a loop which is executed n times to sum up the areas of all the rectangles.
    The function should return the total area of all the summed up rectangles.
    
  • Part 2 - within the loop (1 point):
    For each rectangle, the area should be computed using the midpoint to calculate the rectangle's height, and the width of the entire interval divided by n to calculate the rectangle's width.
• correct implementation of integrate_mc 2 points
  • Part 1 - overall (1 point):
    The function should have the prototype:
    double integrate_mc(double a, double b, double max, double F(double), int n);
    The student should not have changed the prototype.
    The function should contain a integer counter variable initialized to 0.
    The function should contain a loop which is executed n times to count how many random points out of n fall below the function F. The function should return the total area of the bounding box
    times the counter value divided by n.
  • Part 2 - within the loop (1 point):
    For each iteration, a point has to be chosen at random. Make sure the random x-coordinate will be a real number in the range a to b, and make sure the random y-coordinate will be a real number in the range 0 to max. Also be sure that the random x-coordinate is chosen independently of the random y-coordinate. If the random point is under the curve (if randomy <= F(randomx)), the counter should be incremented.
• correct implementation of function F 1 point
  • The prototype should be: double F(double); The function should be passed to the two integration methods and used in each of their function bodies to perform the two area
    approximations.
• Program execution -- 2 points total
  • Make a copy of the program. Edit your copy and change the following hard-coded variables:
    F = x
    a = 1.0
    b = 3.0
    max = 3.0
    • If the student didn't make use of #defines, make a style comment to them but do not deduct points. Compile your copy of the program.
    • Run "% a.out": The output should contain area approximations shown in HW4.output1 with some margin of error for the monte carlo approximation depending on how the student decided to seed his or her program. Note that the simpson area, the monte-carlo area, and the difference between the two are the only values which are required to be printed to the screen.
      The area estimates should get closer and closer to the real area as n gets bigger. The real area for this function and interval is 4.0.
      
      
  • Do the same as above except change the hard-coded variables to:
    F = 3*x*x - 2*x
    a = 2.5
    b = 4.0
    max = 50.0
    
    
  • The real area for this function and interval is 38.625. The output should be similar to HW4.output2
    
    
  • Do the same as above except change the hard-coded variables to:
    F = exp(x) (needs #include <math.h>)
    a = 0.3
    b = 3.14159
    max = 27.0
    
    The real area for this function and interval is 21.80026202...
    The output should be similar to HW4.output3

Deductions

Deduct 3 points if the program does not compile.

(use "gcc <filename.c>" or "gcc <filename.c> -lm")

• Make sure you understand the background behind this assignment. The goal is to find the area under a curve defined by some function F(x). The Simpson Rule and the Monte Carlo Method are two techniques for arriving at an estimate for the area. Finding this area is the same as integrating using calculus.

• For help with passing a function as an argument, here's a sample program , and here's the corresponding output .

• For help with using rand(), check out the mini-program from hw3. You will have use the random number generator for implementing the Monte Carlo method. Note that you will be using the value returned by rand() differently than in hw3.
  

              rand()/((double)(LONG_MAX))
         

Before dividing rand() by LONG_MAX we have to cast LONG_MAX to a double so that the division is a floating point division and not an integer division. (Recall that an integer division will truncate any fractional part, so in this case if we did not do the cast we would get either 0 or 1, but nothing in between. Casting is mentioned on p. 240 of the book.)
Now that you have a value between 0.0 and 1.0, you can manipulate this value with a bit of math to pick a point in the desired interval.

FAQ

• Hard-coding F into the program
• What functions to test and output
• Random limit on CATS machines
• Hard-coding revisited
• Printing out your function, F
• Monte Carlo method and probability

• here's my version of hw4.