Loading [MathJax]/jax/output/HTML-CSS/jax.js
Skip to main content
Library homepage
 

Text Color

Text Size

 

Margin Size

 

Font Type

Enable Dyslexic Font
Physics LibreTexts

1.2: Numerical Integration

( \newcommand{\kernel}{\mathrm{null}\,}\)

There are many occasions when one may wish to integrate an expression numerically rather than analytically. Sometimes one cannot find an analytical expression for an integral, or, if one can, it is so complicated that it is just as quick to integrate numerically as it is to tabulate the analytical expression. Or one may have a table of numbers to integrate rather than an analytical Equation. Many computers and programmable calculators have internal routines for integration, which one can call upon (at risk) without having any idea how they work. It is assumed that the reader of this chapter, however, wants to be able to carry out a numerical integration without calling upon an existing routine that has been written by somebody else.

There are many different methods of numerical integration, but the one known as Simpson's Rule is easy to program, rapid to perform and usually very accurate. (Thomas Simpson, 1710 - 1761, was an English mathematician, author of A New Treatise on Fluxions.)

Suppose we have a function y(x) that we wish to integrate between two limits. We calculate the value of the function at the two limits and halfway between, so we now know three points on the curve. We then fit a parabola to these three points and find the area under that.

In the figure I.1, y(x) is the function we wish to integrate between the limits x2δx and x2+δx. In other words, we wish to calculate the area under the curve. y1, y2 and y3 are the values of

alt
FIGURE I.1 Simpson's Rule gives us the area under the parabola (dashed curve) that passes through three points on the curve y=y(x). This is approximately equal to the area under y=y(x).

the function at x2δx, x2 and x2+δx, and y=a+bx+cx2 is the parabola passing through the points (x2δx,y1), (x2,y2) and (x2+δx,y3).

If the parabola is to pass through these three points, we must have

y1=a+b(x2δx)+c(x2δx)2

y2=a+bx+cx2

y3=a+b(x2+δx)+c(x2+δx)2

We can solve these Equations to find the values of a, b and c. These are

a=y2x2(y3y1)2δx+x22(y32y2+y1)2(δx)2

b=y3y12δxx2(y32y2+y1)(δx)2

c=y32y2+y12(δx)2

Now the area under the parabola (which is taken to be approximately the area under y(x)) is

x2+δxx2δx(a+bx+cx2)dx=2[a+bx2+cx22+13c(δx)2]δx

On substituting the values of a, b and c, we obtain for the area under the parabola

13(y1+4y2+y3)δx

and this is the formula known as Simpson's Rule.

For an example, let us evaluate π/20sinxdx.

We shall evaluate the function at the lower and upper limits and halfway between. Thus

x=0,y=0x=π/4,y=1/2x=π/2,y=1

The interval between consecutive values of x is δx=π/4.

Hence Simpson's Rule gives for the area

13(0+42+1)π4

which, to three significant figures, is 1.00. Graphs of sinx and a+bx+cx2 are shown in figure I.2a. The values of a, b and c, obtained from the formulas above, are

a=0,b=322π,c=8128π2

alt
FIGURE I.2a

The result we have just obtained is quite spectacular, and we are not always so lucky. Not all functions can be approximated so well by a parabola. But of course the interval δx=π/4 was ridiculously coarse. In practice we subdivide the interval into numerous very small intervals. For example, consider the integral

π/40cos322xsinxdx.

Let us subdivide the interval 0 to π/4 into ten intervals of width π/40 each. We shall evaluate the function at the end points and the nine points between, thus:

xcos32xsinxdx0y1=0.000 000 000π/40y2=0.077 014 6222π/40y3=0.145 091 4863π/40y4=0.196 339 0024π/40y5=0.224 863 4305π/40y6=0.227 544 9306π/40y7=0.204 585 4737π/40y8=0.159 828 8778π/40y9=0.100 969 9719π/40y10=0.040 183 06610π/40y11=0.000 000 000

The integral from 0 to 2π/40 is 13(y1+4y2+y3)δx, δx being the interval π/40. The integral from 3π/40 to 4π/40 is 13(y3+4y4+y5)δx. And so on, until we reach the integral from 8π/40 to 10π/40. When we add all of these up, we obtain for the integral from 0 to π/4,

13(y1+4y2+2y3+4y4+2y5+......+4y10+y11)δx

=13[y1+y11+(y2+y4+y6+y8+y10)+2(y3+y5+y7+y9)]δx,

which comes to 0.108 768 816.

We see that the calculation is rather quick, and it is easily programmable (try it!). But how good is the answer? Is it good to three significant figures? Four? Five?

Since it is fairly easy to program the procedure for a computer, my practice is to subdivide the interval successively into 10, 100, 1000 subintervals, and see whether the result converges. In the present example, with N subintervals, I found the following results:

Nintegral100.108 768 8161000.108 709 62110000.108 709 466100000.108 709 465

This shows that, even with a course division into ten intervals, a fairly good result is obtained, but you do have to work for more significant figures. I was using a mainframe computer when I did the calculation with 10000 intervals, and the answer was displayed on my screen in what I would estimate was about one fifth of a second.

There are two more lessons to be learned from this example. One is that sometimes a change of variable will make things very much faster. For example, if one makes one of the (fairly obvious?) trial substitutions y=cosx, y=cos2x or y2=cos2x, the integral becomes

11/2(2y21)3/2dy,10y38(1+y)dyor10y42(1+y2)dy.

Not only is it very much faster to calculate any of these integrands than the original trigonometric expression, but I found the answer 0.108 709 465 by Simpson's rule on the third of these with only 100 intervals rather than 10,000, the answer appearing on the screen apparently instantaneously. (The first two required a few more intervals.)

To gain about one fifth of a second may appear to be of small moment, but in truth the computation went faster by a factor of several hundred. One sometimes hears of very large computations involving massive amounts of data requiring overnight computer runs of eight hours or so. If the programming speed and efficiency could be increased by a factor of a few hundred, as in this example, the entire computation could be completed in less than a minute.

The other lesson to be learned is that the integral does, after all, have an explicit algebraic form. You should try to find it, not only for integration practice, but to convince yourself that there are indeed occasions when a numerical solution can be found faster than an analytic one! The answer, by the way, is 18ln(1+2)216.

You might now like to perform the following integration numerically, either by hand calculator or by computer.

20x2dx2x

At first glance, this may look like just another routine exercise, but you will very soon find a small difficulty and wonder what to do about it. The difficulty is that, at the upper limit of integration, the integrand becomes infinite. This sort of difficulty, which is not infrequent, can often be overcome by means of a change of variable. For example, let x=2sin2θ, and the integral becomes

82π/20sin5θdθ

and the difficulty has gone. The reader should try to integrate this numerically by Simpson's rule, though it may also be noted that it has an exact analytic answer, namely 8192/15.

Here is another example. It can be shown that the period of oscillation of a simple pendulum of length l swinging through 90 on either side of the vertical is

P=8lgπ/20secθdθ.

As in the previous example, the integrand becomes infinite at the upper limit. I leave it to the reader to find a suitable change of variable such that the integrand is finite at both limits, and then to integrate it numerically. (If you give up, see Section 1.13.) Unlike the last example, this one has no simple analytic solution in terms of elementary functions. It can be written in terms of special functions (elliptic integrals) but they have to be evaluated numerically in any case, so that is of little help. I make the answer

P=2.3607πlg.

For another example, consider

0dxx5(e1/x1)

This integral occurs in the theory of blackbody radiation. To help you to visualize the integrand, it and its first derivative are zero at x=0 and x= and it reaches a maximum value of 21.201435 at x=0.201405. The difficulty this time is the infinite upper limit. But, as in the previous two examples, we can overcome the difficulty by making a change of variable. For example, if we let x=tanθ, the integral becomes

π/20c3(c2+1)dθec1,where c=cotθ=1/x.

The integrand is zero at both limits and is easily calculable between, and the value of the integral can now be calculated by Simpson's rule in a straightforward way. It also has an exact analytic solution, namely π4/15, though it is hard to say whether it is easier to arrive at this by analysis or by numerical integration.

Here’s another:

0x2dx(x2+9)(x2+4)2

The immediate difficulty is the infinite upper limit, but that is easily dealt with by making a change of variable: x=tanθ. The integral then becomes

π/2θ=0t(t+1)dθ(t+9)(t+4)2

in which t=tan2θ. The upper limit is now finite, and the integrand is easy to compute - except, perhaps, at the upper limit. However, after some initial hesitation the reader will probably agree that the integrand is zero at the upper limit. The integrand looks like this:

alt

It reaches a maximum of 0.029 5917 at θ=71.789 962. Simpson’s rule easily gave me an answer of 0.015 708. The integral has an analytic solution (try it) of π/200.

There are, of course, methods of numerical integration other than Simpson’s rule. I describe one here without proof. I call it “seven-point integration”. It may seem complicated, but once you have successfully programmed it for a computer, you can forget the details, and it is often even faster and more accurate than Simpson’s rule. You evaluate the function at 6n+1 points, where n is an integer, so that there are 6n intervals. If, for example, n=4, you evaluate the function at 25 points, including the lower and upper limits of integration. The integral is then:

baf(x)dx=0.3×(Σ1+2Σ2+5Σ3+6Σ4)δx,

where δx is the size of the interval, and

Σ1=f1+f3+f5+f9+f11+f15+f17+f21+f23+f25,

Σ2=f7+f13+f19,

Σ3=f2+f6+f8+f12+f14+f18+f20+f24

and Σ4=f4+f10+f16+f22.

Here, of course, f1=f(a) and f25=f(b). You can try this on the functions we have already integrated by Simpson’s rule, and see whether it is faster.

Let us try one last integration before moving to the next section. Let us try

10011+8x3dx.

This can easily (!) be integrated analytically, and you might like to show that it is

112ln147127+112tan1507+π432=0.6039748.

However, our purpose in this section is to learn some skills of numerical integration. Using Simpson’s rule, I obtained the above answer to seven decimal places with 544 intervals. With seven-point integration, however, I used only 162 intervals to achieve the same precision, a reduction of 70. Either way, the calculation on a fast computer was almost instantaneous. However, had it been a really lengthy integration, the greater efficiency of the seven point integration might have saved hours. It is also worth noting that x×x×x is faster to compute than x3. Also, if we make the substitution y = 2x, the integral becomes

2000.51+y3dy.

This reduces the number of multiplications to be done from 489 to 326 – i.e. a further reduction of one third. But we have still not done the best we could do. Let us look at the function 0.51+y3, in figure I.2b:


FIGURE I2b
alt

We see that beyond y=6, our efforts have been largely wasted. We don’t need such fine intervals of integration. I find that I can obtain the same level of precision − i.e. an answer of 0.6039748 − using 48 intervals from y=0 to 6 and 24 intervals from y=6 to 20. Thus, by various means we have reduced the number of times that the function had to be evaluated from our original 545 to 72, as well as reducing the number of multiplications each time by a third, a reduction of computing time by 91. This last example shows that it often advantageous to use fine intervals of integration only when the function is rapidly changing (i.e. has a large slope), and to revert to coarser intervals where the function is changing only slowly.

The Gaussian quadrature method of numerical integration is described in Sections 1.15 and 1.16.


This page titled 1.2: Numerical Integration is shared under a CC BY-NC 4.0 license and was authored, remixed, and/or curated by Jeremy Tatum via source content that was edited to the style and standards of the LibreTexts platform.

Support Center

How can we help?