{VERSION 4 0 "IBM INTEL NT" "4.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Output" 2 20 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "Helvetica" 1 14 128 0 0 1 0 0 2 0 0 0 0 0 0 1 } {CSTYLE "" -1 257 "" 0 24 0 0 255 1 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "Helvetica" 0 1 128 0 0 1 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "Helvetica" 0 1 128 0 0 1 0 0 0 0 0 0 0 0 0 1 }{PSTYLE "Normal" -1 0 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 }{PSTYLE "Maple Output" -1 11 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 3 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Maple Output" -1 12 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 3 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Title" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 18 0 0 0 1 2 1 1 2 2 2 1 1 1 1 }1 1 0 0 12 12 1 0 1 0 2 2 19 1 }} {SECT 0 {EXCHG {PARA 256 "" 0 "" {TEXT 257 16 "Gauss Quadrature" }} {PARA 0 "" 0 "" {TEXT 256 27 "Mth 351 Aug 5 2002 Maple 6" }}{PARA 0 " " 0 "" {TEXT 258 16 "Bent E. Petersen" }}{PARA 0 "" 0 "" {TEXT 259 39 "Filename: 351u2002_gauss_quadrature.mws" }}{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 140 "This wor ksheet introduces some Maple procedures that allow you to experiment w ith Gauss quadrature on any interval for any number of nodes. " }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 154 "The node s in Gauss quadrature are the roots of a Legendre polynomial. These ca n be found to any desired degree of accuracy by using the orthopoly pa ckage." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 316 "The weights in Gauss quadrature are found by solving a system of \+ linear equations, after the nodes are found. Since Gauss quadrature wi th n nodes is exact for polynomials of degree up to 2n-1 we can also f ind the weights by integrating suitable interpolation polynomials. The second procedure is the one we will use." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 296 "We begin by writing a procedure t o calculate the nodes and weights for the interval [-1,1]. Since Maple remembers previously evaluated expressions we will not worry about mu ltiple evaluations in our code. This simplifies it somewhat. The resul t is returned as a list of lists, [Gnodes, Gweights]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 258 "Note we use an option al parameter to specify the precision for intermediate calculations. T he answer however is rounded to \"Digits\" decimal digits. Note also i f Digits is changed inside a procedure it reverts to its previous valu e when the procedure returns." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "gauss:=proc(n) # n is numbe r of nodes" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "local Gnodes,Gweights ,uu,k,u,q,w;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "if nargs > 1 then" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 " Digits:=args[2];" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 3 "fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "Gnode s:=[fsolve(orthopoly[P](n,x),x)]; # P = Legendre" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "Gweights:=[];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "u u:=[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 53 " u:=subsop(k=1.0,uu); \+ # make kth ordinate 1.0" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 56 " \+ q:=interp(Gnodes,u,x); # interpolation polynomial" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 55 " w:=int(q,x=-1..1); # integrate to g et weight" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 57 " Gweights:=[op(Gweigh ts),w]; # push onto list of weights" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "evalf([Gnodes,Gweights]); " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "[Gnodes,Gweights];" }}{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 routine before moving on." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "gauss(5);" }}{PARA 12 "" 1 "" {XPPMATH 20 "6#7$7'$!+f% )zh!*!#5$!+,Jp%Q&F'$\"\"!F+$\"+,Jp%Q&F'$\"+f%)zh!*F'7'$\"+e)o#pBF'$\"+ (o'G'y%F'$\"+\"*)))))o&F'$\"+*p'G'y%F'$\"+^)o#pBF'" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 55 "That looks alrigh t! Let's try specifying the precision." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "gauss(5,4);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#7$7'$!%i!*!\"%$!%&Q&F'$\"\"!F+$\"%&Q&F'$\"%i !*F'7'$\"%'Q#F'$\"%vZF'$\"%*o&F'$\"%$y%F'$\"%rBF'" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 185 "We now write a ro utine gint() to do the actual Gauss quadrature to estimate the integra l of a function on an interval [a,b]. Again we use an optional paramet er to specify the precision." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "gint:=proc(f,R,gnw) # gnw l ist of Gauss nodes and weights" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 41 " \+ local a,b,ords,g,Gnodes,Gweights,k,ans;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 " a:=lhs(R); b:=rhs(R);" }}{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 17 "if nargs > 3 then" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 " Digits:=args[4];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "fi;" }}{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 20 "ords:=map(g,Gnodes);" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 7 "ans:=0;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "for k \+ from 1 to nops(Gweights) do" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 " an s:=ans + ords[k]*Gweights[k];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 3 "od; " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "evalf(ans);" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 20 "Here's a quick test." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "Int(sin(x), x=0..Pi):% = evalf(%); Gauss = gint(sin,0..Pi,gauss(4));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%$IntG6$-%$sinG6#%\"xG/F*;\"\"!%#PiG$\"\"#F- " }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/%&GaussG$\"+GU)***>!\"*" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 42 "Th at's pretty impressive for just 4 nodes." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 154 "We now compare Gauss meth od with Simpson's rule. In order to have Simpson's rule available we m ust load the student package, or at least simpson() command" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "w ith(student,simpson):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 268 "Here is a test routine to ease doing our comparisons. Note n is the number of nodes in the Gauss quadrature \+ and n is the number of subintervals in Simpson's rule so we really hav e n+1 nodes and n must be even. We round the errors to 6 significa nt decimal digits." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "gtest:=proc(f,R,n)" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "local integ,serr,gerr,gwn;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "integ:=evalf(Int(f(x),x=R)): # Maple's numeric quadra ture " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "serr:=integ-evalf(simpson (f(x),x=R,n)):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "gwn:=gauss(n,Digi ts+2); gerr:=integ-gint(f,R,gwn):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 " `integral` = evalf(integ,6), `Simpson error` = evalf(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 31 "Digits:=16: gtest(sin,0..Pi,2);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'++?!\"&/%.Simpson~errorG$!'^R%*!\"(/%,G auss~errorG$\"'/=kF," }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "gte st(sin,0..Pi,4);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'++ ?!\"&/%.Simpson~errorG$!'vfX!\")/%,Gauss~errorG$\"':x:!#5" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "gtest(sin,0..Pi,6);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'++?!\"&/%.Simpson~errorG$!'!>j)!\" */%,Gauss~errorG$\"'HF_!#:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "Digits:=24: gtest(sin,0..Pi,8);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6 %/%)integralG$\"'++?!\"&/%.Simpson~errorG$!'q\"p#!\"*/%,Gauss~errorG$ \"'dRY!#?" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "gtest(sin,0..P i,10);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'++?!\"&/%.Si mpson~errorG$!'<&4\"!\"*/%,Gauss~errorG$\"%H:!#B" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "Digits:=30: gtest(sin,0..Pi,12);" }}{PARA 11 "" 1 "" {XPPMATH 20 "6%/%)integralG$\"'++?!\"&/%.Simpson~errorG$!'Vi_! #5/%,Gauss~errorG$\"%1B!#H" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 291 "We see that Gauss quadrature is remarkly accurate (in the cases we tested above). However in order to obtain t he potential accuracy we must use very high precision in computing the nodes and weights if we have many nodes. Otherwise we may see a great deal of loss of significance (roundoff)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 137 "Error estimates for Gauss quadrat ure suggest that smooth (many times continuously differentiable) integ rands will yield the best results." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 148 "Gauss quadrature also works when we have integrable end-point singularities. In this case the error tends to b e larger, but still surprisingly good." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "Int(log(x),x=0..1): %=v alue(%);" }}{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): g int(log,0,1,gnw); # 6 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$!1SMi- @\"*\\)*!#;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "gnw:=gauss(1 0,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 103 "Note we can not use Simpson's rule (directly) her e since the logarithm has a singularity at the origin." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 23 "Here's another examp le." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "Digits:=16:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "f:=x->x^(-1/2): Int(f(x),x=0..2): %=value(%); rhs(%)= evalf(lhs( %));" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/-%$IntG6$*&\"\"\"F(*$-%%sqrtG 6#%\"xGF(!\"\"/F-;\"\"!\"\"#,$*$-F+6#F2F(F2" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#/,$*$-%%sqrtG6#\"\"#\"\"\"F)$\"?U[uP.w4!>YZ7F%GG!#H" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "gnw:=gauss(10): gint(f,0..2 ,gnw); # 10 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"?96it9qpI`sEy8 6F!#H" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "gnw:=gauss(20): gi nt(f,0..2,gnw); # 20 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"?N\"[ DMM,T*>Wt\"f$oF!#H" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "gnw:= gauss(40): gint(f,0..2,gnw); # 40 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\"?\\0#4RJ=^$3***HB!)z#!#H" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "Digits:=8:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "gnw:=gauss(10): gint(f,0..2,gnw); # 10 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$\")?:6F!\"(" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "gnw:=gauss(20): gint(f,0..2,gnw); # 20 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$!)@#)Gf!\")" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "gnw:=gauss(40): gint(f,0..2,gnw); # 40 nodes" }}{PARA 11 "" 1 "" {XPPMATH 20 "6#$!)$y1\"=\"#6" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 332 " The last two examples above are a disaster. They illustrates how a lar ge number of nodes together with low precision can lead to total nonse nse - we are off by 19 orders of magnitude in the last example! The er ror here probably comes from the calculation of the weights - it might be fun to track it down. Consider that a challenge!" }}{PARA 0 "" 0 " " {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 319 "Of course one has to c ompute the nodes and weights only once, for a given number of nodes. H owever, if you want to use Gauss quadrature once in a while, and just \+ compute the nodes and weights on the fly, then you had better restrict yourself to 10 or 20 nodes at most if you are using normal precision \+ (say 16 digits). " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 200 "It should now be clear why, in spite of the capabilities of modern computers, high precision tables of Gauss nodes and weights are still published (e.g., by the National Bureau of Standards) and u sed." }}}}{MARK "39 4 0" 93 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 0 1 2 33 1 1 }