{VERSION 6 0 "IBM INTEL NT" "6.0" }
{USTYLETAB {CSTYLE "2D Output" -1 20 "Times" 0 1 0 0 255 1 0 0 2 2 2 
2 0 0 0 1 }{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 "_cstyle6" -1 211 "Times" 1 14 0 0 255 1 2 2 2 2 2 2 0 0 
0 1 }{CSTYLE "_cstyle7" -1 212 "Times" 1 14 0 0 0 1 2 2 2 2 2 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 
}{CSTYLE "_cstyle9" -1 214 "Courier" 1 10 0 0 255 1 2 2 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 0 0 1 }3 1 0 0 12 12 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle2
" -1 205 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 "_pstyle3" -1 206 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 "_pstyle4" -1 207 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 "_pstyl
e5" -1 208 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 "_pstyle6" -1 209 1 
{CSTYLE "" -1 -1 "Times" 0 1 0 0 255 1 0 0 2 2 2 2 1 0 0 1 }3 3 0 0 0 
0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle7" -1 210 1 {CSTYLE "" -1 -1 "Time
s" 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 "_pstyle8" -1 211 1 {CSTYLE "" -1 -1 "Times" 0 1 0 0 255 1 0 
0 2 2 2 2 1 0 0 1 }1 3 0 0 0 0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle9" 
-1 212 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 2 2 2 2 2 1 1 0 0 
1 }1 1 0 0 0 0 2 0 2 0 2 2 -1 1 }{PSTYLE "_pstyle10" -1 213 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 204 "" 0 "" {TEXT 206 40 "Two Dimensional Adaptiv
e Open Quadrature" }{TEXT 206 0 "" }}{PARA 205 "" 0 "" {TEXT 207 12 "W
ei-Chi Yang" }{TEXT 207 0 "" }}{PARA 206 "" 0 "" {TEXT 208 23 "e-mail:
 wyang@runet.edu" }{TEXT 208 0 "" }}{PARA 206 "" 0 "" {TEXT 208 32 "UR
L: http://www.runet.edu/~wyang" }{TEXT 208 0 "" }}{PARA 207 "" 0 "" 
{TEXT 209 88 "We would like to use a quadrature to find the double int
egral of the following function." }{TEXT 209 0 "" }}{PARA 208 "> " 0 "
" {MPLTEXT 1 210 31 "f:=proc(x,y) 1/(sqrt(x*y)) end;" }{MPLTEXT 1 210 
0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#>%\"fGf*6$%\"xG%\"yG6\"F)F)*&
\"\"\"F+-%%sqrtG6#*&9$F+9%F+!\"\"F)F)F)" }{TEXT 211 0 "" }}}{EXCHG 
{PARA 208 "> " 0 "" {MPLTEXT 1 210 47 "trueans:=evalf(int(int(f(x,y),x
=0..1),y=0..1));" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 
20 "6#>%(trueansG$\"\"%\"\"!" }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> \+
" 0 "" {MPLTEXT 1 210 40 "plot3d(f(x,y),x=0..1,y=0..1,axes=boxed);" }
{MPLTEXT 1 210 0 "" }}{PARA 210 "" 1 "" {GLPLOT3D 620 620 620 
{PLOTDATA 3 "6&-%%GRIDG6%;$\"\"!F($\"\"\"F(F&7;7;%*undefinedGF-F-F-F-F
-F-F-F-F-F-F-F-F-F-F-F-F-F-F-F-F-F-F-F-7;F-$\"+++++C!\")$\"+vi0(p\"F1$
\"+Y1k&Q\"F1$\"+++++7F1$\"+HEJt5F1$\"+t*ezz*!\"*$\"+\\t9r!*F<$\"+u8G&[
)F<$\"+++++!)F<$\"+$Qm%*e(F<$\"+qAFOsF<$\"+DK?GpF<$\"+_BScmF<$\"+-)pUT
'F<$\"+]Lx'>'F<$\"+********fF<$\"+**\\&3#eF<$\"+YU&ol&F<$\"+5w(f]&F<$
\"+WJcm`F<$\"+k$HsB&F<$\"+!><o6&F<$\"+MfM/]F<$\"+%[z*)*[F<7;F-F2F6F:F?
FC$\"+IK?GpF<FKFOFSFWFenFin$\"+Tsy1ZF<$\"+uOdNXF<$\"+e/y\"Q%F<$\"+&oSE
C%F<$\"+Vg'f6%F<$\"+)*******RF<$\"+0TJ$*QF<$\"+!>LZz$F<$\"+)R!G.PF<$\"
+Mh8=OF<$\"+ZpgQNF<$\"+9;5kMF<7;F-F4F:FAF\\o$\"+aLx'>'F<$\"+\\U&ol&F<F
YFin$\"+_@!)=YF<Fbo$\"+UP'y<%F<Fho$\"+!pvI%QF<F^p$\"+j(3xd$F<Fdp$\"+,A
ngLF<$\"+Bj)fE$F<$\"+cw()yJF<$\"+xmQ)4$F<$\"+$y:P-$F<$\"+$y&>aHF<$\"+t
/E*)GF<$\"+CrUGGF<7;F-F6F?F\\oFOFWFinF`oFdoFhoF\\pF`pFdp$\"+w6?GLF<$\"
+,\\82KF<$\"+vmQ)4$F<$\"+)*******HF<$\"+)\\F/\"HF<Far$\"+0)))Hv#F<$\"+
s:G$o#F<$\"+#o9'=EF<$\"+'f3%eDF<$\"+nH<-DF<$\"+T(*[\\CF<7;F-F8FCFgpFW$
\"+)******z%F<Fbo$\"+@/ucSF<F\\pFaq$\"+[D6%R$F<$\"+R$fhB$F<Fhr$\"+HO$o
(HF<$\"+h'[&oGF<$\"+!H\"GrFF<F`s$\"+$\\lJg#F<$\"+E@#)HDF<$\"+W![BY#F<$
\"+*******R#F<$\"+u,;UBF<$\"+8-J)G#F<$\"+Q:,QAF<$\"+H-*3>#F<7;F-F:F\\o
FipFinFboFho$\"+(R!G.PF<Fdp$\"+Aj)fE$F<FhrF]rFar$\"+!)[Y<FF<FbsF[u$\"+
U(*[\\CF<$\"+-TNwBF<$\"+w5S4BF<$\"+Yf!yC#F<Fgu$\"+M**3Q@F<$\"+q=$*)3#F
<$\"+tc,V?F<$\"+)*******>F<7;F-F=FKFYF`oF]tFju$\"+F9dGMF<FfrF[rFet$\"+
MQ0NFF<Fbs$\"+2O)e^#F<$\"+1hOCCF<Fau$\"+PoynAF<$\"+#Rw+?#F<Fhv$\"+ER1
\"3#F<$\"+5-PG?F<$\"+Om[z>F<$\"+8^(R$>F<$\"+u\\Y\"*=F<$\"+)>S;&=F<7;F-
F?FOFinFdoF\\pFdpFfrFjrFarF`sFdsFhs$\"+?OR`BF<FiwFgu$\"+U.K@@F<$\"+?I)
z0#F<$\"+********>F<$\"+_qlY>F<$\"+&fmt*=F<Fgx$\"+m!o!4=F<$\"+tMIp<F<$
\"+130K<F<7;F-FAFSF[qFhoFaqF\\vF[rFar$\"+lmmmEF<F[u$\"+b247CF<Fdv$\"+%
y+)=AF<Fhv$\"+<6fl?F<F`y$\"+**\\GS>F<$\"+#3=c)=F<$\"+qeKN=F<$\"+\"Qa))
y\"F<$\"+@JuX<F<$\"+Idg0<F<$\"+7`6o;F<$\"+hJ*Hj\"F<7;F-FCFWFboF\\pF_tF
hrFetF`sF[uF_uFcuFgu$\"+X#R\\5#F<F_x$\"+$z\"ff>F<Fdy$\"+sfrS=F<F[[l$\"
+**H9T<F<$\"+ui0(p\"F<$\"+Tt:c;F<$\"+p'z!=;F<$\"+NJ^#e\"F<$\"+PL>\\:F<
7;F-FEFenF]qF`pFatF]rFcwFdsF_zFcu$\"+\"====#F<Fjv$\"+B3)p+#F<Fcx$\"+lu
Ro=F<Ffy$\"+<Q0b<F<F_[l$\"+7Z6g;F<Fb\\l$\"+nS3z:F<$\"+J%yFa\"F<$\"+rq'
)3:F<$\"+!*y4x9F<7;F-FGFinFhoFdpFhrFarFbsFhsFdvFguFjvF`y$\"+Wy`@>F<Fgx
F[[lFjy$\"++hL!o\"F<Fc[l$\"+F)Q%*e\"F<Ff\\l$\"+\"*y&=^\"F<$\"+\"*y4x9F
<$\"+O-jW9F<$\"+hN@99F<7;F-FIF^oF_qFdrFctF^vFewFjxFazFf[lF[]lF\\^l$\"+
XQ:Y=F<$\"+f$)**y<F<$\"+(pv'=<F<$\"+)e+Tm\"F<$\"+0ZT9;F<$\"+!3H*o:F<$
\"+IM3F:F<$\"+9oT)[\"F<$\"+xga_9F<$\"+](\\\">9F<$\"+Oe&zQ\"F<$\"+SCte8
F<7;F-FKF`oF^pFfrFetFbsFgwFiwFhvF_xFcxFgxF]_l$\"+8dG9<F<F`\\l$\"+]uc.;
F<$\"+N*)ob:F<Fb^l$\"+9W`r9F<$\"+IVFM9F<$\"+B%3(*R\"F<$\"+<p_n8F<$\"+4
xYP8F<$\"+StI48F<7;F-FMFboFaqFhrFgtF[uFauFguFczFh[lF]]lF[[lF__lF`\\l$
\"+*******f\"F<Ff\\l$\"+)HQH]\"F<$\"+'[$fg9F<$\"+1tj@9F<$\"+X1k&Q\"F<$
\"+2oC_8F<$\"+<l:@8F<$\"+kl6#H\"F<$\"+j5\"\\E\"F<7;F-FOFdoFdpFjrF`sF`v
FiwF\\yF`yFdyFfyFjyFa_lFf`lFf\\l$\"+*******\\\"F<$\"+\\P@b9F<Fh^l$\"+-
W\\w8F<$\"+'yS;M\"F<$\"+TtI48F<$\"+(H/#z7F<$\"+%['3^7F<$\"+q[uC7F<7;F-
FQFfoFcqF\\sFitFbvF[xF^yFezFj[lF_]lF^^lFc_lFh`lFialF\\cl$\"+0Zw69F<$\"
+!o))>P\"F<$\"+@dRN8F<$\"+YFe,8F<$\"+(p9-F\"F<$\"+b/,T7F<$\"+p@t87F<$
\"+^q<)=\"F<7;F-FSFhoFeqFarF[uFdvFhvF`yFgzF[[lF_[lFc[lFe_lFb^lF[blFh^l
F]dl$\"+KLLL8F<$\"+o8x(H\"F<Fgbl$\"+*zEWB\"F<$\"+y`/17F<$\"+[c`z6F<$\"
+P0qa6F<7;F-FUFjoFgqF^sF]uFfvF]xFbyFizF\\\\lFa]lF`^lFg_lFj`lF]blF^clF_
dlF^el$\"+%*y:j7F<$\"+AS<J7F<$\"+=G],7F<$\"+q$yQ<\"F<$\"+&pv![6F<$\"+t
H!R7\"F<7;F-FWF\\pFiqF`sF_uFguF_xFdyF[[lF^\\lFb\\lFf\\lFi_lF\\alF_blF`
clFadlFgblF[fl$\"+*******>\"F<$\"+(3!3r6F<$\"+1^:W6F<$\"+pd+>6F<$\"+9^
W&4\"F<7;F-FYF^pF[rFbsFauFhvFaxFgxF][lF`\\lFc]lFb^lF[`lF^alFablFbclFcd
lF`elF]flFhfl$\"+Ur&G9\"F<$\"+_5e;6F<$\"+7y.#4\"F<$\"+n\\/p5F<7;F-FenF
`pF]rFdsFcuFjvFcxFfyF_[lFb\\lFe]lFd^lF]`lF`alFcblFdclFedlFbelF_flFjflF
cgl$\"+!4444\"F<$\"+Q,$p1\"F<$\"+NfYW5F<7;F-FgnFbpF_rFfsFeuF\\wFexFhyF
a[lFd\\lFg]lFf^lF_`lFbalFeblFfclFgdlFdelFaflF\\glFeglF\\hl$\"+g#yM/\"F
<$\"+Oy]@5F<7;F-FinFdpFarFhsFguF^wFgxFjyFc[lFf\\lFi]lFh^lFa`lFdalFgblF
hclFidlFfelFcflF^glFgglF^hlFchl$\"+!*********!#5-%*AXESSTYLEG6#%$BOXG-
%+AXESLABELSG6%Q\"x6\"Q\"yFailQ!Fail-%*LINESTYLEG6#F(" 1 2 5 0 10 1 2 
1 1 2 2 1.000000 59.000000 -98.000000 1 0 "Curve 1" }}{TEXT 212 0 "" }
}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 45 "Now we define two uniformly r
egular matrices." }{TEXT 209 0 "" }}}{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;" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#>%$ank
Gf*6&%\"aG%\"bG%\"nG%\"kG6\"F+F+,$*.\"\"'\"\"\",&9%F/9$!\"\"F/9'\"\"#9
&F3,&F6F/F/F/F3,&*&F5F/F6F/F/F/F/F3F/F+F+F+" }{TEXT 211 0 "" }}}
{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;" }{MPLTEXT 1 210 0 "" }}{PARA 209 "
" 1 "" {XPPMATH 20 "6#>%$bnkGf*6&%\"cG%\"dG%\"nG%\"kG6\"F+F+,$*.\"\"'
\"\"\",&9%F/9$!\"\"F/9'\"\"#9&F3,&F6F/F/F/F3,&*&F5F/F6F/F/F/F/F3F/F+F+
F+" }{TEXT 211 0 "" }}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 55 "We defin
e the right end and left end evaluation points." }{TEXT 209 111 " The \+
rx (below) is for right end point and r2x is used for Romberg quadratu
re 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;" }
{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#>%#rxGf*6'%\"aG
%\"bG%\"jG%\"kG%\"nG6\"F,F,,&9$\"\"\"-%$sumG6$-%$ankG6&F.9%9(9&/F8;F/9
'F/F,F,F," }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 
210 54 "r2x:=proc(a,b,j,k,n) a+sum(2*ank(a,b,n,j),j=1..k) end;" }
{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#>%$r2xGf*6'%\"a
G%\"bG%\"jG%\"kG%\"nG6\"F,F,,&9$\"\"\"-%$sumG6$,$*&\"\"#F/-%$ankG6&F.9
%9(9&F/F//F;;F/9'F/F,F,F," }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 
"" {MPLTEXT 1 210 54 "lx:=proc(a,b,j,k,n) a+sum(ank(a,b,n,j), j=0..k-1
) end;" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#>%#lxG
f*6'%\"aG%\"bG%\"jG%\"kG%\"nG6\"F,F,,&9$\"\"\"-%$sumG6$-%$ankG6&F.9%9(
9&/F8;\"\"!,&9'F/F/!\"\"F/F,F,F," }{TEXT 211 0 "" }}}{EXCHG {PARA 208 
"> " 0 "" {MPLTEXT 1 210 57 "l2x:=proc(a,b,j,k,n)  a+sum(2*ank(a,b,n,j
),j=0..k-1) end;" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 
20 "6#>%$l2xGf*6'%\"aG%\"bG%\"jG%\"kG%\"nG6\"F,F,,&9$\"\"\"-%$sumG6$,$
*&\"\"#F/-%$ankG6&F.9%9(9&F/F//F;;\"\"!,&9'F/F/!\"\"F/F,F,F," }{TEXT 
211 0 "" }}}{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;" }{MPLTEXT 1 210 0 "" }}{PARA 
209 "" 1 "" {XPPMATH 20 "6#>%#ryGf*6'%\"cG%\"dG%\"jG%\"lG%\"nG6\"F,F,,
&9$\"\"\"-%$sumG6$-%$bnkG6&F.9%9(9&/F8;F/9'F/F,F,F," }{TEXT 211 0 "" }
}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 54 "r2y:=proc(c,d,j,l,n) c
+sum(2*bnk(c,d,n,j),j=1..l) end;" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 
1 "" {XPPMATH 20 "6#>%$r2yGf*6'%\"cG%\"dG%\"jG%\"lG%\"nG6\"F,F,,&9$\"
\"\"-%$sumG6$,$*&\"\"#F/-%$bnkG6&F.9%9(9&F/F//F;;F/9'F/F,F,F," }{TEXT 
211 0 "" }}}{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;" }{MPLTEXT 1 210 0 "" }}
{PARA 209 "" 1 "" {XPPMATH 20 "6#>%#lyGf*6'%\"cG%\"dG%\"jG%\"lG%\"nG6
\"F,F,,&9$\"\"\"-%$sumG6$-%$bnkG6&F.9%9(9&/F8;\"\"!,&9'F/F/!\"\"F/F,F,
F," }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 57 "l
2y:=proc(c,d,j,l,n)  c+sum(2*bnk(c,d,n,j),j=0..l-1) end;" }{MPLTEXT 1 
210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#>%$l2yGf*6'%\"cG%\"dG%\"jG
%\"lG%\"nG6\"F,F,,&9$\"\"\"-%$sumG6$,$*&\"\"#F/-%$bnkG6&F.9%9(9&F/F//F
;;\"\"!,&9'F/F/!\"\"F/F,F,F," }{TEXT 211 0 "" }}}{EXCHG {PARA 207 "" 
0 "" {TEXT 209 39 "Here is our closed adaptive quadrature." }{TEXT 
209 0 "" }}{PARA 208 "> " 0 "" {MPLTEXT 1 210 371 "trap2closed:=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))/2),k=1..n)+sum(an
k(a,b,n,1)*bnk(c,d,n,l)*(f(rx(a,b,j,1,n),ry(c,d,j,l,n))/2),l=1..n)  en
d;" }{MPLTEXT 1 210 0 "" }}{PARA 211 "" 1 "" {XPPMATH 20 "6#>%,trap2cl
osedGf*6'%\"aG%\"bG%\"cG%\"dG%\"nG6\"F,F,,(-%$sumG6$*&-%$ankG6&9$9%9(%
\"kG\"\"\"-F/6$*&-%$bnkG6&9&9'F7%\"lGF9,**&#F9\"\"%F9-%\"fG6$-%#rxG6'F
5F6%\"jGF8F7-%#ryG6'F@FAFMFBF7F9F9*&FEF9-FH6$-%#lxGFL-%#lyGFPF9F9*&FEF
9-FH6$FTFNF9F9*&FEF9-FH6$FJFVF9F9F9/FB;\"\"#F7F9/F8FinF9-F/6$,$*&#F9Fj
nF9*(F2F9-F>6&F@FAF7F9F9-FH6$FJ-FO6'F@FAFMF9F7F9F9F9/F8;F9F7F9-F/6$,$*
&F`oF9*(-F36&F5F6F7F9F9F=F9-FH6$-FK6'F5F6FMF9F7FNF9F9F9/FBFioF9F,F,F,
" }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 32 "eva
lf(trap2closed(0,1,0,1,100));" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "
" {XPPMATH 20 "6#$\"+kw\")**R!\"*" }{TEXT 211 0 "" }}}{EXCHG {PARA 
208 "> " 0 "" {MPLTEXT 1 210 396 "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;" }{MPLTEXT 1 210 0 "" }}{PARA 211 "" 1 "" {XPPMATH 20 "6#>%
-trap2hclosedGf*6'%\"aG%\"bG%\"cG%\"dG%\"nG6\"F,F,,(-%$sumG6$,$*(\"\"#
\"\"\"-%$ankG6&9$9%9(%\"kGF4-F/6$,$*(F3F4-%$bnkG6&9&9'F:%\"lGF4,**&#F4
\"\"%F4-%\"fG6$-%$r2xG6'F8F9%\"jGF;F:-%$r2yG6'FCFDFPFEF:F4F4*&FHF4-FK6
$-%$l2xGFO-%$l2yGFSF4F4*&FHF4-FK6$FWFQF4F4*&FHF4-FK6$FMFYF4F4F4F4/FE;F
3F:F4F4/F;F\\oF4-F/6$,$**F3F4F5F4-FA6&FCFDF:F4F4-FK6$FM-FR6'FCFDFPF4F:
F4F4/F;;F4F:F4-F/6$,$**F3F4-F66&F8F9F:F4F4F@F4-FK6$-FN6'F8F9FPF4F:FQF4
F4/FEFioF4F,F,F," }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" 
{MPLTEXT 1 210 31 "evalf(trap2closed(0,1,0,1,30));" }{MPLTEXT 1 210 0 
"" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#$\"+WRl\"*R!\"*" }{TEXT 211 0 "
" }}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 12 "Here is our " }{TEXT 213 
4 "open" }{TEXT 209 21 " adaptive quadrature." }{TEXT 209 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;" }{MPLTEXT 1 210 0 "
" }}{PARA 207 "" 0 "" {TEXT 209 0 "" }}{PARA 211 "" 1 "" {XPPMATH 20 "
6#>%&trap2Gf*6'%\"aG%\"bG%\"cG%\"dG%\"nG6\"F,F,-%$sumG6$*&-%$ankG6&9$9
%9(%\"kG\"\"\"-F.6$*&-%$bnkG6&9&9'F6%\"lGF8,**&#F8\"\"%F8-%\"fG6$-%#rx
G6'F4F5%\"jGF7F6-%#ryG6'F?F@FLFAF6F8F8*&FDF8-FG6$-%#lxGFK-%#lyGFOF8F8*
&FDF8-FG6$FSFMF8F8*&FDF8-FG6$FIFUF8F8F8/FA;\"\"#F6F8/F7FhnF,F,F," }
{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 25 "evalf(
trap2(0,1,0,1,30));" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" 
{XPPMATH 20 "6#$\"+^eEsR!\"*" }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> \+
" 0 "" {MPLTEXT 1 210 25 "evalf(trap2(0,1,0,1,80));" }{MPLTEXT 1 210 
0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#$\"+#yC[*R!\"*" }{TEXT 211 0 "
" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 26 "evalf(trap2(0,1,0,1
,100));" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#$\"+O
kW'*R!\"*" }{TEXT 211 0 "" }}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 70 "*
*Don't try evalf(trap2(0,1,0,1,250)) below; it drained out my memory.
" }{TEXT 209 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 26 "eva
lf(trap2(0,1,0,1,250));" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" 
{XPPMATH 20 "6#$\"+%H@#**R!\"*" }{TEXT 211 0 "" }}}{EXCHG {PARA 207 "
" 0 "" {TEXT 209 0 "" }}{PARA 207 "" 0 "" {TEXT 209 61 "Now we define \+
the Richardson Extrapolation in two dimensions." }{TEXT 209 0 "" }}
{PARA 208 "> " 0 "" {MPLTEXT 1 210 92 "richardclosed:=proc(a,b,c,d,n)(
1/3)*(4*trap2closed(a,b,c,d,n)-trap2closed(a,b,c,d,n/2)) end;" }
{MPLTEXT 1 210 0 "" }}{PARA 211 "" 1 "" {XPPMATH 20 "6#>%.richardclose
dGf*6'%\"aG%\"bG%\"cG%\"dG%\"nG6\"F,F,,&*&#\"\"%\"\"$\"\"\"-%,trap2clo
sedG6'9$9%9&9'9(F2F2*&#F2F1F2-F46'F6F7F8F9,$*&#F2\"\"#F2F:F2F2F2!\"\"F
,F,F," }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 
67 "evalf(richardclosed(0,1,0,1,25)); evalf(richardclosed(0,1,0,1,30))
;" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#$\"+n%o2:%!
\"*" }{TEXT 211 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#$\"+8$\\E+%!\"
*" }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 33 "ev
alf(richardclosed(0,1,0,1,50));" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 
1 "" {XPPMATH 20 "6#$\"+C#f8+%!\"*" }{TEXT 211 0 "" }}}{EXCHG {PARA 
208 "> " 0 "" {MPLTEXT 1 210 34 "evalf(richardclosed(0,1,0,1,100));" }
{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#$\"+`X\\+S!\"*
" }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 56 "err
orrichard:=evalf(richardclosed(0,1,0,1,100))-trueans;" }{MPLTEXT 1 
210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#>%-errorrichardG$\"'aX\\!
\"*" }{TEXT 211 0 "" }}}{EXCHG {PARA 207 "" 0 "" {TEXT 209 61 "Clearly
, these results are much faster than those from trap2." }{TEXT 209 0 "
" }}{PARA 207 "" 0 "" {TEXT 209 0 "" }}{PARA 207 "" 0 "" {TEXT 209 29 
"** Romberg Integration in 2D." }{TEXT 209 0 "" }}{PARA 208 "> " 0 "" 
{MPLTEXT 1 210 201 "Romberg:=proc(a,b,c,d,m,n) if m<2 or n<2 then trap
2closed(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;" }{MPLTEXT 1 210 0 "" }}{PARA 207 "" 0 "" {TEXT 209 0 "" }
}{PARA 211 "" 1 "" {XPPMATH 20 "6#>%(RombergGf*6(%\"aG%\"bG%\"cG%\"dG%
\"mG%\"nG6\"F-F-@%529(\"\"#29)F2-%,trap2closedG6(9$9%9&9'F1F4,&-F66(F8
F9F:F;,&F1\"\"\"F@!\"\",&F4F@F@FAF@*(,&F=F@-%-trap2hclosedGF>FAF@,&)\"
\"%F?F@F@FAFA,&)FIFBF@F@FAFAF@F-F-F-" }{TEXT 211 0 "" }}}{EXCHG {PARA 
208 "> " 0 "" {MPLTEXT 1 210 30 "evalf(Romberg(0,1,0,1,50,50));" }
{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#$\"+;'ew*R!\"*
" }{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 32 "eva
lf(Romberg(0,1,0,1,100,100));" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "
" {XPPMATH 20 "6#$\"+;v!)**R!\"*" }{TEXT 211 0 "" }}}{EXCHG {PARA 208 
"> " 0 "" {MPLTEXT 1 210 54 "errorRomberg:=evalf(Romberg(0,1,0,1,100,1
00))-trueans;" }{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6
#>%-errorRombergG$!'$[#>!\"*" }{TEXT 211 0 "" }}}{EXCHG {PARA 207 "" 
0 "" {TEXT 209 83 "Remark: Romberg closed type quadrature is better th
an the Richardson's closed type." }{TEXT 209 0 "" }}}{EXCHG {PARA 208 
"> " 0 "" {MPLTEXT 1 210 32 "evalf(Romberg(0,1,0,1,150,150));" }
{MPLTEXT 1 210 0 "" }}{PARA 209 "" 1 "" {XPPMATH 20 "6#$\"+kM/+S!\"*" 
}{TEXT 211 0 "" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 32 "evalf
(Romberg(0,1,0,1,200,200));" }{MPLTEXT 1 210 0 "" }}{PARA 212 "" 1 "" 
{TEXT 214 32 "Warning, computation interrupted" }{TEXT 214 0 "" }
{TEXT 214 1 "\n" }}}{EXCHG {PARA 208 "> " 0 "" {MPLTEXT 1 210 0 "" }}}
{PARA 213 "" 0 "" {TEXT -1 0 "" }}}{MARK "20 0 1" 0 }{VIEWOPTS 1 1 0 
3 2 1804 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }