# Lesson 18. Numerical Integration # --------------------------------- # # Maple uses some quite sophisticated methods for numerical # approximation of integrals. We'll look at some rather less # sophisticated ones, to get an idea of what is behind some of # these. # # We'll start with some that you may remember from Math 101: the # Midpoint Rule, Trapezoid Rule and Simpson's Rule. # # Suppose we want to approximate int(f(x), x=a .. b). # For the Midpoint Rule, we divide the interval [a,b] into n # equal subintervals and approximate f(x) on each subinterval by # its value at the midpoint of the subinterval. We'll usually # have a, b and f fixed, and look at different n's, so we write # everything as a function of n. > h:= n ->(b-a)/n: > X:= (i,n) -> a + i*(b-a)/n: > midrule:= n -> > h(n)*sum(f((X(i,n)+X(i+1,n))/2), i = 0 .. n-1); /n - 1 \ |----- | | \ | midrule := n -> h(n) | ) f(1/2 X(i, n) + 1/2 X(i + 1, n))| | / | |----- | \i = 0 / > a:= 0: b:= 1: midrule(3); 1/3 f(1/6) + 1/3 f(1/2) + 1/3 f(5/6) # For the Trapezoid Rule, we approximate f(x) on the subinterval # by the average of its values at the two endpoints. > traprule:= n -> > h(n) * sum((f(X(i,n))+f(X(i+1,n)))/2, i=0..n-1); /n - 1 \ |----- | | \ | traprule := n -> h(n) | ) (1/2 f(X(i, n)) + 1/2 f(X(i + 1, n)))| | / | |----- | \i = 0 / > traprule(3); 1/6 f(0) + 1/3 f(1/3) + 1/3 f(2/3) + 1/6 f(1) # Let's try a typical function f for which we know the true # value of the integral, and look at the error (the difference # between the true value and the approximation) for Midpoint and # Trapezoid Rules with different values of n. > f:= x -> 1/(x+1): J:= int(f(x), x=a..b); J := ln(2) > seq([nn, evalf(J-midrule(nn)), evalf(J - traprule(nn))], nn= 1 .. 5); [1, .0264805139, -.0568528194], [2, .0074328949, -.0151861527], [3, .0033924908, -.0068528194], [4, .0019272894, -.0038766289], [5, .0012392949, -.0024877400] # The errors decrease, of course, as n increases. It turns out # the errors in both methods are approximately proportional to # 1/n^2. > plot({[seq([nn, nn^2*evalf(J-midrule(nn))], > nn=1 .. 20)], > [seq([nn, nn^2*evalf(J-traprule(nn))], > nn=1 .. 20)]},style=point); ** Maple V Graphics ** # The theoretical result is that the error in the Midpoint Rule # is f''(t)(b-a)^3/(24n^2) for some t between a and b, and the # error in the Trapezoid Rule is -f''(t)(b-a)^3/(12n^2) (for a # different t). This assumes that f'' is continuous on the # interval [a,b]. If not, the error might not be proportional # to 1/n^2. > f:= x -> sqrt(x): J:= int(f(x), x=a..b): > plot({[seq([nn, nn^2*evalf(J-midrule(nn))], > nn=1 .. 20)], > [seq([nn, nn^2*evalf(J-traprule(nn))], > nn=1 .. 20)]},style=point); ** Maple V Graphics ** # When f'' is continuous,the Midpoint Rule's error is # approximately -1/2 times the Trapezoid Rule's error. This # might suggest that a combination of the two rules would cancel # out the error (at least approximately) and produce a much # better approximation. The appropriate combination is # 2/3*Midpoint + 1/3*Trapezoid. This is called Simpson's Rule. > f:= 'f': 2/3*midrule(3)+1/3*traprule(3); 2/9 f(1/6) + 2/9 f(1/2) + 2/9 f(5/6) + 1/18 f(0) + 1/9 f(1/3) + 1/9 f(2/3) + 1/18 f(1) # We'd call this the Simpson's Rule approximation for n=6 # (rather than n=3): it uses the same 7 equally-spaced points as # the Trapezoid Rule for n=6. We only use simpson(n) when n is # even. > simpson:= n -> > h(n) * sum(1/3*f(X(2*i,n))+4/3*f(X(2*i+1,n))+1/3*f(X(2*i+2,n)), i=0..n/2-1); simpson := n -> h(n) /1/2 n - 1 \ | ----- | | \ | | ) (1/3 f(X(2 i, n)) + 4/3 f(X(2 i + 1, n)) + 1/3 f(X(2 i + 2, n)))| | / | | ----- | \ i = 0 / > simpson(6); 2/9 f(1/6) + 2/9 f(1/2) + 2/9 f(5/6) + 1/18 f(0) + 1/9 f(1/3) + 1/9 f(2/3) + 1/18 f(1) # The error for Simpson's Rule is (D^4)f(t)(b-a)^5/(180n^4). > f:= x -> 1/(x+1): J:= ln(2): > plot([seq([2*nn, (2*nn)^4*evalf(J-simpson(2*nn))], nn=1 .. 20)],style=point); ** Maple V Graphics ** # To get an idea of how much better this is # than Midpoint or Trapezoid, let's see what n would be needed # to have an error less than 10^(-10) in absolute value. # For the Trapezoid Rule, the error was approximately -.06/n^2. # So n would have to be greater than > (.06*10^10)^(1/2); 24494.89743 # For the Midpoint Rule, the error was approximately .03/n^2. > (.03*10^10)^(1/2); 17320.50808 # For Simpson's Rule, it was approximately # -.03/n^4. Change that .03 to .031 for safety. > (.031*10^10)^(1/4); 132.6906811 # Let's try it. I'll increase Digits to # reduce roundoff error. > Digits:= 15: evalf(J-simpson(134)); Digits:= 10: -10 -.96911*10 # Now let's look at it from a different point of view. The # Midpoint and Trapezoid rules give the correct answers for # polynomials of degree 1, but in general not for higher # degrees. > poly:= n -> unapply(sum(c[j]*x^j, j=0..n), x); n ----- \ j poly := n -> unapply( ) c[j] x , x) / ----- j = 0 > f:=poly(5); 2 3 4 5 f := x -> c[0] + c[1] x + c[2] x + c[3] x + c[4] x + c[5] x > J:= f -> int(f(x), x=a..b): J(f) - midrule(1); 11 13 1/12 c[2] + 1/8 c[3] + ---- c[4] + ---- c[5] 80 96 > J(f)-traprule(1); - 1/6 c[2] - 1/4 c[3] - 3/10 c[4] - 1/3 c[5] # Simpson's Rule is exact for polynomials of degree 3. > J(f)-simpson(2); - 1/120 c[4] - 1/48 c[5] # We could have used this to derive Simpson's Rule in another # way. Suppose we didn't know the coefficients, but we knew the # general form. > rule:= n -> > h(n) * sum(d[0]*f(X(2*i,n))+d[1]*f(X(2*i+1,n))+d[2]*f(X(2*i+2,n)), i=0..n/2-1); rule := n -> h(n) /1/2 n - 1 \ | ----- | | \ | | ) (d[0] f(X(2 i, n)) + d[1] f(X(2 i + 1, n)) + d[2] f(X(2 i + 2, n)))| | / | | ----- | \ i = 0 / > J(f)-rule(2); c[0] + 1/2 c[1] + 1/3 c[2] + 1/4 c[3] + 1/5 c[4] + 1/6 c[5] - 1/2 d[0] c[0] - 1/2 d[1] (c[0] + 1/2 c[1] + 1/4 c[2] + 1/8 c[3] + 1/16 c[4] + 1/32 c[5]) - 1/2 d[2] (c[0] + c[1] + c[2] + c[3] + c[4] + c[5]) > eqns:= { seq(coeff(expand("),c[jj],1), jj=0..2) }; eqns := {1/2 - 1/4 d[1] - 1/2 d[2], 1/3 - 1/8 d[1] - 1/2 d[2], 1 - 1/2 d[0] - 1/2 d[1] - 1/2 d[2]} # I'm only doing this for j up to 2, not 3: # there are only three parameters, d[0] to d[2], so I want three # equations. The fact that it is exact for x^3 too is an added # bonus. > solve(eqns); {d[0] = 1/3, d[2] = 1/3, d[1] = 4/3} # This can be generalized. > rule:= (n,k) -> > h(n) * sum(sum(d[j]*f(X(k*i+j,n)), j=0..k), i=0..n/k-1); /n/k - 1 / k \\ | ----- |----- || | \ | \ || rule := (n,k) -> h(n) | ) | ) d[j] f(X(k i + j, n))|| | / | / || | ----- |----- || \ i = 0 \j = 0 // > f:= 'f': rule(2,2); 1/2 d[0] f(0) + 1/2 d[1] f(1/2) + 1/2 d[2] f(1) > eqns:= k -> eval({ seq(subs(f=unapply(x^ii,x), rule(k,k)=J(f)), ii=0..k) }); Warning, `ii` is implicitly declared local eqns := k -> local ii; ii eval({seq(subs(f = unapply(x , x), rule(k, k) = J(f)), ii = 0 .. k)}) > eqns(2); {1/4 d[1] + 1/2 d[2] = 1/2, 1/8 d[1] + 1/2 d[2] = 1/3, 1/2 d[0] + 1/2 d[1] + 1/2 d[2] = 1} > eqns(3); {1/27 d[1] + 4/27 d[2] + 1/3 d[3] = 1/3, 1/81 d[1] + 8/81 d[2] + 1/3 d[3] = 1/4, 1/9 d[1] + 2/9 d[2] + 1/3 d[3] = 1/2, 1/3 d[0] + 1/3 d[1] + 1/3 d[2] + 1/3 d[3] = 1} > solve("); {d[0] = 3/8, d[3] = 3/8, d[2] = 9/8, d[1] = 9/8} > subs(", rule(3,3)); 1/8 f(0) + 3/8 f(1/3) + 3/8 f(2/3) + 1/8 f(1) > eval(subs(f= poly(4),J(f)-")); - 1/270 c[4] > subs(solve(eqns(4)), rule(4,4)); 16 16 7/90 f(0) + ---- f(1/4) + 2/15 f(1/2) + ---- f(3/4) + 7/90 f(1) 45 45 # The rule for k=3 is exact for polynomials of degree 3 but not # of degree 4, just like Simpson's rule. But the k=4 rule will # be exact for polynomials of degree 5. In general, when k is # even the rule will be exact for polynomials of degree k+1. # These rules are called Newton-Cotes rules. > for kk from 2 to 12 do > ncrule[kk]:= unapply(subs(solve(eqns(kk)), rule(n,kk)), n) od: > ncrule[4](4); 16 16 7/90 f(0) + ---- f(1/4) + 2/15 f(1/2) + ---- f(3/4) + 7/90 f(1) 45 45 > eval(subs(f=poly(6),J(f)-")); - 1/2688 c[6] > eval(subs(f=poly(14),J(f)-ncrule[12](12))); 251 - ----------- c[14] 12899450880 # The error in the Newton-Cotes rule ncrule[k](n) is # approximately proportional to 1/n^k if k is odd, or 1/n^(k+1) # if k is even. Let's look at this for k=2,3 and 4 with # f(x)=1/(1+x). The k=4 rule is so accurate that I need to # increase Digits (otherwise roundoff error will overwhelm the # real error). > f:= x -> 1/(1+x): Digits:= 20: > plot({ [ > seq([2*jj,(ln(2)-ncrule[2](2*jj))*(2*jj)^4], jj=1..20)],[ > seq([3*jj,(ln(2)-ncrule[3](3*jj))*(3*jj)^4], jj=1..13)],[ > seq([4*jj,(ln(2)-ncrule[4](4*jj))*(4*jj)^6], jj=1..10)]}, style=point); ** Maple V Graphics ** # You might think that the higher the order, # the better. That is, if one method has error estimate A/n^p # and a second has B/n^q with q > p, then the second will be # more accurate when n is sufficiently large. # This isn't necessarily true for a fixed n, however, especially # for a function whose higher-order derivatives may grow # rapidly. > f:= x -> 1/(x^2 + 1/20): J(f); 1/2 1/2 2 5 arctan(2 5 ) > map(k -> evalf(sqrt(20)*arctan(sqrt(20))-ncrule[k](36)), [2,3,4,6,9,12]); -7 -6 -5 [.619560396912*10 , .2226565484415*10 , -.10063093342562*10 , -5 -5 -6 .35623850744073*10 , .11208195265217*10 , -.5031750236487*10 ] >