Numerical Integration

Numerical integration is the process of numerically computing a definite integral. There are many occasions where numerical integration is important. Examples include calculating the area of a shape or the area under a graph, and solving differential equations.

For the purpose of this discussion we'll calculate the area of half a circle of radius 1. We already know this area to be re/2. So in a sense, calculating the area of half a circle is equivalent to calculating the numerical value of re.

First, we create two vectors: x and y. These two vectors satisfy the circle equation x2 + y2 = 1:

>>> x = linspace(-l, 1, N) >>> y = sqrt(l-x**2)

I chose the variable N arbitrarily; N is the number of points in the vectors x and y.

To visualize the numerical integration, I plot rectangles that approximate the area of the circle:

... rect = Rectangle((x[i], o), dx, 0.5*(y[i]+y[i+l])) ... gca().add_patch(rect)

>>> title('Approximating the area of half a circle') >>> axis('equal')

The area under the curve, that is, the integral, is approximately the sum of these squares. Each square's area is 0.5*(y[i]+y[i+1])*dx, so the total sum can be written as follows:

>>> dx*(sum(y[0:-l]+y[l:])) 2.9175533787759904

I've multiplied the result by 2 so we can compare with re instead of re /2. Obviously, the bigger N is, the closer this number will be to re:

... print "N=%d, estimated pi is %f" % (N, est_pi)

N=5, estimated pi is 2.732051 N=10, estimated pi is 3.019232 N=20, estimated pi is 3.101560 N=100, estimated pi is 3.138218

As you can see, for N = 100, the accuracy is about 1 percent. Figure 8-1 captures this visually.

Figure 8-1. Calculating the area of a circle

In calculating the area of the circle, I chose values that are evenly spaced. In case you'd like to use non-evenly spaced values, the implementation is more complex. Also, the method uses rectangles to approximate the area under the curve, but in this particular example (and many others), trapezoidals are probably better suited, which brings us to the function trapz(y, x). The function accepts vectors y and x and returns the numerical integral. The following performs numerical integration of non-evenly spaced x values using the function trapz():

>>> x = array([-l, -0.9, -0.4, 0.0, 0.4, 0.9, l]) >>> y = sqrt(l-x**2) >>> trapz(y, x)*2 2.9727951234089831

Figure 8-2 shows a visual representation of the trapezoidal integration.

Trapezoidal integration

Figure 8-2. Calculating the area of a circle using the trapezoidal method and non-evenly spaced values

Was this article helpful?

0 0

Post a comment