Integration

In addition to calculating derivatives, the ability to calculate an integral is a fundamental task in physics. As with differentiation, the task is not as straightforward as we would hope. The reason for this lies in the same regime as it did with differentiation, namely the use of a limiting process. Recall that the definition of an integral is

,

(6.1)

where x₀ = a and x_N = b. The integral is found by taking the limits of the number of strips going to infinity while the strip width goes to zero. Notice that there is no requirement that the strip width be the same for each strip. Also recall that by this definition, the integral is simply the area enclosed by the boundaries x = a, x = b, y = 0, and y = f(x). So, if we can find methods that allow us to calculate this area fairly accurately, we can determine the integral without having to deal with the limit.

Quadrature Integration

One of the easiest methods to calculate an integral is to evaluate the function at a number of specific locations and then use the results of these evaluations to approximate the integral. We accomplish this be assigning a weight to each result and then summing all of these weighted results together

,

(6.2)

where x_i are the evaluation points, w_i is the weight given to the i^th point, and there are N+1 points evaluated. For now, let’s consider the case where the evaluation points are evenly spaced with a separation of h.

In order to use equation (6.2) on a wide range of integrals, we want to require that equation (6.2) be exact as f(x) is sequentially approximated by higher order polynomials

.

(6.3)

Substituting these into equation (6.2) gives a series of simultaneous equations in w_i, which can be solved for the various weights.

Example:

The simplest solution of equation (6.2) is found by using the first two formulae in equation (6.3). Thus, for f(x) = 1 we get

,

while f(x) = x yields

.

Solving these two equations for w₀ and w₁, we find

.

Setting the step size, h = b – a, the two point approximation for equation (5.20) becomes

,

(6.4)

where the truncation error is shown in the last two lines. This is known as the trapezoid rule.

Extending this to three formulae, we get for

f(x) = 1:	,
f(x) = x:	,
f(x) = x^2:	,

Where a step size of h = (b – a)/2 has been used. This is another set of (more complicated) simultaneous equations, which can be solved to yield

.

(6.5)

This is called Simpson’s rule.

This can be continued for four and five formulae, yielding Simpson’s three eighths rule,

,

(6.6)

and Milne’s rule,

.

(6.7)

It’s important to realize that for each rule the step size is defined to be , where N is the number of points used in the evaluation!

From this example, we see that we are able to get better accuracy with more points, as we would expect. Unfortunately, the solution of the simultaneous equations become more complicated as the number of points increases, since the number of equations grows linearly with the number of points. Thus, while we could, in theory, use this method to solve any integral, the computation involved to go beyond Milne’s rule is often not worth the effort.

Comparing the different algorithms, we see that they are simply different ways of choosing the points and weights. In general, the precision increases as N gets larger, until round-off error eventually limits the process. However, because the “best” approximation depends upon the specific behavior of f(x), there isn’t a universally “best” approximation.

When the intervals are kept constant, the quadrature integration methods are known as Newton-Cotes methods. However, there is no requirement that the steps remain equal. If this assumption is relaxed, we can improve on our previous results. These Gaussian quadrature methods enable us to integrate exactly (within the limits of round-off error, of course) the product of a function times a (2N – 1) degree polynomial, with only N function evaluations. To see this, first consider equation (6.2) again, rewriting it as

,

(6.8)

where the sum is now started at m = 1 so that N refers to the number of function evaluations being made. As in the Newton-Cotes methods, the w_m values are unknown, but now the x_m are also unknown. This means that we now have 2N unknowns, and thus require 2N equations. We can again obtain these equations by requiring that the quadrature be exact for the lowest order polynomials, f(x) = 1, x, x², …, x^2N-1. Unfortunately, these equations are nonlinear and extremely difficult to solve.

Example:

Let N = 2 and require the integration formula be exact for f(x) = 1, x, x², and x³. Then equation (6.8) yields the equations

.

In order to solve these equations, note that any finite interval [a,b] can be mapped onto the interval [-1, 1] by a change of variable

.

(6.9)

Thus, we only need to consider this normalized integration region, in terms of which the nonlinear equations become

.

(6.10)

Combining the second and fourth equations yields . Since we need to have four independent variables, each point must be distinct, which implies that y₁ = -y₂. Substituting this into the second equation shows that w₁ = w₂, which, from the first equation, shows that w₁ = 1. Substituting these into the third equation then shows that .

As can be seen, N = 2 is complicated enough, applying this to more points gets even harder very quickly. This is where the (2N – 1) degree polynomial comes in, but before we can use it, we need to understand a little about orthogonal polynomials.

Orthogonal Polynomials

The basic idea behind the theory of orthogonal polynomials is that there exists (or if it doesn’t that we can build) a set of polynomials f_m(x) such that

,

(6.11)

where d_mn is the Kronecker delta and C_m is a constant. Equation (6.11) tells us that the set of polynomials form an independent set of functions which can be used to describe the N independent coordinate axes associated with our N dimensional function space.

Example:

The Gram-Schmidt orthogonalization procedure can be used to build up our orthogonal functions. In general, the m^th function can be found by forming the sum

,

(6.12)

where the a_mi are determined by the orthogonality condition (6.11), and f₀(x) is chosen based on the values of W(x), a and b.

Consider the polynomials u_m = x^m. Setting W(x) = 1, and letting a = -1 and b = 1, chose f₀(x) be constant. Then if we normalize f₀(x), we get

.

(6.13)

The next term would be

f₁(x) = u₁ + a₁₀ f₀(x).

Using equation (6.11) to find a₁₀ yields

,

(6.14)

which implies a₁₀ = 0 and thus f₁(x) = x. This is then normalized by requiring

,

so that .

In a similar manner, the third function is given by

f₂(x) = u₂ + a₂₁ f₁(x) + a₂₀ f₀(x).

Now we require that f₂(x) satisfy the three equations

,

,

and

.

This procedure can be repeated for higher order terms until all 2N orthogonal functions can be found.

Gaussian Quadrature

We are now ready to apply orthogonal functions to problem of performing numerical integration. Begin by rewriting equation (6.8) as

,

(6.15)

where W(x) is a positive definite weighting function. As before, neither the weights, w_m, nor the evaluation points, x_m, are known, so there are 2N values to be determined. However, this time we will require that this expansion be exact for polynomials of order 2N – 1 and less, and use this to determine the weights and evaluation points.

This can be done by letting f(x) be a polynomial of degree 2N – 1 or less, and f_N(x) be a specific orthogonal function of order N such that for the particular weighting function W(x) and limits of integration a and b,

.

(6.16)

To find the evaluation points, divide f(x) by f_N(x). The leading term in the quotient function, as well as the remainder function, will be of order N – 1. In terms of this quotient and remainder, the original function can be written as

,

(6.17)

where both the quotient function, q_N-1(x), and the remainder function, r_N-1(x), are polynomials of order N – 1.

Substituting equation (6.17) into equation (6.15) yields

.

(6.18)

We can simplify the first integral by using the fact that the functions f_N(x), being an orthogonal set of functions, span the N dimensional function space. Thus, we can write q_N-1(x) in terms of these functions as

,

(6.19)

where the q_i are constants. The first integral in equation (6.18) then becomes

,

(6.20)

where the last equality follows from the fact that the summation only runs to N – 1.

Now we reach the key step in this process. Recall that we started by requiring that equation (6.15) be exact for f(x), which was an arbitrary function of order 2N – 1. Since f(x) is arbitrary, it must also hold for the polynomial q_N-1(x)f_N(x), which is also of order 2N – 1. Thus, equation (6.15) tells us that

.

(6.21)

But equation (6.20) shows that the left hand side of equation (6.21) is zero. Since f(x) is an arbitrary function, the quotient function q_N-1(x) is also arbitrary, so the sum on the right hand side of equation (6.21) can only be guaranteed to vanish if the f_N(x_m) are all zero. Thus, through the correct choice of f_N(x), all of the evaluation points x_m are determined.

We still need to find the weights, w_m. Since equation (6.15) is required to be exact for polynomials of order 2N – 1, it must also be exact for polynomials of lesser order as well. In particular, it must be true for the (N – 1)^th order polynomial l_i,N(x), defined as

,

(6.22)

with j running from 1 to N. One of the properties of this function is that

.

Thus,

.

(6.23)

So the weights w_i can be determined by performing the integration in equation (6.23).

Example:

What would be the evaluation points and weights for the integral , if we limit ourselves to two function evaluations?

For this integral a = -1, b = 1 and W(x) = 1. The orthogonal functions associated with these values was calculated in a previous example to be

.

If we only want to do two function evaluations, the evaluation points are determined by finding the zeros of f₂

.

(6.24)

The weights are found from equation (6.23) to be

(6.25)

and

.

(6.26)

In general, results found via Gaussian quadrature are better than the results obtained via Newton-Cotes methods, so long as there are no singularities in the integrand or in its derivative. When there are Simpson’s rule can be used to avoid problems, but in general, we are better off removing singularities analytically from the integration before attempting any numerical method. This can be done by breaking the interval into subintervals, so the singularity is at an endpoint where an evaluation point never falls, or by changing the integration variable. Similarly, if the integrand has a region in which it only has slow variations, a change of variable can compress the region so that there are less evaluation points in that region. Conversely, if the integrand varies rapidly in some region, a change of variables can expand that region so that more evaluation points are included.

Euler-McClaurin Integration

Let’s return to equation (6.1). If we expand f(x) as a Taylor series about x = a, we get

.

(6.27)

Similarly, if we expand about x = b,

.

(6.28)

where h = b – a. We can combine these two equations by multiplying each by ½ and adding them together. In this case

.

(6.29)

This equation is useful if we eliminate the even order derivatives. To accomplish this, recall that a Taylor series expansion of f’(x = b) about x = a yields

,

(6.30)

while the expansion of f’(x = a) about x = b is

.

(6.31)

Substituting equation (6.31) into equation (6.30) yields

.

(6.32)

In a similar way, f^(2N)(a) + f^(2N)(b) , with N = 1, 2, 3, …, can be found by expanding f^(2N‑1)(x) about x =a and x = b. Substituting these solutions back into equation (6.29) eliminates the even derivatives to yield

.

(6.33)

The important thing to note is that the derivatives show up as differences. If we divide the integration region into many smaller segments and apply equation (6.33) to each of these segments, then with the exception of points a and b the derivatives contribute to the two adjacent segments in such a way that their contribution to the total integral vanishes. Thus, the Euler-McClaurin rule can be written as

.

(6.34)

Notice that the first part is simply the trapezoid rule, but now we have correction terms, so we can determine the truncation error associated with this rule. When the integrand is easily differentiated, the Euler-McClaurin scheme is vastly superior to other integration schemes. It is particularly powerful when the derivatives vanish at the limits of integration, so that the “simple” trapezoid rule can yield surprisingly accurate results.