{VERSION 6 0 "IBM INTEL NT" "6.0" }
{USTYLETAB {CSTYLE "_cstyle1" -1 206 "Times" 1 18 0 0 0 1 2 1 1 2 2 2 
0 0 0 1 }{CSTYLE "_cstyle2" -1 207 "Times" 1 14 0 0 0 1 2 2 2 2 2 2 0 
0 0 1 }{CSTYLE "_cstyle3" -1 208 "Times" 1 14 0 0 0 1 2 2 2 2 2 2 0 0 
0 1 }{CSTYLE "_cstyle4" -1 209 "Times" 1 14 0 0 0 1 2 2 2 2 2 2 0 0 0 
1 }{CSTYLE "_cstyle5" -1 210 "Courier" 1 14 255 0 0 1 2 1 2 2 1 2 0 0 
0 1 }{CSTYLE "_cstyle8" -1 213 "Times" 1 14 0 0 0 1 2 1 2 2 2 2 0 0 0 
1 }{PSTYLE "_pstyle1" -1 204 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 
2 1 1 2 2 2 1 1 1 1 }3 1 0 0 12 12 2 0 2 0 2 2 0 1 }{PSTYLE "_pstyle2
" -1 205 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }
3 1 0 0 8 8 2 0 2 0 2 2 0 1 }{PSTYLE "_pstyle3" -1 206 1 {CSTYLE "" 
-1 -1 "Times" 1 14 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 2 0 2 0 2 
2 0 1 }{PSTYLE "_pstyle4" -1 207 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 
0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 2 0 2 0 2 2 0 1 }{PSTYLE "_pstyle
5" -1 208 1 {CSTYLE "" -1 -1 "Courier" 1 12 255 0 0 1 2 1 2 2 1 2 1 1 
1 1 }1 1 0 0 0 0 2 0 2 0 2 2 0 1 }{PSTYLE "_pstyle10" -1 213 1 
{CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 
0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 204 "" 0 "" {TEXT 206 52 "Two Dimensional Adaptiv
e Closed and Open Quadratures" }}{PARA 205 "" 0 "" {TEXT 207 12 "Wei-C
hi Yang" }}{PARA 206 "" 0 "" {TEXT 208 23 "e-mail: wyang@runet.edu" }}
{PARA 206 "" 0 "" {TEXT 208 32 "URL: http://www.runet.edu/~wyang" }}
{PARA 207 "" 0 "" {TEXT 209 88 "We would like to use a quadrature to f
ind the double integral of the following function." }}{PARA 208 "> " 
0 "" {MPLTEXT 1 210 25 "f:=(x,y)-> 1/(sqrt(x*y));" }}}{EXCHG {PARA 
208 "> " 0 "" {MPLTEXT 1 210 47 "trueans:=evalf(int(int(f(x,y),x=0..1)
,y=0..1));" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 40 "plot3d(f(
x,y),x=0..1,y=0..1,axes=boxed);" }}}{EXCHG {PARA 207 "" 0 "" {TEXT 
209 45 "Now we define two uniformly regular matrices." }}}{EXCHG 
{PARA 208 "> " 0 "" {MPLTEXT 1 210 57 "ank:=proc(a,b,n,k) (6*(b-a)*(k^
2))/(n*(n+1)*(2*n+1)) end;" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 
210 57 "bnk:=proc(c,d,n,k) (6*(d-c)*(k^2))/(n*(n+1)*(2*n+1)) end;" }}}
{EXCHG {PARA 207 "" 0 "" {TEXT 209 166 "We define the right end and le
ft end evaluation points. The rx (below) is for right end point and r2
x is used for Romberg quadrature later; same for lx and l2x below." }}
{PARA 208 "> " 0 "" {MPLTEXT 1 210 54 "rx:=proc(a,b,j,k,n) a + sum(ank
(a,b,n,j), j=1..k) end;" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 
55 "#r2x:=proc(a,b,j,k,n) a+sum(2*ank(a,b,n,j),j=1..k) end;" }}}
{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 54 "lx:=proc(a,b,j,k,n) a+su
m(ank(a,b,n,j), j=0..k-1) end;" }}}{EXCHG {PARA 208 "> " 0 "" 
{MPLTEXT 1 210 58 "#l2x:=proc(a,b,j,k,n)  a+sum(2*ank(a,b,n,j),j=0..k-
1) end;" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 51 "ry:=proc(c,d
,j,l,n) c+sum(bnk(c,d,n,j),j=1..l) end;" }}}{EXCHG {PARA 208 "> " 0 "
" {MPLTEXT 1 210 55 "#r2y:=proc(c,d,j,l,n) c+sum(2*bnk(c,d,n,j),j=1..l
) end;" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 53 "ly:=proc(c,d,
j,l,n) c+sum(bnk(c,d,n,j),j=0..l-1) end;" }}}{EXCHG {PARA 208 "> " 0 "
" {MPLTEXT 1 210 58 "#l2y:=proc(c,d,j,l,n)  c+sum(2*bnk(c,d,n,j),j=0..
l-1) end;" }}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 39 "Here is our close
d adaptive quadrature." }}{PARA 208 "> " 0 "" {MPLTEXT 1 210 434 "trap
2closed:=proc(a,b,c,d,n) sum(ank(a,b,n,k)*sum(bnk(c,d,n,l)*((f(rx(a,b,
j,k,n),ry(c,d,j,l,n))+ f(lx(a,b,j,k,n), ly(c,d,j,l,n))+f(lx(a,b,j,k,n)
, ry(c,d,j,l,n)) + f(rx(a,b,j,k,n), ly(c,d,j,l,n)))/4), l = 2..n), k =
2..n)+sum(ank(a,b,n,k)*bnk(c,d,n,1)*(f(rx(a,b,j,k,n),ry(c,d,j,1,n))/4)
,k=2..n)+sum(ank(a,b,n,1)*bnk(c,d,n,l)*(f(rx(a,b,j,1,n),ry(c,d,j,l,n))
/4),l=2..n)+(1/4)*ank(a,b,n,1)*bnk(a,b,n,1)*f(rx(a,b,j,1,n),ry(a,b,j,1
,n))  end;" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 40 "evalf(tra
p2closed(0,1,0,1,100))-trueans;" }}}{EXCHG {PARA 208 "> " 0 "" 
{MPLTEXT 1 210 397 "#trap2hclosed:=proc(a,b,c,d,n) sum(2*ank(a,b,n,k)*
sum(2*bnk(c,d,n,l)*((f(r2x(a,b,j,k,n),r2y(c,d,j,l,n))+ f(l2x(a,b,j,k,n
), l2y(c,d,j,l,n))+f(l2x(a,b,j,k,n), r2y(c,d,j,l,n)) + f(r2x(a,b,j,k,n
), l2y(c,d,j,l,n)))/4), l = 2..n), k =2..n)+sum(2*ank(a,b,n,k)*2*bnk(c
,d,n,1)*(f(r2x(a,b,j,k,n),r2y(c,d,j,1,n))/2),k=1..n)+sum(2*ank(a,b,n,1
)*2*bnk(c,d,n,l)*(f(r2x(a,b,j,1,n),r2y(c,d,j,l,n))/2),l=1..n)  end;" }
}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 12 "Here is our " }{TEXT 213 4 "o
pen" }{TEXT 209 21 " adaptive quadrature." }{TEXT -1 0 "" }}{PARA 208 
"> " 0 "" {MPLTEXT 1 210 218 "trap2:=proc(a,b,c,d,n) sum(ank(a,b,n,k)*
sum(bnk(c,d,n,l)*((f(rx(a,b,j,k,n),ry(c,d,j,l,n))+ f(lx(a,b,j,k,n), ly
(c,d,j,l,n))+f(lx(a,b,j,k,n), ry(c,d,j,l,n)) + f(rx(a,b,j,k,n), ly(c,d
,j,l,n)))/4), l = 2..n), k =2..n) end;" }}{PARA 207 "" 0 "" {TEXT 209 
0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 34 "evalf(trap2(0,1,
0,1,100))-trueans;" }}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 70 "**Don't \+
try evalf(trap2(0,1,0,1,250)) below; it drained out my memory." }}}
{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 27 "#evalf(trap2(0,1,0,1,250
));" }}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 0 "" }}{PARA 207 "" 0 "" 
{TEXT 209 61 "Now we define the Richardson Extrapolation in two dimens
ions." }}{PARA 208 "> " 0 "" {MPLTEXT 1 210 93 "#richardclosed:=proc(a
,b,c,d,n)(1/3)*(4*trap2closed(a,b,c,d,n)-trap2closed(a,b,c,d,n/2)) end
;" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 69 "#evalf(richardclos
ed(0,1,0,1,25)); #evalf(richardclosed(0,1,0,1,30));" }}}{EXCHG {PARA 
208 "> " 0 "" {MPLTEXT 1 210 34 "#evalf(richardclosed(0,1,0,1,50));" }
}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 35 "#evalf(richardclosed(0
,1,0,1,100));" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 57 "#error
richard:=evalf(richardclosed(0,1,0,1,100))-trueans;" }}}{EXCHG {PARA 
207 "" 0 "" {TEXT 209 61 "Clearly, these results are much faster than \+
those from trap2." }}{PARA 207 "" 0 "" {TEXT 209 0 "" }}{PARA 207 "" 
0 "" {TEXT 209 29 "** Romberg Integration in 2D." }}{PARA 208 "> " 0 "
" {MPLTEXT 1 210 202 "#Romberg:=proc(a,b,c,d,m,n) if m<2 or n<2 then t
rap2closed(a,b,c,d,m,n) else trap2closed(a,b,c,d,m-1,n-1)+(trap2closed
(a,b,c,d,m-1,n-1)-trap2hclosed(a,b,c,d,m-1,n-1))/((4^(m-1)-1)*(4^(n-1)
-1)) end; end;" }}{PARA 207 "" 0 "" {TEXT 209 0 "" }}}{EXCHG {PARA 
208 "> " 0 "" {MPLTEXT 1 210 31 "#evalf(Romberg(0,1,0,1,50,50));" }}}
{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 33 "#evalf(Romberg(0,1,0,1,1
00,100));" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 55 "#errorRomb
erg:=evalf(Romberg(0,1,0,1,100,100))-trueans;" }}}{EXCHG {PARA 207 "" 
0 "" {TEXT 209 83 "Remark: Romberg closed type quadrature is better th
an the Richardson's closed type." }}}{EXCHG {PARA 208 "> " 0 "" 
{MPLTEXT 1 210 33 "#evalf(Romberg(0,1,0,1,150,150));" }}}{EXCHG {PARA 
208 "> " 0 "" {MPLTEXT 1 210 33 "#evalf(Romberg(0,1,0,1,200,200));" }}
}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 0 "" }}}{PARA 213 "" 0 "" 
{TEXT -1 0 "" }}}{MARK "35" 0 }{VIEWOPTS 1 1 0 3 2 1804 1 1 1 1 }
{PAGENUMBERS 0 1 2 33 1 1 }