{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 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 "He
lvetica" 0 1 128 0 0 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 "Text Output" -1 2 1 {CSTYLE "" -1 -1 "Courier" 
1 10 0 0 255 1 0 0 0 0 0 1 3 0 3 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }
{PSTYLE "Warning" 2 7 1 {CSTYLE "" -1 -1 "" 0 1 0 0 255 1 0 0 0 0 0 0 
1 0 0 }0 0 0 -1 -1 -1 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 "" 11 12 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 }{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 16 "Gauss Quadrature" }}
{PARA 0 "" 0 "" {TEXT 256 18 "Mth 351 Aug 5 1999" }}{PARA 0 "" 0 "" 
{TEXT 258 16 "Bent E. Petersen" }}{PARA 0 "" 0 "" {TEXT 259 35 "Filena
me: gauss_quadrature_proc.mws" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 140 "This worksheet introduces some Maple pro
cedures that allow you to experiment with Gauss quadrature on any inte
rval for any number of nodes. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 154 "The nodes in Gauss quadrature are
 the roots of a Legendre polynomial. These can be found to any desired
 degree of accuracy by using the orthopoly package." }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 493 "The weights in Gauss q
uadrature are found by solving a system of linear equations, after the
 nodes are found. This system is not well conditioned however. Since G
auss quadrature with n nodes is exact for polynomials of degree up to \+
2n-1 we can find the weights by integrating suitable interpolation pol
ynomials. This method may give better results than solving the system \+
of linear equations directly, but keep in mind finding interpolation p
olynomials is not all that well conditioned either." }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 296 "We begin by writing a \+
procedure to calculate the nodes and weights for the interval [-1,1]. \+
Since Maple remembers previously evaluated expressions we will not wor
ry about multiple evaluations in our code. This simplifies it somewhat
. The result is returned as a list of lists, [Gnodes, Gweights]." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 
0 59 "gauss:=proc(n,d) # n = number of nodes, d = extra precision" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "local Gnodes,Gweights,uu,k,u,q,w;" 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "Digits:=Digits+d;" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 53 "Gnodes:=[fsolve(orthopoly[P](n,x),x)]; # P = L
egendre" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "Gweights:=[]; # initiall
y empty" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "uu:=[seq(0,k=1..n)];" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "for k from 1 to n do" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 42 "  u:=subsop(k=1,uu); # make kth ordinate 1" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "  q:=interp(Gnodes,u,x); # interpol
ation polynomial" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "  w:=int(q,x=-1
..1); # integrate to get weight" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "
  Gweights:=[op(Gweights),w]; # push onto list of weights" }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 3 "od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "Dig
its:=Digits-d;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "evalf([Gnodes,Gwe
ights]);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 
"" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 41 "Let's test this r
outine before moving on." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "gauss(5,2);" }}{PARA 12 "" 1 "" 
{XPPMATH 20 "6#7$7'$!+f%)zh!*!#5$!+,Jp%Q&F'\"\"!$\"+,Jp%Q&F'$\"+f%)zh!
*F'7'$\"+^)o#pBF'$\"+0nG'y%F'$\"+*))))))o&F'F2F0" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 185 "That looks alrigh
t! Note if you have many nodes you may want to use a larger value for \+
d in calculating the nodes and weights since calculating the weights i
s not very well conditioned." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 
0 "" 0 "" {TEXT -1 265 "We now write a routine gint() to do the actual
 Gauss quadrature to estimate the integral of a function on an interva
l [a,b]. Note gint() will generally not produce much error unless the \+
nodes and weights are way off, or the function to be integrated is pre
tty wild." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 59 "gint:=proc(f,a,b,gnw) # gnw list of Gauss nodes and
 weights" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "local ords,g,Gnodes,Gwe
ights;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "g:=t->((b-a)/2)*f(((b-a)*
t+(b+a))/2); # rescale to [-1,1]" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 
"Gnodes:=gnw[1];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "Gweights:=gnw[2
];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "ords:=map(g,Gnodes); # ordina
tes" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "evalf(linalg[dotprod](ords,G
weights,`orthogonal`));" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 124 "W
e now compare Gauss method with Simpson's rule. In order to have Simps
on's rule available we must load the student package." }}{PARA 0 "" 0 
"" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "with(stu
dent):" }}{PARA 7 "" 1 "" {TEXT -1 29 "Warning, new definition for D" 
}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 
106 "Here is a test routine to ease doing our comparisons. We round th
e errors to 6 significant decimal digits." }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "gtest:=proc(f,a,b,n)
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "local integ,serr,gerr,gwn;" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "integ:=evalf(Int(f(x),x=a..b)): # M
aple's numeric quadrature  " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "serr
:=integ-evalf(simpson(f(x),x=a..b,n)):" }}{PARA 0 "> " 0 "" {MPLTEXT 
1 0 45 "gwn:=gauss(n,2); gerr:=integ-gint(f,a,b,gwn):" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 62 "`integral` = evalf(integ,6), `Simpson error` = e
valf(serr,6), " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 " `Gauss error` = \+
evalf(gerr,6);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 31 "Now let's
 do some actual tests." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Digits:=16: gtest(sin,0,Pi,4);" }}
{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'++?!\"&/%.Simpson~erro
rG$!'vfX!\")/%,Gauss~errorG$\"':x:!#5" }}}{EXCHG {PARA 0 "> " 0 "" 
{MPLTEXT 1 0 18 "gtest(exp,-1,5,6);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6
%/%)integralG$\"'X![\"!\"$/%.Simpson~errorG$!'()Qt!\"'/%,Gauss~errorG$
\"'&y7#!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "gtest(exp,-1,
5,6);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'X![\"!\"$/%.S
impson~errorG$!'()Qt!\"'/%,Gauss~errorG$\"'&y7#!#5" }}}{EXCHG {PARA 0 
"> " 0 "" {MPLTEXT 1 0 18 "gtest(exp,-1,5,8);" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6%/%)integralG$\"'X![\"!\"$/%.Simpson~errorG$!'$yV#!\"'/%
,Gauss~errorG$\"&OS#!#8" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "
gtest(exp,-1,5,12);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"
'X![\"!\"$/%.Simpson~errorG$!'R\"*\\!\"(/%,Gauss~errorG$!\"\"!#8" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "Digits:=20: gtest(exp,-1,5,1
2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'X![\"!\"$/%.Sim
pson~errorG$!'R\"*\\!\"(/%,Gauss~errorG$\"\"*!#<" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "f:=x
->1/(1+x^8);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#>%\"fGR6#%\"xG6\"6$%)o
peratorG%&arrowGF(*&\"\"\"F-,&\"\"\"F/*$)9$\"\")F-F/!\"\"F(F(F(" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "Digits:=16: gtest(f,-1,1,2);
" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'I\\=!\"&/%.Simpson
~errorG$\"'PE=!\"'/%,Gauss~errorG$!'1j7F," }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 16 "gtest(f,-1,1,4);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6
%/%)integralG$\"'I\\=!\"&/%.Simpson~errorG$\"'\"e6#!\"(/%,Gauss~errorG
$\"'*o5\"F," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "gtest(f,-1,1
,8);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'I\\=!\"&/%.Sim
pson~errorG$!'!yQ&!\")/%,Gauss~errorG$!',g5!#5" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 17 "gtest(f,-1,1,12);" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6%/%)integralG$\"'I\\=!\"&/%.Simpson~errorG$!'1s6!\")/%,G
auss~errorG$!'Io6!#6" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "gte
st(f,-1,1,20);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'I\\=
!\"&/%.Simpson~errorG$!'B*)z!#5/%,Gauss~errorG$\"#n!#:" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 208 "Note Map
le is of the opinion that the last Gauss error has only 2 significant \+
figures. Estimates of significant figures are not all that reliable, s
o the error could be bogus. Let's check at higher precision." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "Digits:=17: gtest(f,-1,1,20)
;" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'I\\=!\"&/%.Simpso
n~errorG$!'B*)z!#5/%,Gauss~errorG$!')**[\"!#;" }}}{EXCHG {PARA 0 "> " 
0 "" {MPLTEXT 1 0 29 "Digits:=18: gtest(f,-1,1,20);" }}{PARA 11 "" 1 "
" {XPPMATH 20 "6%/%)integralG$\"'I\\=!\"&/%.Simpson~errorG$!'B*)z!#5/%
,Gauss~errorG$!'dQI!#;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "D
igits:=19: gtest(f,-1,1,20);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)int
egralG$\"'I\\=!\"&/%.Simpson~errorG$!'B*)z!#5/%,Gauss~errorG$!'`vG!#;
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "Digits:=28: gtest(f,-1,
1,20);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'I\\=!\"&/%.S
impson~errorG$!'B*)z!#5/%,Gauss~errorG$!'!e)G!#;" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 229 "The conclusion we
 make is that Gauss quadrature is remarkly accurate (in the cases we t
ested) but in order to obtain the potential accuracy we must use very \+
high precision in computing the nodes and weights if we have many node
s." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 148 "Ga
uss quadrature also works when we have integrable end-point singularit
ies. In this case the error tends to be larger, but still surprisingly
 good." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "Int(log(x),x=0..1
): %=value(%);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%$IntG6$-%#lnG6#%
\"xG/F*;\"\"!\"\"\"!\"\"" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 
"Digits:=16:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "gnw:=gauss(
6,2): gint(log,0,1,gnw); # 6 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#
$!1SMi-@\"*\\)*!#;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "gnw:=
gauss(10,2): gint(log,0,1,gnw); # 10 nodes" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#$!1=@;AqjU**!#;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 23 "Here's another example." }}{PARA 0 "
" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "Digi
ts:=16:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "f:=x->x^(-1/2): \+
Int(f(x),x=0..2): %=value(%); evalf(lhs(%));" }}{PARA 11 "" 1 "" 
{XPPMATH 20 "6#/-%$IntG6$*&\"\"\"F(*$-%%sqrtG6#%\"xGF(!\"\"/F-;\"\"!\"
\"#,$*$-F+6#F2F(F2" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"1!>YZ7F%GG!#:
" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "gnw:=gauss(10,2): gint(
f,0,2,gnw); # 10 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"1E`sEy86F
!#:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "gnw:=gauss(20,2): gi
nt(f,0,2,gnw); # 20 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"1mNAt
\"f$oF!#:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "gnw:=gauss(40,
2): gint(f,0,2,gnw); # 40 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"
1m)\\p6On8#!#:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "Digits:=8
:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "gnw:=gauss(10,2): gint
(f,0,2,gnw); # 10 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\")$Q6r#!
\"(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "gnw:=gauss(20,2): gi
nt(f,0,2,gnw); # 20 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\")3\"[y
#!\"(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "gnw:=gauss(40,2): \+
gint(f,0,2,gnw); # 40 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\")\\2
:(*\"\"(" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 304 "The last example abo
ve is a disaster. It illustrates how a large number of nodes together \+
with low precision can lead to total nonsense - we are off by 15 order
s of magnitude! The error here probably comes from the calculation of \+
the weights - it might be fun to track it down. Consider that a challe
nge!" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 318 "
Of course one has to compute the nodes and weights only once, for a gi
ven number of nodes. However, if you want to use Gauss quadrature once
 in a while, and just compute the nodes and weights on the fly, then y
ou had better restrict yourself to 10 or 20 nodes at most if you are u
sing normal precision (say 16 digits)." }}{PARA 0 "" 0 "" {TEXT -1 0 "
" }}{PARA 0 "" 0 "" {TEXT -1 185 "It should now be clear why, in spite
 of the capabilities of modern computers, tables of Gauss nodes and we
ights are still published (e.g., by the National Bureau of Standards) \+
and used." }}}}{MARK "45 4 0" 185 }{VIEWOPTS 1 1 0 1 1 1803 }