V.4.1 Introduction

In many simulations one needs to integrate a function over a specified interval. From standard calculus it is known that only a limited number of functions can be integrated analytically.  In this chapter we will study both numerical integration schemes and the so-called stochastic integration methods (Monte Carlo methods) and compare them, both in terms of number of operations and in terms of accuracy.

V.4.2 Numerical Approximations

Introduction

Many numerical methods exist to calculate integrals. In this subsection we will study some of them to integrate the function

You will compare the numerical methods, and, in the next chapter, also compare numerical integration schemes with a method based on random numbers (the hit and miss integration scheme).

REQUIRED
V.4.2a  Plot the function f[x] and calculate the indefinite integral of the function, and the definite integral in the interval (0,1).

The numerical integration methods studied in this section are based on interpolation of the integrand by polynomial functions between sample points. The integral is computed by summing up the exact integral for the approximate polynomial functions, which can be solved analytically. The main assumption is that the integrand is smooth.

Straight Summation

The simplest approximation of the integral is to sample the function on equidistant points over distance and sum up the area of the rectangles formed by the function value on the left sampling point and the sampling length . This computation uses the summation command with no refinements.

ADVANCED
V.4.2b Show graphically for example using a bar chart how this summation works (by e.g. plotting a function, and plotting the rectangles that are used in the approximation; only use 4 sample points).

The Mathematica implementation is as follows.

OPTIONAL
V.4.2c  Measure and plot the error in the straight summation as a function of .

Midpoint Summation

A very simple refinement of the previous approximation is to calculate the integral by using the function value in the middle of two sample points. This method is the so-called midpoint integration scheme.

OPTIONAL
V.4.2d Show graphically how this summation works.

The Mathematica implementation is as follows.

REQUIRED
V.4.2e  Measure and plot the error in the midpoint summation as a function of .

Trapezoidal Summation

The next refinement is the so-called trapezoidal method. In this method the integrand is approximated by a linear function between the sample points. The integral is now approximated by summing the areas of the resulting trapezia.

In Mathematica the trapezoidal rule looks as follows,

REQUIRED
V.4.2f  Measure and plot the error in the trapezoidal summation as a function of .

Simpson's Summation

A more sophisticated method is the Simpson's quadrature rule. In this method the integrand is approximated by a quadratic function between the sample points. It turns out that Simpson's rule can be expressed as,  

REQUIRED
V.4.2g  Measure the error in Simpson's as a function of .

ADVANCED
V.4.2h  Compare all schemes in terms of complexity and accuracy.

V.4.3 'Hit and Miss' Integration     

Introduction

Another way to approximate integrals is by using stochastic techniques. These techniques are especially useful to evaluate integrals of high dimensional functions. To illustrate the method we will integrate in a one-dimensional method. The idea of Hit and Miss integration is to randomly select points and count all points that are 'below' the function (the hits) and neglect the points that are 'above' the function (the misses). The ratio between the hits and the total number of trials gives an approximation of the integral.

The program evaluates definite integrals by creating a "box" with dimensions defined by the limits of integration and the maximum and minimum values of the function on the selected interval. The computer randomly generates the coordinates of a point within the box. The function is evaluated at the randomly selected value. If the random point falls within the area of the integral, it is counted by Mathematica. After performing the procedure for the specified number of points (iterations), the computer calculates the final result by multiplying the total area of the box by the percentage of points that fell within the bounds of the integral. Of course, many functions have regions that fall below the line . If the program encounters a point that is within the area of the integral but below the -axis, it subtracts one from the total number of points counted.
One would expect that the more points generated for each evaluation of the integral, the more accurate the answer will be. Since random numbers are used, however, there will always be variation from one determination to the next. Because of the lack of precision, the program is designed to perform many calculations with a varying number of iterations. Doing so makes it possible to "zero in" on the actual answer.

The Program

The following is a listing of the entire Hit and Miss program. The next section explains some of the details of the program.

Details of the program:

The example of using the function:

First parameter is a function, the limits of integration are specified in the next parameter.  The last parameter is the number of iterations.
The following lines determine the maximum and minimum value of the function on the specified interval.  The program uses these results as boundaries between which random values are selected.  The values are also necessary to determine the area of the "box."

The next part of the program generates random coordinates and determines how each point should be counted. If the random value is less than or equal to the function evaluated at the random value , it is counted.

The next lines add the remaining part of the integral.

For example consider integration of the function illustrated in the figure below. The bounding box boundaries are determined by =0.05 and =0.3. The results obtained using the code above have to be corrected by adding the size of the grey area ))

The code performing correction is presented below:

Example 1

In this example one can see the standard deviation of the estimated results decrease with number of points. This is done by creating the table of results of Hit and Miss function for increasing number of points (from 100 to 2000 with step 100). For each step, the algorithm is repeated 10 times and the results are used for calculating a standard deviation. (Standard deviation gives the idea about the error of measured value of the integral. It is calculated by comparing each of the results with their average value.)

Example 2

In this example we run the program for a function that can be easily integrated. Below, the actual value of the definite integral is plotted along with 'Hit and Miss' results.

The graph demonstrates that the accuracy and precision of the results change as the number of points used in each calculation increases. Generally, the more points that are used, the more accurate and precise the output will be.

REQUIRED
V.4.3a Investigate the error in the approximation of the following integral with the Hit and Miss method: , for different number of trials and for λ = 0.1, 1, 5.

ADVANCED
V.4.3b Investigate the error in the approximation of the following integral with the Simpson numerical integration and the Hit and Miss method: . Compare both methods.

Stochastic integration is very well suited to integrate high dimensional functions. Now you can understand why, by investigating how the Hit and Miss method depends on the dimensionality of the Integrand, and by comparing this with numerical methods (such as Simpson's rule).