Fun with numerical integration « Alasdair's musings

May 2009
M	T	W	T	F	S	S
« Apr		Jun »
	1	2	3
4	5	6	7	8	9	10
11	12	13	14	15	16	17
18	19	20	21	22	23	24
25	26	27	28	29	30	31

Fun with numerical integration

Posted on May 23, 2009 by amca01

Numerical integration, or quadrature (to use the old-fashioned, but charming, term) is the business of finding approximate values to definite integrals; integrals which are either impossible to be computed analytically, or for which the antiderivative is too complex to be used. In the first category we might have

$\int\sin(x+x^3)\,dx$

At least, Wolfram|Alpha can’t compute this, so I assume it can’t be done analytically; and in the second category we might have

$\int\sin(1+x^3)\,dx$

for which the antiderivative is a horrific expression involving the incomplete gamma function.

In an elementary calculus course, students might be introduced to some very simple numerical techniques: the trapezoidal rule or Simpson’s rule, but unless the course is more specifically geared to numerical computations, not much further.

However, there are lots of different methods for numerical integration, some of which are great fun to play with.

Newton-Cotes integration

All these methods are based on one simple idea: sample the function values at equidistant points, find the interpolating polynomial through those points, and integrate the polynomial. Since polynomials are easy to integrate, this method is easy to apply. For example, let’s generate the fourth-order Newton-Cotes rule, which uses a quartic polynomial, and hence requires five points.

Using Maxima:

(%i1) load(interpol); (%i2) pts:[[a,f1],[a+h,f2],[a+2*h,f3],[a+3*h,f4],[a+4*h,f5]]; (%i3) integrate(lagrange(pts),x,a,a+4*h); (%i4) ratsimp(%); (%o4) (14 f5 + 64 f4 + 24 f3 + 64 f2 + 14 f1) h/45

In other words,

$\displaystyle{\begin{array}{rcl}\int_a^{a+4h}f(x)\,dx&\approx&\frac{2h}{45}(7f(a)+32f(a+h)+12f(a+2h)\\&&\qquad+32f(a+3h)+7f(a+4h))\end{array}}$ .

and this is known as Boole’s rule.

The trouble with single Newton-Cotes rules is that as the number of points get larger, the interpolating polynomial may develop more “wiggles” and the resulting integral approximation may get worse. For example, take the integral

$\displaystyle\int^5_{-5}\frac{1}{1+x^2}\,dx$

for which the exact value is

$2\tan^{-1}{5}\approx 2.746801533890031$

However, let’s try his with a tenth-order rule:

(%i5) f(x):=1/(1+x^2); (%i6) pts:makelist([-5+n,f(-5+n)],n,0,10); (%i7) integrate(lagrange(pts),x,-5,5),numer; (%o7) 4.673300563481597

We can see why this is so poor by plotting the original function and its interpolant together (click on image to show full-sized image):

For this reason, better results can be obtained by chopping the integral up into small segments, and applying a low order rule (say, Simpson’s) to each segment.

Newton-Cotes with adjustments

If you look up old texts, or search about on the web, you’ll find rules similar to the Newton-Cotes rules, but with much simpler coefficients.

For example, the sixth-order rule is

$\displaystyle{\begin{array}{rcl}\int_a^{a+6h}f(x)\,dx&\approx&\frac{h}{140}(41f(a)+216f(a+h)+27f(a+2h)\\&&\qquad+272f(a+3h)+27f(a+4h)\\&&\qquad+216f(a+5h)+41f(a+6h))\end{array}}$ .

As you can see, the coefficients are already getting cumbersome.

One way of adjusting a Newton-Cotes formula is by adding some function differences. If we denote $f(a+nh)$ by $f_n$ these are defined as

$\Delta f_n=f_{n+1}-f_n$

and all higher differences can be obtained recursively for $k\ge 1$ :

$\Delta^{k+1} f_n=\Delta^kf_{n+1}-\Delta^kf_n$ .

It is not hard to show that

$\Delta^kf_n=\sum_{p=0}^k(-1)^p{n\choose p}f_{n+p}$

So, for example

$\Delta^6f_0=f_6-6f_5+15f_4-20f_3+15f_2-6f_1+f_0$ .

The idea is that for most functions, differences tend to get smaller, so by adding some small multiple of a difference, we shouldn’t affect the formula too much, but we may simplify its coefficients. Let’s try this with the Newton-Cotes sixth-order rule, which for simplicity we shall integrate over the range $[1,7]$ :

(%i8) [a,h]:[1,1]; (%i9) remfunction(f); (%i10) pts:makelist([a+n*h,f(a+n*h)],n,0,6); (%i11) integrate(lagrange(pts),x,1,7); (%i12) nc6:ratsimp(%);

This just produces the sixth-order rule seen above. But now!

(%i13) d6:sum(binomial(6,k)*(-1)^k*f(a+k*h),k,0,6); (%o13) f(7) - 6 f(6) + 15 f(5) - 20 f(4) + 15 f(3) - 6 f(2) + f(1) (%i14) ratsimp(nc6+d6/140); (%o14) (3 f(7) + 15 f(6) + 3 f(5) + 18 f(4) + 3 f(3) + 15 f(2) + 3 f(1))/10

This is Weddle’s rule; more commonly written as

$\displaystyle{\begin{array}{rcl}\int_a^{a+6h}f(x)\,dx&\approx&\frac{3h}{10}(f(a)+5f(a+h)+f(a+2h)\\&&\qquad+6f(a+3h)+f(a+4h)\\&&\qquad+5f(a+5h)+f(a+6h))\end{array}}$ .

Notice how much simpler the coefficients are! Another rule, Hardy’s rule, can also be obtained by an adjustment:

(%i15) ratsimp(nc6-d6*9/700); (%o15) (14 f(7) + 81 f(6) + 110 f(4) + 81 f(2) + 14 f(1))/50

or in integral form as

$\displaystyle{\begin{array}{rcl}\int_a^{a+6h}f(x)\,dx&\approx&\frac{h}{50}(14f(a)+81f(a+h)+110(a+4h)\\&&\qquad+81f(a+5h)+14f(a+6h)\end{array}}$ .

Not only is this simpler than the sixth-order Newton-Cotes rule, but two of the coefficients have vanished.

Averaging different rules

You and your students can invent your own rules by taking some Newton-Cotes rules and forming weighted averages of them. For example, Simpson’s rule (the second order rule) over a difference of $2h$ is

$\displaystyle{\int_a^{a+4h}f(x)\,dx=\frac{2h}{3}(f(a)+4f(a+2h)+f(a+4h))}.$

If we, for example, take

$\displaystyle{\frac{5}{3}N-\frac{2}{3}S}$

where $N$ is the fourth order Newton-Cotes rule, and $S$ is the Simpson rule over $2h$ , we obtain

$\displaystyle{\int_a^{a+4h}f(x)\,dx\approx\frac{2h}{9}(f(a)+8f(a+h)+8f(a+3h)+f(a+4h))}$

which certainly has nice coefficients.

Such a formula may not in fact give very accurate results, but students are always excited to discover something “for themselves”, and inventing their own integration rules (and testing them), would be a fascinating exercise.

Filed under: Computation, Maths teaching, Maxima

« A weekend with Wolfram Alpha The “Mathematics Change Plan” »

Alasdair's musings

Recent Posts

Categories

Blog Stats

Visitors

Pages

Archives

Blogroll

Meta

Fun with numerical integration

2 Responses

Leave a Reply

Calendar

Blogroll