- Math: Extension of the mean value theorem
- Computing: Loops
- Science: Radioactive decay and heat transfer
General Information: This section is a continuation of the previous section. We study the second order accurate trapezoid rule and numerical derivatives.
Contact Person: R. E. White, NCSU, white@math.ncsu.edu
Revision Date: 7-29-94
Copyright ©1994, R. E. White. All rights reserved.
2.3 Numerical Integration and the Trapezoid Rule
2.3.1. Introduction.
In this section we will give an extension of the mean value theorem which will allow us to establish second order convergence of the trapezoid algorithm for numerical integration. Second order convergence means the error is proportional to the square of the mesh, h. We will also use this theorem for numerical derivative approximation. In later sections we will also see how these results can be applied to radioactive decay problems and to a continuum model of heat diffusion.
The second form of the mean value theorem's conclusion is
f(x) = f(a) + f'(c) (x - a) (1)
f(x) = f(a) + f'(a) (x - a) + f''(c)(x - a)2/2 (2)
for some c between a and x, provided f''(x) is continuous on [a,b]. We will use these two representations of f(x) in the analysis of the following applications.
We present three applications. The first two are a continuation of those in the previous lecture, but here we seek more accurate numerical approximations. The third application involves numerical derivatives and we see how this is related to the previous discrete models for the diffusion of heat.
Application to Radioactive Decay. The integral definition of the ln(x) was approximated by using rectangular regions. But, the graph of f(t) = 1/t indicates that a much more accurate approximation can be gotten by using trapezoids, see the figure below. The area of a single trapezoid is given by the base times the average of left and right heights. And, the corresponding approximation of ln(3) with four trapezoids is
- 1/2 (1/1 + 1/(3/2))/2 + 1/2 (1/(3/2) + 1/2)/2 + 1/2 (1/2 + 1/(5/2))/2 + 1/2 (1/(2/5) + 1/3))/2 = 134/120.
 |
| Figure: Trapezoid Approximation |
Application to Queen's Gold Corrugated Roof. Consider the roof problem in the previous section. Suppose the roof is to be made from gold so that the size of the error that is tolerable decreases to 10-6 . One can take advantage of the symmetry of the sine curve and integrate over one fourth of its cycle. The surface area is then given by
Another move to reduce the amount of calculations is the use the trapezoid sum algorithm in place of the Riemann sum algorithm. Here we will approximate the area by trapezoids and not rectangles.
Application to Heat Diffusion and Numerical
Derivatives. The discrete form of Fourier's heat law states that the
heat flowing through a surface A is proportional to the product of A, change
in time and
Application of (2) gives the first order approximation of the derivative
- f(x + h) = f(x) + f'(x)h + f''(c) h2/2 and
- f(x - h) = f(x) - f'(x)h + f''(C) h2/2.
This then gives
A second order estimate of the derivative can be obtained, provided the function has a continuous third derivative. From the above note that
- (f(x + h) - f(x - h))/(2h) = f'(x) + (f''(c) - f''(C))h/4.
Now apply the mean value theorem to f''(x) to obtain
- (f(x + h) - f(x - h))/(2h) = f'(x) + f'''(C')(c - C) h/4.
Thus,
- |(f(x + h) - f(x - h))/(2h) - f'(x)| < (M3 /2)h2
The trapezoid sum algorithm for numerical integration is based on the approximation of the area by looking at trapezoids associated with portions of the area, provided the function is nonnegative. This approximation has the form
The following figure depicts this. One should note the small error for the trapezoid approximation, and contrast this with the much larger error for the rectangle approximation in the Riemann sum algorithm.
 |
| Figure: Trapezoid Sums |
The trapezoid sum algorithm for an integral is also easily programmed. Below we have made an adjustment to avoid computation of f(x(i)) twice. This minor change and the smaller errors for each trapezoid approximation makes the trapezoid algorithm much more efficient than the Riemann sum algorithm. There are many other numerical integration methods which have higher orders of convergence or are for special functions such as those with large oscillations.
h = (b - a)/n, x(i) = a, and sumt = f(x(i))/2
for i = 1, n - 1
x(i) = x(i) + h
sumt = sumt + f(x(i))
endloop
sumt = (sumt + f(x(i))/2)*h.
The following Maple code is a slight variation of the Maple code for the Riemann sum algorithm. However, note the more rapid convergence to the solution.
Maple Code for the Trapezoid Sum Algorithm.
Input Data:Next we implement the Matlab code and apply it to the numerical approximation in the surface area integral for the above roof problem.> f:=x->1/x;
> a:=1: > b:=3: > n:=4: Procedure TSUM is Defined: > TSUM:=proc(f,a,b,n) > h:= (b-a)/n: > x:=a: > sumt:= f(a)/2: > for i from 1 to n-1 do > x:= x + h: > sumt:= sumt + f(x): > od; > sumt:= sumt + f(b)/2: > sumt:= evalf(h*sumt); > end: Output Data: > evalf(ln(3)); 1.09861 > TSUM(f,a,b,4); 1.1667 > TSUM(f,a,b,8); 1.10321 > TSUM(f,a,b,16); 1.09977 > TSUM(f,a,b,32); 1.09890
Matlab Code for the Trapezoid Sums.
% Define f(x)=(1 + ((8*pi/12)*cos(8*pi*x))^2)^.5 in the m-file f.m
n = 8;
a = 0;
b = 1/16;
h = (b-a)/n;
x = a;
sumt = f(a)/2;
for i = 1:n-1
x = x + h;
sumt = sumt + f(x);
end
sumt = sumt + f(b)/2;
n
sumt = sumt*h
Output:
n =
2
sumt =
0.10773144259152
n =
4
sumt =
0.10792279645513
n =
8
sumt =
0.10792438244445
The second form of the mean value theorem's conclusion is
- f(x) = f(a) + f'(c) (x - a)
for some c between a and x, provided f'(x) is continuous on [a,b]. The extended mean value theorem will conclude that
- f(x) = f(a) + f'(a) (x - a) + f''(c)(x - a)2 /2
for some c between a and x, provided f''(x) is continuous on [a,b].
The extended mean value result follows in a natural manner from integration by parts and the integral form of the mean value theorem, see the homework problems of the previous section.
Extended Mean Value Theorem. If
f:[a,b]
R has two continuous derivatives on [a,b],
then there is a c between a and x such
- f(x) = f(a) + f'(a) (x - a) + f''(c)(x - a)2 /2
Here x is in (a,b) and c will depend on the choice of x.
We first analyze the trapezoid error on a single interval
Next we use the mean value theorem for x = a + h and (1) with c replaced by C.
Now apply the mean value theorem with f(x) replaced by f'(x) to get
Let M2 = max |f''(x)| where x is in
.(3)
Trapezoid Sum Error Theorem. If
f:[a,b]
R has two continuous derivatives on [a,b],
then the error for the trapezoid sum algorithm is bounded as indicated in
(4).
Here we can replace the 8/12 in (4) by 1/12 (see Burden and Faires, Numerical Analysis).
Proof. We may apply (3) to each interval in the trapezoid sum where
Example. Consider the integral over [0,1] for
Since x is in the interval [0,1], we have
So, the following are error estimates for both the Riemann sums and the
trapezoid sums.
Consequently, if the calculation is required to be within 0.001 of the exact
integral, we may require
Since our estimates for M1 and M2 are very crude, we can
expect these values for n to be a little large. This is evident in the
following table where this is reflected by the near constant third and fourth
columns. Also, observe the nice behavior of order two convergence.
(4)
n
RS Error
(RS Error)/h
TS Error
(TS Error)/(h2 )
10
.024373889
.2437389
.000626087
.0626087
20
.012343705
.2468741
.000156342
.0625371
40
.006211161
.2484465
.000039100
.0625610
80
.003115177
.2492142
.000009715
.0621795
160
.001560032
.2496052
.000002562
.0656127