{VERSION 3 0 "IBM INTEL LINUX" "3.0" }
{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 Output" 2 20 "" 0 1 0 0 255 1 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 8 4 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 8 2 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading 3
" 4 5 1 {CSTYLE "" -1 -1 "" 1 12 0 0 0 0 1 0 0 0 0 0 0 0 0 }0 0 0 -1 
0 0 0 0 0 0 0 0 -1 0 }{PSTYLE "Maple Output" 0 11 1 {CSTYLE "" -1 -1 "
" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 3 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 4 "" 0 "" {TEXT -1 26 "MA332 Homework 3 Solutions
" }}{PARA 0 "" 0 "" {TEXT -1 71 "Fall 1998                            \+
                   Allen B. Downey" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 39 "f:=x->evalf(x^2*exp(-x)):  a:=0:  b:=2:" }}}{EXCHG 
{PARA 5 "" 0 "" {TEXT -1 16 "Basic Quadrature" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 34 "midpoint_estimate := (b-a) * f(1);" }}{PARA 11 "
" 1 "" {XPPMATH 20 "6#>%2midpoint_estimateG$\"+C))edt!#5" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "maple_midpoint := evalf(middlesum (
f(x), x=0..2, 1));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%/maple_midpoin
tG$\"+C))edt!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "trapezoi
d_estimate := f(0) + f(2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%3trape
zoid_estimateG$\"+G8T8a!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
54 "maple_trapezoid := evalf(trapezoid (f(x), x=0..2, 1));" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6#>%0maple_trapezoidG$\"+I8T8a!#5" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "simpson_estimate := (f(0) + 4*f(1) \+
+ f(2))/3;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%1simpson_estimateG$\"+
f'H&4n!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "maple_simpson \+
:= evalf(simpson (f(x), x=0..2, 2));" }}{PARA 11 "" 1 "" {XPPMATH 20 "
6#>%.maple_simpsonG$\"+e'H&4n!#5" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 34 "exact_answer := int(f(x), x=0..2);" }}{PARA 11 "" 1 "
" {XPPMATH 20 "6#>%-exact_answerG$\"*orkY'!\"*" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 51 "midpoint_error := midpoint_estimate - exact_an
swer;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%/midpoint_errorG$\"*Wr6\"*)
!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "trapezoid_error := t
rapezoid_estimate - exact_answer;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>
%0trapezoid_errorG$!+_.1`5!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 50 "simpson_error := midpoint_estimate - exact_answer;" }}{PARA 11 "
" 1 "" {XPPMATH 20 "6#>%.simpson_errorG$\"*Wr6\"*)!#5" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "fprime := D(f);" }}{PARA 11 "" 1 "
" {XPPMATH 20 "6#>%'fprimeGR6#%\"xG6\"6$%)operatorG%&arrowGF(-%&evalfG
6#,&*&9$\"\"\"-%$expG6#,$F1!\"\"F2\"\"#*&)F1F8\"\"\"F3F;F7F(F(F(" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "f2prime := D(fprime);" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(f2primeGR6#%\"xG6\"6$%)operatorG%&a
rrowGF(-%&evalfG6#,(-%$expG6#,$9$!\"\"\"\"#*&F4\"\"\"F0F8!\"%*&)F4F6\"
\"\"F0F<F8F(F(F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "f3prime
 := D(f2prime);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%(f3primeGR6#%\"xG
6\"6$%)operatorG%&arrowGF(-%&evalfG6#,(-%$expG6#,$9$!\"\"!\"'*&F4\"\"
\"F0F8\"\"'*&)F4\"\"#\"\"\"F0F=F5F(F(F(" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "f4prime := D(f3prime);" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#>%(f4primeGR6#%\"xG6\"6$%)operatorG%&arrowGF(-%&evalfG6#,(-%$exp
G6#,$9$!\"\"\"#7*&F4\"\"\"F0F8!\")*&)F4\"\"#\"\"\"F0F=F8F(F(F(" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "f2prime_max := evalf(maximiz
e (f2prime(x),\{x\},\{x=0..2\}));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>
%,f2prime_maxG$\"\"#\"\"!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
57 "f4prime_max := evalf(maximize (f4prime(x),\{x\},\{x=0..2\}));" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%,f4prime_maxG$\"#7\"\"!" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "midpoint_error_bound := f2prime_max
 / 24 * 8;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%5midpoint_error_boundG
$\"+nmmmm!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "trapezoid_e
rror_bound := f2prime_max / 12 * 8;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6
#>%6trapezoid_error_boundG$\"+LLLL8!\"*" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 47 "simpson_error_bound := f4prime_max / 2880 * 32;" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%4simpson_error_boundG$\"+LLLL8!#5" }
}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "All errors are well within the e
rror bounds." }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 43 "The answer to Que
stion 9 on page 115 is 1/2" }}{PARA 0 "" 0 "" {TEXT -1 45 "The answer \+
to Question 10 on page 115 is 13/3" }}}{EXCHG {PARA 5 "" 0 "" {TEXT 
-1 20 "Composite Quadrature" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
69 "intervals := 8:  n:= 2 * (intervals-1):  h := (b-a)/(n+2):  Tot :=
 0:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for j from 0 to n/2 \+
do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 " xj := a + (2*j+1)*h;" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 " Tot := Tot + f(xj);" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
21 "midpoint2 := 2*h*Tot;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%*midpoi
nt2G$\"+__vOl!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "midpoin
t4 := 2*h*Tot;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%*midpoint4G$\"+5%Q
6Z'!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "midpoint8 := 2*h*
Tot;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%*midpoint8G$\"+Qywmk!#5" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "intervals := 8:  n := interv
als:  h:= (b-a)/n:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "Tot:=
 (f(a) + f(b))/2:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "for j \+
from 1 to n-1 do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "  xj := a + j*h
;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "  Tot := Tot + f(xj);" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "trapezoid2 := h * Tot;" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#>%+trapezoid2G$\"+w+]&Q'!#5" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "trapezoid4 := h * Tot;" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#>%+trapezoid4G$\"+lw7hk!#5" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 22 "trapezoid8 := h * Tot;" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#>%+trapezoid8G$\"+NI8mk!#5" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 55 "maple_trapezoid2 := evalf(trapezoid (f(x), x=0..2, 2)
);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%1maple_trapezoid2G$\"+w+]&Q'!#
5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "maple_trapezoid4 := ev
alf(trapezoid (f(x), x=0..2, 4));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>
%1maple_trapezoid2G$\"+iw7hk!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 55 "maple_trapezoid8 := evalf(trapezoid (f(x), x=0..2, 8));" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%1maple_trapezoid2G$\"+OI8mk!#5" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 84 "For n intervals, the midpoint rule
 evaluates the function at n places; the trapezoid" }}{PARA 0 "" 0 "" 
{TEXT -1 79 "rule at n+1.  Nevertheless, the midpoint estimates are be
tter across the board." }}{PARA 0 "" 0 "" {TEXT -1 85 "The only price \+
to pay is that it appears to be a little more complicated to implement
" }}{PARA 0 "" 0 "" {TEXT -1 18 "the midpoint rule." }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 45 "p := x-> evalf(exp (-x^2 / 2)):  a:=0:  b
:=2:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "intervals := 256:  \+
n:= 2 * (intervals-1):  h := (b-a)/(n+2):  Tot := 0:" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 70 "for j from 0 to n/2 do \n xj := a + (2*j+
1)*h;\n Tot := Tot + p(xj);\nod:" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 24 "integral := 2 * h * Tot:" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 44 "answer := evalf(2*integral / sqrt (2 * Pi));" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'answerG$\"+/G+X&*!#5" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "exact_integral := int (p(x), x=a..b
):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "exact_answer := evalf
(2*exact_integral / sqrt (2 * Pi));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6
#>%-exact_answerG$\"+jt*\\a*!#5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 
61 "It takes about 256 subintervals to get the required accuracy." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "p := x-> evalf(exp (-x^2 / 2
)):  a:=0:  b:=3:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "interv
als := 64:  n:= 2 * (intervals-1):  h := (b-a)/(n+2):  Tot := 0:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "for j from 0 to n/2 do \n xj
 := a + (2*j+1)*h;\n Tot := Tot + p(xj);\nod:" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 24 "integral := 2 * h * Tot:" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 44 "answer := evalf(2*integral / sqrt (2 * Pi));" 
}}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%'answerG$\"+dj-t**!#5" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "exact_integral := int (p(x), x=a..b
):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "exact_answer := evalf
(2*exact_integral / sqrt (2 * Pi));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6
#>%-exact_answerG$\"+Q?+t**!#5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 82 
"Interestingly, it takes fewer subintervals to get the integral over a
 wider range." }}}{EXCHG {PARA 5 "" 0 "" {TEXT -1 19 "Romberg Integrat
ion" }}{PARA 0 "" 0 "" {TEXT -1 71 "Since the trapezoid rule has O(h^2
) error, the denominator is 2^2-1 = 3" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 27 "bad_estimate := trapezoid4:" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 28 "good_estimate := trapezoid8:" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 63 "richardson := good_estimate + (good_estimat
e - bad_estimate)/3;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%+richardsonG
$\"+\"\\,yY'!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "exact_an
swer := int(f(x), x=0..2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%-exact
_answerG$\"*orkY'!\"*" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "er
ror := richardson - exact_answer;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>
%&errorG$\"(6)H8!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "old_
error := good_estimate - exact_answer;" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#>%*old_errorG$!'X'Q$!#5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 84 "
In this case the extrapolation actually makes things worse!  That's po
ssible because" }}{PARA 0 "" 0 "" {TEXT -1 84 "the Richardson extrapol
ation is based on the assumption that the leading-order error" }}
{PARA 0 "" 0 "" {TEXT -1 89 "term is significantly bigger than the sum
 of the other error terms, and apparently that's" }}{PARA 0 "" 0 "" 
{TEXT -1 16 "not always true." }}}{EXCHG {PARA 5 "" 0 "" {TEXT -1 19 "
Gaussian quadrature" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "x1 :
= sqrt(3)/3:  x2 := -x1:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 
"t1 := x1+1:  t2 := x2+1:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
22 "est2 := f(t1) + f(t2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%est2G
$\"+%fA*3j!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "error := e
st2 - exact_answer;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%&errorG$!*'3
\\v:!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "c1 := evalf(5/9)
;  c3 := c1;  c2 := evalf(8/9);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#
c1G$\"+cbbbb!#5" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#c3G$\"+cbbbb!#5
" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#c2G$\"+*)))))))))!#5" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "x1 := 0.7745966692;  x3 := -
x1;  x2 := 0;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x1G$\"+#pmfu(!#5" 
}}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%#x3G$!+#pmfu(!#5" }}{PARA 11 "" 1 
"" {XPPMATH 20 "6#>%#x2G\"\"!" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 36 "t1 := x1+1:  t2 := x2+1:  t3:= x3+1:" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 45 "est3 := c1 * f(t1) + c2 * f(t2) + c3 * f(t3);" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%%est3G$\"+&>M<Y'!#5" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "error := est3 - exact_answer;" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%&errorG$!(&[PZ!#5" }}}{EXCHG {PARA 
0 "" 0 "" {TEXT -1 72 "The three point Gaussian  is better than either
 midpoint4 or trapezoid4." }}}}{MARK "38 0 0" 16 }{VIEWOPTS 1 1 0 1 1 
1803 }