{VERSION 3 0 "IBM INTEL NT" "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 Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 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 }{CSTYLE "
" -1 256 "Helvetica" 1 14 128 0 0 1 0 0 2 0 0 0 0 0 0 }{CSTYLE "" -1 
257 "" 0 24 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "Helvetica
" 0 1 128 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "Helvetica" 0 1 
128 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 255 1 0 0 0 
0 0 0 0 0 0 }{CSTYLE "" -1 261 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 262 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 
263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 0 1 0 0 
255 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 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 "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 }{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 "" 18 256 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 
0 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT 257 18 "Romberg Quadrature" }}
{PARA 0 "" 0 "" {TEXT 256 19 "Mth 351 Winter 1999" }}{PARA 0 "" 0 "" 
{TEXT 258 16 "Bent E. Petersen" }}{PARA 0 "" 0 "" {TEXT 259 32 "Filena
me: romberg_quadrature.mws" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 
"" 0 "" {TEXT -1 109 "This worksheet illustrates how Romberg quadratur
e is carried out by extrapolation from the trapezoidal rule. " }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 351 "Maple's \+
student package contains procedures to carry out trapezoidal and Simps
on's quadrature. Rather than using these procedures we will define our
 own since it is not difficult to do so, and it is instructive. The pr
ocedures in the student package are designed to work with expressions \+
(in some variable). We will design ours to work with functions." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 8 "restart;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 
"" {TEXT -1 36 "First here is the trapezoidal rule T" }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "T:=proc(f
,r::range,n)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "  local k,s,h,a,b;
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "  if type(n, numeric) then" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "    if not type(n,posint) then ERRO
R(\"n must be a positive integer\"); fi;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 5 "  fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "  a:=op(1,
r); b:=op(2,r);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "  h:=(b-a)/n;" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "  s:=f(a)+f(b)+2*sum(f(a+k*h),k=1.
.n-1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "  s*h/2;" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 250 "Note we check  n  carefully if it is num
eric, but we do not check it if it is symbolic. The parameter  f  shou
ld be a function, of course, but we do not check it either, because we
 want to be able to do symbolic manipulation with unspecified symbols.
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 90 "Here'
s our Simpson's rule procedure. Note  n  must be a positive even integ
er or symbolic." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 21 "S:=proc(f,r::range,n)" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 20 "  local k,m,s,h,a,b;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 24 "if type(n, numeric) then" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 " \+
   if not type(n,posint) then ERROR(\"n must be a positive integer\");
 fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "    if irem(n,2)=1 then ERR
OR(\"n must be an even integer\"); fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 3 "fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "m:=n/2; a:=op(1,r); b
:=op(2,r);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "h:=(b-a)/n;" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 71 "s:=f(a)+f(b)+4*sum(f(a+(2*k+1)*h),k=0..m-
1)+2*sum(f(a+2*k*h),k=1..m-1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "s*
h/3;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 113 "The error in the tr
apezoid rule can be expanded in an asymptotic series in positive even \+
powers of the step size " }{XPPEDIT 18 0 "h = (b-a)/n;" "6#/%\"hG*&,&%
\"bG\"\"\"%\"aG!\"\"F(%\"nGF*" }{TEXT -1 71 ". In Richardson extrapola
tion we use two different step sizes (usually " }{XPPEDIT 18 0 "h;" "6
#%\"hG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "2*h;" "6#*&\"\"#\"\"\"%\"h
GF%" }{TEXT -1 267 ") to eliminate the top order term in an error expa
nsion. In the case of the trapezoidal rule, the result is an order 4 r
ule which turns out to be Simpson's rule. If we extrapolate again we o
btain an order 6 rule, and so on. If we start with the trapezoidal rul
e with " }{XPPEDIT 18 0 "2^m;" "6#)\"\"#%\"mG" }{TEXT -1 30 " subinter
vals and extrapolate " }{XPPEDIT 18 0 "m;" "6#%\"mG" }{TEXT -1 21 " ti
mes we obtain the " }{XPPEDIT 18 0 "m^th;" "6#)%\"mG%#thG" }{TEXT -1 
97 " Romberg estimate of the integral. The calculation of course invol
ves the trapezoidal rules with " }{XPPEDIT 18 0 "2^k;" "6#)\"\"#%\"kG
" }{TEXT -1 16 "  subintervals, " }{XPPEDIT 18 0 "k = 0;" "6#/%\"kG\"
\"!" }{TEXT -1 6 ",1,..." }{XPPEDIT 18 0 "m;" "6#%\"mG" }{TEXT -1 1 ".
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 11 "Here'
s the " }{XPPEDIT 18 0 "m^th;" "6#)%\"mG%#thG" }{TEXT -1 18 " Romberg \+
estimate:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 21 "R:=proc(f,r::range,m)" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 12 "local k,j,R;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "if
 type(m, numeric) then" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "    if no
t type(m,nonnegint) then ERROR(\"m must be a non-negative integer\"); \+
fi; fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for k from 0 to m do" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "  R[0,k]:=T(f,r,2^k); od;" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 20 "for j from 1 to m do" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 22 "  for k from j to m do" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 46 "  R[j,k]:=1/(4^j-1)*(4^j*R[j-1,k]-R[j-1,k-1]);" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "od; od;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 7 "R[m,m];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 9 "Not
e the " }{XPPEDIT 18 0 "0^th;" "6#)\"\"!%#thG" }{TEXT -1 69 " estimate
 is the trapezoidal rule with one interval (as it should be)" }}{PARA 
0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "R
(f,a..b,0);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,$*&,&-%\"fG6#%\"aG\"\"
\"-F'6#%\"bGF*F*,&F-F*F)!\"\"F*#F*\"\"#" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 152 "Note if we had typed  f  as a function in the trapezoida
l code above we would not be able to obtain this simple expression wit
hout actually supplying a " }{TEXT 263 8 "function" }{TEXT -1 8 " for \+
 f." }}{PARA 0 "" 0 "" {TEXT -1 4 "The " }{XPPEDIT 18 0 "1^st;" "6#)\"
\"\"%#stG" }{TEXT -1 69 "  estimate is Simpson's rule with two subinte
rvals, again as expected" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "simplify(R(f,a..b,1)-S(f,a..b,2));
" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#\"\"!" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 72 "Our confidence bolstered \+
by these two checks let's try some comparisons." }}}{EXCHG {PARA 0 "" 
0 "" {TEXT 260 0 "" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 350 "The
 code for  R(f, a..b, m)  is particularly susceptible to loss of signi
ficance error (due to roundoff). Maple computes  R(f, a..b, m)  exactl
y symbolically, but a sum of very many terms of roughly the same magni
tude (but varying signs)  results. Converting the large sum that resul
ts to floating point turns out to be quite a challenge for evalf()." }
}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 296 "In the \+
examples below the relative error is first computed symbolically. Then
 it is converted to floating point. Thus it is possible to control the
 roundoff error in computing the error simply by increasing the precis
ion of conversion. With high enough precision we see only the truncati
on error." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 4 "" 0 "" {TEXT 
261 9 "Example 1" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 ">
 " 0 "" {MPLTEXT 1 0 37 "rng1:=0..Pi: ex1:=Int(sin(x),x=rng1);" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$ex1G-%$IntG6$-%$sinG6#%\"xG/F+;\"\"
!%#PiG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 119 "`trapezoidal 32`
; ex:=T(sin,rng1,32):evalf(ex); rel_error = evalf((ex1-ex)/ex1,10); re
l_error = evalf((ex1-ex)/ex1,24);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%
/trapezoidal~32G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"+hLR)*>!\"*" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$\"++]>L!)!#8" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$\"9++c,nqF\\^>L!)!#F" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 117 "`Simpson's 32`; ex:=S(sin,r
ng1,32):evalf(ex); rel_error = evalf((ex1-ex)/ex1,10); rel_error = eva
lf((ex1-ex)/ex1,24);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%-Simpson's~32
G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"+M5++?!\"*" }}{PARA 11 "" 1 "
" {XPPMATH 20 "6#/%*rel_errorG$!++++q^!#;" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#/%*rel_errorG$!9+++NbM+lq%o;&!#I" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 157 "`Romberg 5 (32)`; ex:=R(sin,rng1,5):evalf(ex,
10); rel_error = evalf((ex1-ex)/ex1,10); rel_error = evalf((ex1-ex)/ex
1,24); rel_error = evalf((ex1-ex)/ex1,36);" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#%/Romberg~5~(32)G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$
\"+********>!\"*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$\"++
+lEt!#?" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!9+++!**e0)G2
A0m!#O" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!E+++Ir'H3f2jq
%yG2A0m!#[" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 120 "Only the Romberg estimates showed much sensitivity to po
ssible roundoff here. On the other hand is is far more accurate." }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 4 "" 0 "" {TEXT 262 9 "Ex
ample 2" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 37 "rng2:=8..12: ex2:=Int(exp(x),x=rng2);" }}{PARA 11 "" 
1 "" {XPPMATH 20 "6#>%$ex2G-%$IntG6$-%$expG6#%\"xG/F+;\"\")\"#7" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 117 "`trapezoidal 32`; ex:=T(exp
,rng2,32):evalf(ex); rel_error = evalf((ex2-ex)/ex2,10); rel_error=eva
lf((ex2-ex)/ex2,24);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%/trapezoidal~
32G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"+\"==)*f\"!\"%" }}{PARA 11 "
" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!+JSu,8!#7" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#/%*rel_errorG$!9e\"4.R3I_PW<I\"!#E" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 117 "`Simpson's 32`; ex:=S(exp,rng2,32):evalf(ex
); rel_error = evalf((ex2-ex)/ex2,10); rel_error = evalf((ex2-ex)/ex2,
24);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%-Simpson's~32G" }}{PARA 11 "
" 1 "" {XPPMATH 20 "6#$\"+(\\Sxf\"!\"%" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#/%*rel_errorG$!+2^%RN\"!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%
*rel_errorG$!9/0v>^)R-*z\"QN\"!#H" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 154 "`Romberg 5 (32)`; ex:=R(exp,rng2,5):evalf(ex); rel_e
rror = evalf((ex2-ex)/ex2,10); rel_error = evalf((ex2-ex)/ex2,24); rel
_error = evalf((ex2-ex)/ex2,40);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%/
Romberg~5~(32)G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"+L$Qxf\"!\"%" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$\"+_(*3*>)!#>" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!98\\tV]a]X(*eg^!#N" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!IqOH:!>ZGU$*QwJ*R\\X(*eg^!#^
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
137 "The sensitivity to precision in the Romberg quadrature shows up c
learly here. Lets try to using an unrealistically large number of node
s:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 121 "`trapezoidal 1024`; ex:=T(exp,rng2,1024):evalf(ex); \+
rel_error = evalf((ex2-ex)/ex2,10); rel_error=evalf((ex2-ex)/ex2,24);
" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%1trapezoidal~1024G" }}{PARA 11 "
" 1 "" {XPPMATH 20 "6#$\"+/0u(f\"!\"%" }}{PARA 11 "" 1 "" {XPPMATH 20 
"6#/%*rel_errorG$!+bIoe8!#:" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_
errorG$!9h2GvcD$=Vl:F\"!#H" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
121 "`Simpson's 1024`; ex:=S(exp,rng2,1024):evalf(ex); rel_error = eva
lf((ex2-ex)/ex2,10); rel_error = evalf((ex2-ex)/ex2,24);" }}{PARA 11 "
" 1 "" {XPPMATH 20 "6#%/Simpson's~1024G" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#$\"+I$Qxf\"!\"%" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG
$\"+Q0L<M!#>" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!9Te75mU
@E7]$H\"!#N" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 193 "st:=time():
`Romberg 10 (1024)`; ex:=R(exp,rng2,10):evalf(ex); rel_error = evalf((
ex2-ex)/ex2,10); rel_error = evalf((ex2-ex)/ex2,24); rel_error = evalf
((ex2-ex)/ex2,40); `cpu time` = time()-st;" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#%2Romberg~10~(1024)G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6
#$\"+N&Qxf\"!\"%" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!+:R
\\x7!#;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$\"9\"f%QiFU$[
m.x^*!#Y" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$\"I5$o1uzeLc
6FW?%flH>BF@!#x" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%)cpu~timeG$\"%Ib!
\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 376 "When relative error changes so much with precision it is almos
t certain that we are seeing mostly roundoff error. The Romberg quadra
ture is extremely accurate, but one has to use high precision to achie
ve the accuracy. In any event Romberg is the clear winner here as far \+
truncation error is concerned but more work is required to realize the
 advantages of the Romberg method." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 400 "From a philosophical point of view the d
ifficulty with loss of significance should be expected. Extrapolation \+
depends upon cancellation and so is bound to be susceptible to loss of
 significance from roundoff unless suitable precautions are taken. For
 example, if the major part of the cancellation occurs symbolically ra
ther than in the context of floating point numbers, we would reduce th
e problem." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 382 "A small possible improvement is to re-arrange the iterative pa
rt of the Romberg calculation so that  R[j,k] is obtained by adding a \+
fairly small number to  R[j-1,k]  (at least for large j). It turns out
 this idea has no effect here, probably because Maple carries out the \+
calculation symbolically and utimately obtains the same symbolic expre
ssion as before. Here's a demonstration:" }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "Rb:=proc(f,r::range,
m)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "local k,j,R;" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 24 "if type(m, numeric) then" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 84 "    if not type(m,nonnegint) then ERROR(\"m must be a
 non-negative integer\"); fi; fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
20 "for k from 0 to m do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "  R[0,k
]:=T(f,r,2^k); od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for j from 1 \+
to m do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "  for k from j to m do" 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "  R[j,k]:=R[j-1,k]+(R[j-1,k]-R[j-
1,k-1])/(4^j-1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "od; od;" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "R[m,m];" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
4 "" 0 "" {TEXT 264 9 "Example 3" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 49 "Let's test  Rb to verify there is no impr
ovement:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 37 "rng3:=8..12: ex3:=Int(exp(x),x=rng3);" }}{PARA 11 "" 
1 "" {XPPMATH 20 "6#>%$ex3G-%$IntG6$-%$expG6#%\"xG/F+;\"\")\"#7" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 194 "st:=time():`Romberg b 10 (1
024)`; ex:=Rb(exp,rng3,10):evalf(ex); rel_error = evalf((ex3-ex)/ex3,1
0); rel_error = evalf((ex3-ex)/ex3,24); rel_error = evalf((ex3-ex)/ex3
,40); `cpu time`=time()-st;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#%4Rombe
rg~b~10~(1024)G" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"+N&Qxf\"!\"%" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!+:R\\x7!#;" }}{PARA 11 
"" 1 "" {XPPMATH 20 "6#/%*rel_errorG$\"9\"f%QiFU$[m.x^*!#Y" }}{PARA 
11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$\"I5$o1uzeLc6FW?%flH>BF@!#x" }
}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%)cpu~timeG$\"%'H'!\"$" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 81 "We can tr
y to improve efficiency a bit by noting that the trapezoidal rule with
  " }{XPPEDIT 18 0 "2*m;" "6#*&\"\"#\"\"\"%\"mGF%" }{TEXT -1 65 "  sub
intervals can be calculated from the trapezoidal rule with  " }
{XPPEDIT 18 0 "m;" "6#%\"mG" }{TEXT -1 44 "  subintervals by doing onl
y an additional  " }{XPPEDIT 18 0 "m;" "6#%\"mG" }{TEXT -1 90 "  funct
ion evaluations. Here is the code which produces  the trapezoidal quad
rature with  " }{XPPEDIT 18 0 "n;" "6#%\"nG" }{TEXT -1 34 "  subinterv
als from the one with  " }{XPPEDIT 18 0 "n/2;" "6#*&%\"nG\"\"\"\"\"#!
\"\"" }{TEXT -1 17 "  subintervals,  " }{XPPEDIT 18 0 "n;" "6#%\"nG" }
{TEXT -1 7 "  even." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 62 "TA:=proc(f,r::range,n,T)  # T should be n
/2 trapezoidal result" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "  local k,
h,a,b,m;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "if type(n, numeric) the
n" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "  if not type(n,posint) then E
RROR(\"n must be a positive integer\"); fi;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 61 "  if irem(n,2)=1 then ERROR(\"n must be an even integ
er\"); fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "fi;" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 31 "m:=n/2; a:=op(1,r); b:=op(2,r);" }}{PARA 0 "> " 0 
"" {MPLTEXT 1 0 11 "h:=(b-a)/n;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "
(1/2)*T+h*sum(f(a+(2*k-1)*h),k=1..m);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 38 "Let's introduce this code in Rb above:" }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "Rc:=proc(
f,r::range,m)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "local k,j,R;" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "if type(m, numeric) then" }}{PARA 
0 "> " 0 "" {MPLTEXT 1 0 84 "    if not type(m,nonnegint) then ERROR(
\"m must be a non-negative integer\"); fi; fi;" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 17 "R[0,0]:=T(f,r,1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 
20 "for k from 1 to m do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "  R[0,k
]:=TA(f,r,2^k,R[0,k-1]); od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for
 j from 1 to m do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "  for k from j
 to m do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "  R[j,k]:=R[j-1,k]+(R[j
-1,k]-R[j-1,k-1])/(4^j-1);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "od; od
;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "R[m,m];" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
4 "" 0 "" {TEXT 265 9 "Example 4" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "rng4:=8..12: ex4:=Int(exp(x)
,x=rng4);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$ex4G-%$IntG6$-%$expG6#
%\"xG/F+;\"\")\"#7" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 193 "st:=
time():`Romberg c 10 (1024)`; ex:=Rc(exp,rng4,10):evalf(ex); rel_error
 = evalf((ex4-ex)/ex4,10); rel_error = evalf((ex4-ex)/ex4,24); rel_err
or = evalf((ex4-ex)/ex4,40);`cpu time`=time()-st;" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#%4Romberg~c~10~(1024)G" }}{PARA 11 "" 1 "" {XPPMATH 20 
"6#$\"+M$Qxf\"!\"%" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!+
)G/B2#!#=" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!9`N1BQQZ\"
=bnL%!#Y" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%*rel_errorG$!IK.U)=r8Ey'
=<.FyOg.?h\\!#x" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%)cpu~timeG$\"&Y+
\"!\"$" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 202 "Our quest for efficiency seems headed in the wrong direc
tion! More detailed timing and a lot more information about how Maple \+
works in the background is needed. I'll leave it for you to experiment
 with!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 23 
"We could, for a given  " }{XPPEDIT 18 0 "m;" "6#%\"mG" }{TEXT -1 715 
",  write the Romberg rule as a weighted sum of the evaluations of  f \+
 at equispaced abscissas. This would give us a convenient formula. In \+
applications, the Romberg rule is often applied recursively, with the \+
change between subsequent results used as an estimate of the error. Th
is can be done quite efficiently, since previous calculations are recy
cled. Moreover, the error estimate is often good enough to be used as \+
a stopping condition (by comparing it with a pre-assigned desired prec
ision). We loose these conveniences if we pass to a formula, but on th
e other hand the weights are positive, and loss of precision is not a \+
problem (for positive integrands), or is manageable (for integrands th
at change sign)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" 
{TEXT -1 230 "The weights can be computed by a simple technique. Just \+
have Maple compute the Romberg estimates of the integrals of suitable \+
interpolation polynomials. The polynomials you use are those that are \+
1 at one node and 0 at the others." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 24 "Let's look at the case  " }{XPPEDIT 18 0 
"m = 2;" "6#/%\"mG\"\"#" }{TEXT -1 15 ". Here we have " }{XPPEDIT 18 
0 "2^m+1 = 5;" "6#/,&)\"\"#%\"mG\"\"\"\"\"\"F(\"\"&" }{TEXT -1 19 " no
des. Let's take " }{XPPEDIT 18 0 "a = 0;" "6#/%\"aG\"\"!" }{TEXT -1 5 
" and " }{XPPEDIT 18 0 "b = 1;" "6#/%\"bG\"\"\"" }{TEXT -1 1 "." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 122 "xc:=[0,1/4,1/2,3/4,1]: yc[0]:=[1,0,0,0,0]: yc[1]:=[0,1,0,0,0]: \+
yc[2]:=[0,0,1,0,0]: yc[3]:=[0,0,0,1,0]: yc[4]:=[0,0,0,0,1]:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "for k from 0 to 4 do pc[k]:=
unapply(interp(xc,yc[k],x),x); od;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#
>&%#pcG6#\"\"!R6#%\"xG6\"6$%)operatorG%&arrowGF+,,*$)9$\"\"%\"\"\"#\"#
K\"\"$*$)F2F7F4#!#!)F7*$)F2\"\"#F4#\"#qF7F2#!#DF7\"\"\"FCF+F+F+" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%#pcG6#\"\"\"R6#%\"xG6\"6$%)operator
G%&arrowGF+,**$)9$\"\"%\"\"\"#!$G\"\"\"$*$)F2F7F4\"#'**$)F2\"\"#F4#!$3
#F7F2\"#;F+F+F+" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%#pcG6#\"\"#R6#%
\"xG6\"6$%)operatorG%&arrowGF+,**$)9$\"\"%\"\"\"\"#k*$)F2\"\"$F4!$G\"*
$)F2F'F4\"#wF2!#7F+F+F+" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%#pcG6#\"
\"$R6#%\"xG6\"6$%)operatorG%&arrowGF+,**$)9$\"\"%\"\"\"#!$G\"F'*$)F2F'
F4#\"$C#F'*$)F2\"\"#F4#!$7\"F'F2#\"#;F'F+F+F+" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>&%#pcG6#\"\"%R6#%\"xG6\"6$%)operatorG%&arrowGF+,**$)9$
F'\"\"\"#\"#K\"\"$*$)F2F6F3!#;*$)F2\"\"#F3#\"#AF6F2!\"\"F+F+F+" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "for k from 0 to 4 do r[k]:=R
(pc[k],0..1,2); od;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%\"rG6#\"\"!#
\"\"(\"#!*" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%\"rG6#\"\"\"#\"#;\"#X
" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>&%\"rG6#\"\"##F'\"#:" }}{PARA 11 
"" 1 "" {XPPMATH 20 "6#>&%\"rG6#\"\"$#\"#;\"#X" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#>&%\"rG6#\"\"%#\"\"(\"#!*" }}}{EXCHG {PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 62 "These weights are well-kn
own! They correspond to Boole's rule." }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "RS4:=f->sum(r[k]*f(k/4)
,k=0..4); k:=evaln(k):" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%$RS4GR6#%
\"fG6\"6$%)operatorG%&arrowGF(-%$sumG6$*&&%\"rG6#%\"kG\"\"\"-9$6#,$F3#
F4\"\"%F4/F3;\"\"!F:F(F(F(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
28 "RS4(x->sin(Pi*x)); evalf(%);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#,&
#\"\"#\"#:\"\"\"*$-%%sqrtG6#F%\"\"\"#\"#;\"#X" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#$\"+@#[;O'!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 32 "int(sin(Pi*x),x=0..1); evalf(%);" }}{PARA 11 "" 1 "" {XPPMATH 
20 "6#,$*&\"\"\"F%%#PiG!\"\"\"\"#" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$
\"+Ax>mj!#5" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "
" {TEXT -1 225 "Note, for the interval  [0,1]  the weights will be rat
ional numbers, so in principle we can compute them exactly. In practic
e, of course, there is the risk that the numerators and denominators m
ay become inconveniently large," }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 11 "Experiment!" }}}}{MARK "40 1 0" 175 }
{VIEWOPTS 1 1 0 1 1 1803 }