Previous Contents Next

Chapter 7  Integration

7.1  Numerical Integration:

We will now discuss certain ways to compute definite integrals.

7.2  Newton-Cotes Formulae:

In this technique, we approximate òab f(x) dx using an appropriate Newton Interpolating polynomial.

7.2.1  The Basics:

Consider the problem of integrating the function f on the finite interval [a,b]. Let's assume that the interval [a,b] has been subdivided into points {x1, ..., xm} which are uniformly spaced. We let the function values at these points be given by fi = f(xi). This gives
xi =
a +
i-1
m-1
(b-a)
x1 = a
x2 =
a + 1 × h, where h =
m-1
b-a
··· = ···
xm = a + (m-1) × h

For convenience, let's recall our Newton Interpolant discussions for the case of a four point partition of [a,b], {x1, x2, x3, x4 } as given above. Hence h = b-a/3 here. The Newton interpolating polynomial in this case is given by
p3(x) = c1 + c2 (x-x1) + c3 (x-x1) (x-x2) + c4 (x-x1)(x-x2)(x-x3)

where the coefficients are given by
c1 = f1
c2 =
f2 - f1
x2 - x1
c3 =
f3 - f1
x3 - x1
-
f2 - f1
x2 - x-1
x3-x2


and c4, being a bit more complicated, will be shown in parts: let
y421 =
f4-f1
x4-x1
-
f2-f1
x2-x1
x4-x2
y321 =
f3-f1
x3-x1
-
f2-f1
x2-x1
x3-x2

then
c4 =
y421 - y321
x4-x3

This can be then be rewritten as
c1 = f1
c2 =
f2 - f1
h
c3 =
f3 - f1
2h
-
f2 - f1
h
h
  =
f3 - 2f2 + f1
2h2
y421 =
f4-f1
3h
-
f2-f1
h
2h
=
2f1-3f2+f4
6h2
y321 =
f3-f1
2h
-
f2-f1
h
h
  =
f1-2f2+f3
2h2
c4 =
y421 - y321
h
  =
f4-3f3+3f2-f1
6h3

It is easy to compute the coefficients efficiently via a Matlab loop like 
m = length(x);
for k = 1:m-1;
 f(k+1:m) = (f(k+1:m) - f(k) )./(x(k+1:m) - x(1) );
end
c = f;
We will approximate òab f(x) ds by replacing f by a Newton interpolant on the interval [a,b]. For an m point unifromly spaced partition, we will use the m-1 degree polynomial pm-1 whose coefficients can be computed recursively like we do above. Let's introduce the product notation
Pi=1k-1 (x-xi) =
1, i = 1
(x-x1) i = 2
(x-x1) (x-x2) i = 3
··· ···
(x-x1) (x-x2) ··· (x-xk-1) i = k-1

Then
pm-1(x) =
m
å
k=1
ck Pi=1k-1 (x-xi)

and
ó õ
b


a
f(x) dx
approx
ó õ
b


a
pm-1(x) dx
  =
ó õ
b


a
æ
ç
ç
è
m
å
k=1
ck Pi=1k-1 (x-xi) ö
÷
÷
ø		 dx
=
m
å
k=1
ck æ
ç
ç
ç
ç
è		 ó õ
b


a
Pi=1k-1 (x-xi) ö
÷
÷
÷
÷
ø		 dx

To figure out this approximation, we need to evaluate
ó õ
b


a
Pi=1k-1 (x-xi) dx

Make the change of variable s defined by x = a + sh. Then we find x - xi becomes s - i + 1 and
ó õ
b


a
Pi=1k-1 (x-xi) dx
=
ó õ
m-1


0
hk Pi=1k-1 (s-i+1) ds

We will let
Smk =
ó õ
m-1


0
Pi=1k-1 (s-i+1) ds

which gives us the final form of our approximation:
ó õ
b


a
f(x) dx
approx
k
å
1
ck hk Smk

7.2.2  Evaluating Smk:

The value of m we choose to use gives rise then to what is called an m point rule and given m we can easily evaluate the needed coefficients Sm1 through Smm.

Here are some calculations:
Sm1 =
ó õ
m-1


0
1 ds
  = (m-1)
Sm2 =
ó õ
m-1


0
s ds
  =
(m-1)2
2
Sm3 =
ó õ
m-1


0
s(s-1) ds
  =
(m-1)2
6
(2m - 5)
Sm4 =
ó õ
m-1


0
s(s-1)(s-2) ds
  =
(m-1)2 (m-3)2
4

We will denote this Newton Polynomial approximation to this integral using m points by the symbol QNC(m). Hence, for the case m=4, we have
Sm1 = 3
Sm2 =
9
2
Sm3 =
9
2
Sm4 =
9
4

This leads to the approximation
ó õ
b


a
f(x) dx
approx c1 h S41 + c2 h S42 + c3 h S43 + c4 h S44
  =
3 f1 h +
9
2
f2 - f1
h
h2 +
9
2
f3 - 2f2 + f1
h2
h3 +
9
4
f4-3f3+3f2-f1
6h3
h4
  =
3 f1 h +
9
2
(f2 - f1) h +
9
4
(f3 - 2f2 + f1) h +
9
24
(f4 - 3f3 + 3f2 - f1) h
  =
3 f1 h +
9
2
(f2 - f1) h +
9
4
(f3 - 2f2 + f1) h +
3
8
(f4-3f3+3f2-f1) h
  =
3h
8
( 8f1 + 12f2 - 12f1 6f3 - 12f2 + 6f1 f4-3f3+3f2-f1 )
  =
3h
8
( f1 + 3f2 + 3f3 + f4 )

But h is b-a/3 here, so our final 4 point formula is
ó õ
b


a
f(x) dx
approx
b-a
8
( f1 + 3f2 + 3f3 + f4 )

7.2.3  Exercise:

  1. Show that for m = 2, we get the Trapezoidal Rule
    ó õ
    b


    a
    		f(x) dx
    		approx
    b-a æ
    ç
    ç
    è
    1
    2
    		f1 +
    1
    2
    		f2 ö
    ÷
    ÷
    ø


  2. Show that for m = 3, we get the Simpson Rule
    ó õ
    b


    a
    		f(x) dx
    		approx
    b-a
    6
    		( f1 + 4f2 + f3 )

7.2.4  matlab Implementation:

We store the Newton-Cotes Weight vectors in this short piece of Matlab code

WNC.m 
 1 function w = WNC(m)
 2 %
 3 % m   an integer from 2 to 4
 4 % w   this is the Newton Cotes Weight Vector
 5 %
 6 if m==2
 7   w = [1 1]'/2;
 8 elseif m==3
 9   w = [1 4 1]'/6;
10 elseif m==4
11   w = [1 3 3 1]'/8;
12 else
13   disp('You must use a value of m between 2 and 4');
14   w = [0 0]';
15 end
syntax highlighted by Code2HTML, v. 0.8.8b

The Newton-Cotes implementation is then given by

QNC.m 
 1 function numI = QNC(fname,a,b,m)
 2 %
 3 % fname  the name of the function which is to integrated
 4 % a,b    the intergation interval is [a,b]
 5 % m      the order of the Newton-Cotes Method to use
 6 %        this is the same as the number of points in
 7 %        the partition of [a,b]
 8 %
 9 % numI   the value of the Newton-Cotes approximation to
10 %        the integral
11 %
12 if m>= 2 & m <= 4
13   w = WNC(m);
14   x = linspace(a,b,m)';
15   f = feval(fname,x);
16   numI = (b-a)*(w'*f);
17 else
18   disp('You need to use an order m between 2 and 4');
19 end
syntax highlighted by Code2HTML, v. 0.8.8b

and finally, a script to run the numerical integration routines is as follows:

ShowNC.m 
 1 while input('Another Example? (1=yes, 0=no). ');
 2   fname = input('Enter function name: ');
 3   a = input('Enter left endpoint ');
 4   b = input('Enter right endpoint: ');
 5   s = ['QNC(' fname sprintf(',%6.3f,%6.3f,m )',a,b)];
 6   disp(['  m      ' s])
 7   disp(' ')
 8   for m = 2:4
 9     numI = QNC(fname,a,b,m);
10     disp(sprintf(' %2.0f      %20.16f',m,numI))
11   end
12 end
syntax highlighted by Code2HTML, v. 0.8.8b

7.2.5  Run Time Output:

We use the functions defined in func1.m, (f(x) = sin(x)), and func2.m, (f(x) = exp-x2), which have matlab codes

func1.m 
1 function y = func1(x)
2 y = sin(x);
syntax highlighted by Code2HTML, v. 0.8.8b

func2.m 
1 function y = func2(u)
2 %
3 %
4 %
5 z = -u.*u;
6 y = exp(z);
syntax highlighted by Code2HTML, v. 0.8.8b Here is our runtime

NC.prompts 
 1 Enter function name: 'func1'
 2 Enter left endpoint 0
 3 Enter right endpoint: pi/2
 4   m      QNC(func1, 0.000, 1.571,m )
 5  
 6   2        0.7853981633974483
 7   3        1.0022798774922104
 8   4        1.0010049233142790
 9 Another Example? (1=yes, 0=no). 1
10 Enter function name: 'func2'
11 Enter left endpoint 0
12 Enter right endpoint: 2.0
13   m      QNC(func2, 0.000, 2.000,m )
14  
15   2        1.0183156388887342
16   3        0.8299444678581679
17   4        0.8622241875991991
18 Another Example? (1=yes, 0=no). 0
syntax highlighted by Code2HTML, v. 0.8.8b

Previous Contents Next