{VERSION 6 0 "IBM INTEL NT" "6.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 2 
1 2 0 0 0 1 }{CSTYLE "_cstyle10" -1 204 "Times" 0 1 0 0 0 0 0 0 0 2 2 
2 0 0 0 1 }{CSTYLE "_cstyle11" -1 215 "Times" 1 18 0 0 0 1 2 1 1 2 2 
2 0 0 0 1 }{CSTYLE "_cstyle12" -1 216 "Times" 1 14 0 0 0 1 2 2 2 2 2 
2 0 0 0 1 }{CSTYLE "_cstyle13" -1 217 "Times" 1 14 0 0 0 1 2 2 2 2 2 
2 0 0 0 1 }{CSTYLE "_cstyle14" -1 218 "Times" 1 14 0 0 0 1 2 2 2 2 2 
2 0 0 0 1 }{CSTYLE "_cstyle15" -1 219 "Courier" 1 14 255 0 0 1 2 1 2 
2 1 2 0 0 0 1 }{CSTYLE "_cstyle16" -1 220 "Courier" 0 1 255 0 0 1 0 1 
0 2 1 2 0 0 0 1 }{CSTYLE "_cstyle19" -1 223 "Courier" 1 14 255 0 0 1 
2 1 2 2 1 2 0 0 0 1 }{CSTYLE "_cstyle20" -1 224 "Courier" 0 1 255 0 0 
1 0 1 0 2 1 2 0 0 0 1 }{CSTYLE "_cstyle21" -1 225 "Courier" 1 14 255 
0 0 1 0 1 0 2 1 2 0 0 0 1 }{CSTYLE "_cstyle22" -1 226 "Courier" 1 14 
255 0 0 1 0 1 0 2 1 2 0 0 0 1 }{CSTYLE "_cstyle23" -1 227 "Times" 1 
14 0 0 0 1 2 2 2 2 2 2 0 0 0 1 }{CSTYLE "_cstyle24" -1 228 "Times" 1 
12 0 0 0 1 2 2 2 2 2 2 0 0 0 1 }{CSTYLE "_cstyle25" -1 229 "Times" 0 
1 0 0 0 0 0 0 0 2 2 2 0 0 0 1 }{PSTYLE "_pstyle10" -1 210 1 {CSTYLE "
" -1 -1 "Times" 1 18 0 0 0 1 2 1 1 2 2 2 1 0 0 1 }3 1 0 0 12 12 2 0 2 
0 2 2 -1 1 }{PSTYLE "_pstyle11" -1 211 1 {CSTYLE "" -1 -1 "Times" 1 
14 0 0 0 1 2 2 2 2 2 2 1 0 0 1 }3 1 0 0 8 8 2 0 2 0 2 2 -1 1 }{PSTYLE 
"_pstyle12" -1 212 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 0 1 2 2 2 2 2 
2 1 0 0 1 }3 1 0 0 0 0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle13" -1 213 1 
{CSTYLE "" -1 -1 "Times" 1 14 0 0 0 1 2 2 2 2 2 2 1 0 0 1 }1 1 0 0 0 
0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle14" -1 214 1 {CSTYLE "" -1 -1 "Cou
rier" 1 14 255 0 0 1 2 1 2 2 1 2 1 0 0 1 }1 1 0 0 0 0 2 0 2 0 2 2 -1 
1 }{PSTYLE "_pstyle16" -1 216 1 {CSTYLE "" -1 -1 "Courier" 0 1 255 0 
0 1 0 1 0 2 1 2 1 0 0 1 }1 1 0 0 0 0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyl
e19" -1 219 1 {CSTYLE "" -1 -1 "Courier" 0 1 255 0 0 1 0 1 0 2 1 2 1 
0 0 1 }1 1 0 0 0 0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle20" -1 220 1 
{CSTYLE "" -1 -1 "Courier" 0 1 255 0 0 1 0 1 0 2 1 2 0 0 0 1 }1 1 0 0 
0 0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle21" -1 221 1 {CSTYLE "" -1 -1 "T
imes" 1 12 0 0 0 1 2 2 2 2 2 2 0 0 0 1 }1 1 0 0 0 0 2 0 2 0 2 2 -1 1 }
{PSTYLE "_pstyle22" -1 222 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 0 1 2 
2 2 2 2 2 1 0 0 1 }1 1 0 0 0 0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle23" 
-1 223 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 0 0 1 }1 
1 0 0 0 0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle24" -1 224 1 {CSTYLE "" 
-1 -1 "Times" 0 1 0 0 0 0 0 0 0 2 2 2 0 0 0 1 }0 0 0 -1 -1 -1 2 0 2 0 
2 2 -1 1 }{PSTYLE "_pstyle25" -1 225 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 
0 0 0 0 2 2 2 0 0 0 1 }0 0 0 -1 -1 -1 2 0 2 0 2 2 -1 1 }{PSTYLE "_psty
le26" -1 226 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 2 2 2 0 0 0 1 }0 
0 0 -1 -1 -1 1 0 1 0 2 2 -1 1 }}
{SECT 0 {EXCHG {PARA 210 "" 0 "" {TEXT 215 48 "One Dimensional Closed \+
Adaptive Trapezoidal Rule" }{TEXT 215 0 "" }}{PARA 211 "" 0 "" {TEXT 
216 12 "Wei-Chi Yang" }{TEXT 216 0 "" }}{PARA 212 "" 0 "" {TEXT 217 
23 "e-mail: wyang@runet.edu" }{TEXT 217 0 "" }}{PARA 212 "" 0 "" 
{TEXT 217 32 "URL: http://www.runet.edu/~wyang" }{TEXT 217 0 "" }}
{PARA 213 "" 0 "" {TEXT 218 0 "" }}{PARA 213 "" 0 "" {TEXT 218 117 "We
 would like to use an adaptive trapezoidal rule with a regular matrix \+
to approximate the integral of  f over [0,1]." }{TEXT 218 0 "" }}
{PARA 213 "" 0 "" {TEXT 218 0 "" }}{PARA 214 "> " 0 "" {MPLTEXT 1 219 
26 "f:=proc(x)  1/sqrt(x) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 
216 "> " 0 "" {MPLTEXT 1 220 18 "plot(f(x),x=0..1);" }{MPLTEXT 1 220 
0 "" }}}{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 1 219 24 "evalf(int(f(x),x
=0..1));" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 213 "" 0 "" {TEXT 218 
305 "We shall define the regular matrices amk and bmk. The symbols \"r
ight\" and \"left\" corresopnd to right end evaluation point and left \+
end evaluation point with repect to amk. Similarly, the symbols \"Righ
t\" and \"Left\" corresopnd to right end evaluation point and left end
 evaluation point with repect to bmk. " }{TEXT 218 0 "" }}{PARA 214 ">
 " 0 "" {MPLTEXT 1 219 49 "amk:=proc(a,b,m,k1) (2*(b-a)*(k1))/(m*(m+1)
) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 
220 13 "amk(0,1,4,1);" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 214 "> " 
0 "" {MPLTEXT 1 219 59 "bmk:=proc(a,b,m,k1) (6*(b-a)*(k1^2))/(m*(m+1)*
(2*m+1)) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 214 "> " 0 "" 
{MPLTEXT 1 219 56 "right:=proc(a,b,j,k1,m) a+sum(amk(a,b,m,j),j=1..k1)
 end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 213 "" 0 "" {TEXT 218 0 "
" }}{PARA 214 "> " 0 "" {MPLTEXT 1 219 56 "Right:=proc(a,b,j,k1,m) a+s
um(bmk(a,b,m,j),j=1..k1) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 
214 "> " 0 "" {MPLTEXT 1 219 58 "left:=proc(a,b,j,k1,m)  a+sum(amk(a,b
,m,j),j=0..k1-1) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 214 "> " 
0 "" {MPLTEXT 1 219 58 "Left:=proc(a,b,j,k1,m)  a+sum(bmk(a,b,m,j),j=0
..k1-1) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 213 "" 0 "" {TEXT 
218 75 "Here is our first adaptive quadrature with the unformly regula
r matrix amk." }{TEXT 218 0 "" }}{PARA 213 "" 0 "" {TEXT 218 0 "" }}
{PARA 214 "> " 0 "" {MPLTEXT 1 219 133 "try1:=proc(a,b,m) (amk(a,b,m,1
)*f(right(a,b,j,1,m)))/2+sum(amk(a,b,m,k1)*((f(right(a,b,j,k1,m))+f(le
ft(a,b,j,k1,m)))/2),k1=2..m) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG 
{PARA 213 "" 0 "" {TEXT 218 92 "Next we define the richardson extrapol
ation quadrature, \"richard\",  by starting with \"try0\"." }{TEXT 
218 0 "" }}{PARA 213 "" 0 "" {TEXT 218 0 "" }}{PARA 214 "> " 0 "" 
{MPLTEXT 1 219 18 "try0:=proc(a,b,m) " }{MPLTEXT 1 219 0 "" }}{PARA 
214 "> " 0 "" {MPLTEXT 1 219 78 "sum(amk(a,b,m,k1)*((f(right(a,b,j,k1,
m))+f(left(a,b,j,k1,m)))/2),k1=2..m) end;" }{MPLTEXT 1 219 0 "" }}}
{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 1 219 100 "richard:=proc(a,b,m) (1
/2)*amk(a,b,m,1)*f(right(a,b,j,1,m))+(1/3)*(4*try0(a,b,m)-try0(a,b,m/2
)) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 
1 219 42 "fac:=n->if n<2 then 1 else n*fac(n-1) end;" }{MPLTEXT 1 219 
0 "" }}}{EXCHG {PARA 219 "> " 0 "" {MPLTEXT 1 223 21 "FAC:=proc(n::pos
int);" }{MPLTEXT 1 223 0 "" }{MPLTEXT 1 223 35 "\nif n<2 then 1 else n
*fac(n-1) end;" }{MPLTEXT 1 223 0 "" }{MPLTEXT 1 223 11 "\n      end;
" }{MPLTEXT 1 224 0 "" }{MPLTEXT 1 224 1 "\n" }}}{EXCHG {PARA 214 "> \+
" 0 "" {MPLTEXT 1 219 7 "FAC(5);" }{MPLTEXT 1 219 0 "" }}}{EXCHG 
{PARA 214 "> " 0 "" {MPLTEXT 1 219 121 "Romberg:=proc(a,b,m,n) if n<2 \+
then try1(a,b,m) else try1(a,b,m-1)+(try1(a,b,m-1)-try1(a,b,(m-1)/2))/
(4^(n-1)-1) end; end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 216 "> " 
0 "" {MPLTEXT 1 220 27 "evalf(Romberg(0,1,400,20));" }{MPLTEXT 1 220 
0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 27 "evalf(Romberg(0,
1,430,20));" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" 
{MPLTEXT 1 220 27 "evalf(Romberg(0,1,440,20));" }{MPLTEXT 1 220 0 "" }
}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 27 "evalf(Romberg(0,1,445,
20));" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 
220 26 "evalf(Romberg(0,1,445,4));" }{MPLTEXT 1 220 0 "" }}}{EXCHG 
{PARA 216 "> " 0 "" {MPLTEXT 1 225 120 "Romberg2:=proc(a,b,m,n) if n<2
 then try1(a,b,m) else try1(a,b,m-1)+(try1(a,b,m-1)-try1(a,b,2*m-2))/(
4^(n-1)-1) end; end;" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 225 28 
"evalf(Romberg2(0,1,200,20));" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 
220 "> " 0 "" {MPLTEXT 1 226 28 "evalf(Romberg2(0,1,400,20));" }}}
{EXCHG {PARA 220 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 221 "" 0 
"" {TEXT 227 78 "Remark: Notice the difererent step sizes we used in '
Romberg' and 'Romberg2'. " }}}{EXCHG {PARA 220 "> " 0 "" {MPLTEXT 1 0 
0 "" }}}{EXCHG {PARA 220 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 
220 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 213 "" 0 "" {TEXT 218 
129 "Here is another quadrature, \"Try\", using bmk as the uniformly r
egular matrix. This is to compare two quadratures \"try\" and \"Try\".
" }{TEXT 218 0 "" }}{PARA 213 "" 0 "" {TEXT 218 0 "" }}{PARA 214 "> " 
0 "" {MPLTEXT 1 219 128 "Try:=proc(a,b,m) bmk(a,b,m,1)*f(Right(a,b,j,1
,m))+sum(bmk(a,b,m,k1)*((f(Right(a,b,j,k1,m))+f(Left(a,b,j,k1,m)))/2),
k1=2..m) end;" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 214 "> " 0 "" 
{MPLTEXT 1 219 21 "evalf(try1(0,1,300));" }{MPLTEXT 1 219 0 "" }}}
{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 1 219 21 "evalf(try1(0,1,400));" }
{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 1 219 21 "e
valf(try1(0,1,430));" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 214 "> " 0 
"" {MPLTEXT 1 219 21 "evalf(try1(0,1,530));" }{MPLTEXT 1 219 0 "" }}}
{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 1 219 49 "evalf(richard(0,1,300));
 evalf(richard(0,1,400));" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 214 ">
 " 0 "" {MPLTEXT 1 219 25 "evalf(richard(0,1,430)); " }{MPLTEXT 1 219 
0 "" }}}{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 1 219 24 "evalf(richard(0,
1,530));" }{MPLTEXT 1 219 0 "" }}}{EXCHG {PARA 213 "" 0 "" {TEXT 218 
129 "Notice that if we use bmk in our quadrature for this function, we
 can't pick too many points as we did for the quadrature \"try1\"." }
{TEXT 218 0 "" }}{PARA 214 "> " 0 "" {MPLTEXT 1 219 61 "evalf(Try(0,1,
300)); evalf(Try(0,1,400));evalf(Try(0,1,430));" }{MPLTEXT 1 219 0 "" 
}}}{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 1 219 61 "evalf(Try(0,1,200)); \+
evalf(Try(0,1,250));evalf(Try(0,1,280));" }{MPLTEXT 1 219 0 "" }}}
{EXCHG {PARA 214 "> " 0 "" {MPLTEXT 1 219 61 "evalf(Try(0,1,100)); eva
lf(Try(0,1,150));evalf(Try(0,1,180));" }{MPLTEXT 1 219 0 "" }}}{EXCHG 
{PARA 216 "> " 0 "" {MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 222 "" 0 "" 
{TEXT 217 90 "We conjecture that given a function, there is an optimal
 choice of uniform regular matrix." }{TEXT 217 0 "" }}}{EXCHG {PARA 
219 "> " 0 "" {MPLTEXT 1 224 6 "f:=x->" }{MPLTEXT 1 223 1 "1" }
{MPLTEXT 1 224 9 "/sqrt(x);" }{MPLTEXT 1 224 0 "" }}}{EXCHG {PARA 216 
"> " 0 "" {MPLTEXT 1 220 21 "int(f(x),x=0.001..1);" }{MPLTEXT 1 220 0 
"" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 15 "diff(f(x),x$2);" }
{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 18 "g
:=x->3/4/x^(5/2);" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" 
{MPLTEXT 1 220 34 "plot(g(x),x=0.001..1,y=-100..100);" }{MPLTEXT 1 
220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 13 "M:=g(0.0001)
;" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 
18 "amk(0.0001,1,4,1);" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 
0 "" {MPLTEXT 1 220 59 "err1:=proc(n) sum(7.5*10^9*amk(0.0001,1,n,k)^3
,k=1..n) end;" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" 
{MPLTEXT 1 220 59 "err2:=proc(n) sum(7.5*10^9*bmk(0.0001,1,n,k)^3,k=1.
.n) end;" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 
1 220 14 "err1(1000000);" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> \+
" 0 "" {MPLTEXT 1 220 13 "err2(100000);" }{MPLTEXT 1 220 0 "" }}}
{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 16 "h:=(1-0.0001)/n;" }
{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 223 "" 0 "" {TEXT 228 54 "the foll
owing is the error if trapezoidal rule is used" }{TEXT 228 0 "" }}}
{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 25 "ter:=(1-0.0001)*M*h^2/12
;" }{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 
34 "terr:=proc(n) 624812518.7/n^2 end;" }{MPLTEXT 1 220 0 "" }}}
{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 13 "terr(100000);" }
{MPLTEXT 1 220 0 "" }}}{EXCHG {PARA 216 "> " 0 "" {MPLTEXT 1 220 0 "" 
}}}{EXCHG {PARA 223 "" 0 "" {TEXT 228 0 "" }}}{EXCHG {PARA 223 "" 0 "
" {TEXT 228 0 "" }}}{EXCHG {PARA 213 "" 0 "" {TEXT 218 0 "" }}}{PARA 
224 "" 0 "" {TEXT 204 0 "" }}{PARA 224 "" 0 "" {TEXT 204 0 "" }}{PARA 
225 "" 0 "" {TEXT 229 0 "" }}{PARA 226 "" 0 "" {TEXT -1 0 "" }}}{MARK 
"11 1 0" 0 }{VIEWOPTS 1 1 0 3 2 1804 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 
1 }