CSC270: Numerical Integration




Fundamental Theorem of Calculus

b / | f(x) dx = F(b) - F(a) a / where F(x) is the `antiderivative': d -- F(x) = f(x) dx To compute integrals this way, we need to derive F(x) analytically. For example, if f(x) = x^2, the antiderivative is 1/3 x^3. Consider f(x) = sqrt( 1 + (cos x)^2 ). This is the length of the curve (x, sin x) for x in [0,pi]. There is no analytic expression for the antiderivative of f(x)! Thus, for some functions, we cannot use the fundamental theorem of calculus and we must resort to numerical integration.

Numerical Integration

Recall that the definite integral b / | f(x) dx a / is equal to the area under the curve f(x) in the interval [a,b]. We split the interval [a,b] into n smaller intervals, each of width delta_x: a b |--------------------------| x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 |--|--|--|--|--|--|--|--|--| nine intervals of width delta_x \ / Here, n = 9 \/ width delta_x Then delta_x = (b-a) / n and x_i = a + i * delta_x and the integral on [a,b] is just the sum over all the small intervals: x_{i+1} b / n-1 / | f(x) dx = sum( | f(x) dx ) a / i=0 / x_i n-1 = sum( A_i ) i=0 where A_i is the area under the curve in the interval [ x_i, x_{i+1} ].

We will estimate A_i in each of the intervals and sum up the estimates. That will be our estimate of the integral. Two things affect the quality of our integral: - the interval size: smaller is better - the method of estimating A_i

Riemann Sums

The first method of estimating A_i uses a constant function p(x) = c_i to approximate f(x). The constant function interpolates f(x) in the middle of the interval. Thus, c_i is chosen to be: x_i + x_{i+1} c_i = f( ------------- ) 2 and the area under p(x) approximates the area A_i under f(x): A_i ~ delta_x * c_i ~

Trapezoid Rule

The second method of estimating A_i uses a linear function p(x) = b_i x + c_i to approximate f(x) in the interval [ x_i, x_{i+1} ]. The linear function interpolates f(x) at the endpoints of the interval: p(x_i) = f(x_i) and p(x_{i+1}) = f(x_{i+1}) and the area under p(x) approximates the area A_i under f(x): f(x_i) + f(x_{i+1}) A_i ~ delta_x * ------------------- ~ 2

Simpson's Rule

The third method of estimating A_i uses a quadratic function p(x) = a_i x^2 + b_i x + c_i to approximate f(x) in the interval [ x_i, x_{i+1} ]. Here, the quadratic function interpolates f(x) at three points: p(x_i) = f(x_i) and p(x_{i+1}) = f(x_{i+1}) and p(x_{i+2}) = f(x_{i+2}) The area under p(x), which approximates A_i plus A_{i+1}, is 1 A_i + A_{i+1} ~ delta_x * --- * ( f(x_i) + 4 f(x_{i+1}) + f(x_{i+2}) ) ~ 3

Summary

In each case, we split the interval [a,b] into many smaller intervals of equal width. We estimate the integral of f(x) on [a,b] with the sum of our estimates for the areas A_i of the smaller intervals. For example, if we use Simpson's rule and split the interval [a,b] into n smaller intervals, each of size (b-a)/n, we get b / n-1 | f(x) dx = sum( A_i ) a / i=0 n/2 = sum( A_{2k-2} + A_{2k-1} ) k=1 n/2 ~ sum( delta_x * (1/3) * ( f(x_{2k-2}) + 4 f(x_{2k-1}) + f(x_{2k}) ) ) ~ k=1 delta_x n/2 = ------- sum( f(x_{2k-2}) + 4 f(x_{2k-1}) + f(x_{2k}) ) 3 k=1 delta_x n/2 n/2-1 = ------- ( f(x_0) + 4 sum( f(x_{2k-1}) ) + 2 sum( f(x_{2k}) ) + f(x_n) ) 3 k=1 k=1