{VERSION 2 3 "DEC ALPHA UNIX" "2.3" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 
1 0 0 0 0 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 
0 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }
{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 
0 0 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 1" 0 3 1 
{CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 }1 0 0 0 6 6 0 0 0 
0 0 0 -1 0 }{PSTYLE "Heading 2" 3 4 1 {CSTYLE "" -1 -1 "" 1 14 0 0 0 
0 0 0 0 0 0 0 0 0 0 }0 0 0 -1 4 4 0 0 0 0 0 0 -1 0 }{PSTYLE "Title" 0 
18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 }3 0 0 -1 12 
12 0 0 0 0 0 0 19 0 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 
0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }}
{SECT 0 {PARA 18 "" 0 "" {TEXT -1 21 "Numerical Integration" }}{PARA 
19 "" 0 "" {TEXT -1 10 "Bill Farr " }}{PARA 19 "" 0 "" {TEXT -1 13 "Ja
nuary, 2000" }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 48 "Loading the studen
t package and a simple example" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 211 
"The commands for the trapezoidal rule and Simpson's rule are in the s
tudent package. This package is part of the standard Maple distributio
n, so if you have your own copy of Maple, you already have this packag
e." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 44 "The
 command below loads the student package." }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 14 "with(student);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
68 "For an example, here is a function that is fairly easy to integrat
e." }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "f := x -> x^2*
exp(x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 55 "This computes the inte
gral of our function from 0 to 2." }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 17 "int(f(x),x=0..2);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT 
-1 57 "Using the evalf command provides a decimal approximation." }
{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "evalf(int(f(x),x=0.
.2));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 306 "The command for using t
he trapezoidal rule is trapezoid. The syntax is very similar to that o
f the int command. The last argument specifies the number of subinterv
als to use. In the command below, the number of subintervals is set to
 10, but you should experiment with increasing or decreasing this numb
er." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 72 "No
te that Maple writes out the sum and doesn't evaluate it to a number. \+
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "trapezoid(f(x),x=0..2,10);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 94 "Putting an evalf command on the ou
tside lets you calculate the trapezoidal rule approximation." }{TEXT 
-1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "evalf(trapezoid(f(x),x=0.
.2,10));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 47 "The command for Simps
on's rule is very similar." }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 24 "simpson(f(x),x=0..2,10);" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 31 "evalf(simpson(f(x),x=0..2,10));" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 28 "Two proofs of Simpson's r
ule" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 146 "In this section we use two
 different methods to derive the formula for Simpson's rule over a pai
r of subintervals.  The first subinterval is from " }{XPPEDIT 18 0 "x=
m-h" "/%\"xG,&%\"mG\"\"\"%\"hG!\"\"" }{TEXT -1 4 " to " }{XPPEDIT 18 
0 "x=m" "/%\"xG%\"mG" }{TEXT -1 24 " and the second is from " }
{XPPEDIT 18 0 "x=m" "/%\"xG%\"mG" }{TEXT -1 4 " to " }{XPPEDIT 18 0 "x
=m+h" "/%\"xG,&%\"mG\"\"\"%\"hGF&" }{TEXT -1 49 ". We want to fit a pa
rabola to the three points (" }{XPPEDIT 18 0 "(m-h,y0)" "6$,&%\"mG\"\"
\"%\"hG!\"\"%#y0G" }{TEXT -1 5 " ), (" }{XPPEDIT 18 0 "m,y1" "6$%\"mG%
#y1G" }{TEXT -1 8 "), and (" }{XPPEDIT 18 0 "m+h,y2" "6$,&%\"mG\"\"\"%
\"hGF%%#y2G" }{TEXT -1 87 "). Then, we'll integrate this parabola over
 the two subintervals to get Simpson's rule." }}}{SECT 1 {PARA 4 "" 0 
"" {TEXT -1 33 "Method 1 - a brute force approach" }}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 54 "For this method, we first define a general quadrat
ic. " }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "f1 := x -> a
*x^2+b*x+c;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 110 "The three command
s below set up the three equations that require our parabola to go thr
ough the three points. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "
eq1 := y0 = f1(m-h);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "eq2
 := y1 = f1(m);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "eq3 := y
2 = f1(m+h);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 58 "This command solv
es for the coefficients in the quadratic." }{TEXT -1 0 "" }}{PARA 0 ">
 " 0 "" {MPLTEXT 1 0 37 "sol1 := solve(\{eq1,eq2,eq3\},\{a,b,c\});" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 78 "This command assigns the values o
f a, b, and c to be the ones in our solution." }{TEXT -1 0 "" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 13 "assign(sol1);" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 83 "This shows that the values of a, b, and c have been subst
ituted into the quadratic." }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 6 "f1(x);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 33 "Next, w
e integrate the quadratic." }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}{PARA 0 "" 0 "" {TEXT -1 12 "What a mess." }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 35 "int_of_f1 := int(f1(x),x=m-h..m+h);" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 52 "Using the simplify command gives the result we \+
want." }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "simplify(in
t_of_f1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 114 "Just in case we wan
t to use a, b, or c later on in the worksheet, we clear out the values
 from the assign command." }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 29 "a := 'a'; b := 'b'; c := 'c';" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 28 "Meth
od 2 - a clever approach" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 47 "The se
cond method is more clever, and shorter. " }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 34 "First, set up a general quadratic." }{TEXT -1 0 "" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "f2 := x -> A*x^2+B*x+C;" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 58 "Next, integrate this quadratic over the t
wo subintervals. " }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 
"int2 := int(f2(x),x=m-h..m+h);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
45 "Here, we use the quadratic in Simpson's rule." }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 39 "simp1 := h/3*(f2(m-h)+4*f2(m)+f2(m+h));" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 99 "This command shows that Simpson's rule gi
ves the correct formula for the integral of the quadratic." }{TEXT -1 
0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "simplify(simp1-int2);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{SECT 1 {PARA 3 "" 0 "
" {TEXT -1 11 "Exercise 2a" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
20 "f := x -> x/(1+x^3);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 129 "This
 uses Maple's built-in numerical integration routines, which are much \+
more sophisticated than the ones we have learned about." }{TEXT -1 0 "
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "ans := evalf(int(f(x),x=0..5));
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 64 "Here is the trapezoidal rule \+
approximation with 10 subintervals." }{TEXT -1 0 "" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 41 "trap := evalf(trapezoid(f(x),x=0..5,10));" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 306 "This computes the difference betw
een the integral and the trapezoidal rule approximation. You should in
crease the number of subintervals until the absolute value of this num
ber is less than 0.001. However, you might want to go through the comm
ands just below so you can use the error estimate to guide you." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "ans-trap;" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 248 "To use the error estimate, we need to find the maximum
 of the absolute value of the second derivative of our function.  The \+
easisest way to do this is to plot the second derivative and use the m
ouse to get a good estimate of the maximum or minimum." }{TEXT -1 23 "
 I got a value of 2.38." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "plot(dif
f(f(x),x,x),x=0..5);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 65 "Using the
 formula near the end of the background section, we can " }{TEXT -1 
162 "get an upper bound on the number of subintervals required. The re
sult below says that we should need fewer than 158 subintervals to ach
ieve the desired accuracy. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "n_er
ror := sqrt(5^3*2.38/(12*0.001));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
62 "Here is the Simpson's rule approximation with 10 subintervals." }
{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "simp := evalf(simps
on(f(x),x=0..5,10));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "ans-
simp;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 93 "To get an upper bound fo
r the number of subintervals required, we can use the error estimate.
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 142 "The \+
command below plots the 4th derivative. Using the mouse, I found the m
aximum of the absolute value of the 4th derivative to be about 31.3." 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "plot(diff(f(x),x,x,x,x),x=0..5);
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 192 "Using the maximum from the p
lot and the formula from the Background section, we can calculate the \+
following upper bound for the number of subintervals required to achie
ve the desired accuracy." }{TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 40 "n_error := (5^5*31.3/(180*0.001))^(1/4);" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 
"" }}}}{MARK "7 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 }