CS 2073 -- Program 6

CS 2073, Spring 2006
Program 6
Numerical Integration
    Week 7: Feb 27-Mar 3
Due (on time): 2006-03-10  23:59:59
Due (late):        2006-03-20  23:59:59

Program 6 must be emailed to: nrwagner@cs.utsa.edu
following directions for: running and submitting a C program, with deadlines:

2006-03-10 23:59:59 (that's Friday, 10 March 2006, 11:59:59 pm) for full credit.
2006-03-20 23:59:59 (that's Monday, 20 March 2006, 11:59:59 pm) for 75% credit.

Introduction: For this assignment, we want a program that will do numerical integration. You don't really need to know any calculus, since for us the integral of a function will just be the area under its graph, or its average value. This assignment will look at four numerical integration methods:

the Rectangle method (worked out below),
the Trapezoid method,
Simpson's method, and
a Monte-Carlo method.

I will give a solution to the first method (Rectangle method), and ask you to program the remaining three. I want you to practice using functions, so you must use actual C functions for the functions being integrated. (It makes the formulas look simpler anyway.)

Here are a few links giving write-ups about these methods:

See Random Numbers for information about generating the random number needed for this program.
See Comparisons: while and for loops for examples of loops that add sequences of numbers.

Details of the four methods: We will be finding the value of the integral of a function f(x), for x from a to b. You will also start with an integer n representing the number of intervals to divide the segment from a to b into. The size of each interval is h = (b - a)/n. Given these starting values, here are the four methods:

One version of the rectangle method uses the formula:

Rectangle-Integral = (approximately)
     (h)*[f(a) + f(a+h) + f(a+2*h) + f(a+3*h) + f(a+4*h) + . . . + f(a+i*h) +
                             . . .  + f(a+(n-2)*h) + f(a+(n-1)*h)]

The trapezoid method uses the formula:

Trapezoid-Integral = (approximately)
     (h/2)*[f(a) + 2*f(a+h) + 2*f(a+2*h) + 2*f(a+3*h) + 2*f(a+4*h) + . . . + 2*f(a+i*h) +
                             . . .  + 2*f(a+(n-2)*h) + 2*f(a+(n-1)*h) + f(b)]

Similarly, Simpson's method uses the following formula, where we want n to be even, so that the alternating weights of 1, 4, 2, 4 will come out as shown: ending with 2, 4, 1.

Simpson's-Integral = (approximately)
     (h/3)*[f(a) + 4*f(a+h) + 2*f(a+2*h) + 4*f(a+3*h) + 2*f(a+4*h) + . . .
                                     + 2*f(a+(n-2)*h) + 4*f(a+(n-1)*h) + f(b)]

Finally, a Monte-Carlo method might use:
Here the numbers x₁, x₂, . . . , x_n are randomly chosen from the interval from a to b.

Details about the Rectangle method:

For example, suppose f(x) = x² - x + 1, a = 1 to b = 3, and n = 4. Then h = (b-a)/n = (3-1)/4 = 1/2. The sum given by the rectangle method is:

(1/2)*(f(1) + f(1.5) + f(2) + f(2.5)) =

(1/2)*(1 + 1.75 + 3 + 4.75) = 10.5/2 = 5.25

The exact answer is the integral of f(x) from 1 to 3, which works out to be 20/3 = 6.6666...

Now, lets write a little loop to calculate this. I'll leave n as a variable, so we can put in anything we want for it. Assuming that the function f has been defined correctly, on the left below is what the main part could look like. On the right are 4 runs, starting with 4 different values for n:

   int n = 4;
   double a = 1;
   double b = 3;
   double h = (b - a)/n;
   int i;
   double sum = 0;
   printf("n: %i, a: %.2f, b: %.2f, h: %.4f\n",  
      n, a, b, h);
   
   for (i = 0; i < n; i++)
      sum = sum + h*f(a + i*h);
   printf("sum: %.6f\n", sum);

n: 4, a: 1.00, b: 3.00, h: 0.5000
sum: 5.250000
------------------------------------
n: 10, a: 1.00, b: 3.00, h: 0.2000
sum: 6.080000
------------------------------------
n: 100, a: 1.00, b: 3.00, h: 0.0200
sum: 6.606800
------------------------------------
n: 10000, a: 1.00, b: 3.00, h: 0.0002
sum: 6.666067

Integrals to approximate numerically: You should use the specific function f below, the specific values 0 and 1, and find its approximate integral:

f(x) = 1/(1 + x²), for x from 0 to 1, or:

In the case of f we can determine the exact answer analytically: The indefinite integral is just arctan(x), which is 0 when x = 0, so the exact answer to this part is arctan(1), which is pi/4 = 0.785398163397448. For this function, Simpson's rule is very accurate, while the Trapezoid method is intermediate, and we always expect a Monte Carlo method to be the least accurate.

After determining the integral for f, you should use the same code to find numerical integral approximations for the function g below, with the limits -1 and 1. This function is used in statistics. There is no formula for the exact answer in this case.
g(x) = e^-x², for x from -1 to 1, or:

More details: So write a program that computes the first function f, using each of the three integration methods, and using n = 10, n = 1000 and n = 100000. (There should be 9 answers altogether.)

After that you should change your program to compute the second function g and produce another 9 answers.

The output should be clearly labeled with the function (f or g above), the values of a and b, the value of n, and the integration method.

Hints and suggestions:

Multipliers on terms in the sum: Don't panic when computing the sum above. You have terms for values of i from 0 to n inclusive. With the Trapezoid rule each term is multiplied by 2 except for the 0th and the nth terms. (Notice that there are n+1 terms. The example with the Rectangle method had only n terms.)
With Simpson's rule, the 0th and the nth terms are multiplied by 1. Otherwise, odd-numbered terms are multiplied by 4 and even-numbered terms are multiplied by 2. (Simpson's rule requires that n be an even number.). To check that term number i is odd, you can check:
The fourth method: Monte-Carlo: The Monte-Carlo method is even easier. Just add up n terms of the form f(random point) and then multiply by h = (b - a)/n. (Notice that in the two previous methods, there are n+1 terms, but here there are only n terms.)
The first program at the link above (Random Numbers) has a function rand_double2 that returns a random number between two input parameters x and y. For the different values of n, you generate n values of x between a and b, and then you produce an approximation for the area under the curve, which is the average value of the numbers f(x) (average height) times b - a (the width).
Getting n to take on the proper values: In order to get n to take on the three successive values 10, 1000, and 100000, you can use the for loop:

Extra credit part: For extra credit, try to find an approximate value for the following. After you have found this approximate value, square it and see if you recognize this number. (Hint: the function g(x) gets extremely small as abs(x) increases.)

What you should email: Refer to the submissions directions and to deadlines at the top of this page. The text file that you submit should first have Your Name, the Course Number, and the Program Number. The rest of the file should have the following in it, in the order below, and clearly labeled, including at the beginning the appropriate item letters: a, b, c, etc.

Contents of email submission for Program 6:
Last Name, First Name; Course Number; Program Number.

C source for your program, integral.c. (Or whatever you wish to name it.)
Results of a run of the program used to determine the integral of f three ways, each for three values of n.
Results of another run doing the same as b above for the function g.
(Extra credit) Results of a run computing the approximate integral of g from minus infinity to plus infinity. Then the square of this answer.

Revision date: 2006-03-01. (Please use ISO 8601, the International Standard.)