%this is an example for multiple integrals


display('Now would be a good time to look at doublegauss.m.')
display('This is my code for gaussian quadrature of a double integral taking one set of limits of integration as functions.')
display('The arguments are the integrand, outer limits, inner (functional) limits, # of outer integral nodes, # of inner integral nodes.')
display('doublegauss(fxn, a, b, c(x), d(x), m, n)')
pause
display('The example I did in class on the board:')
fnx=@(x,y) y/x+x^2
pause
cc=@(x) x
pause
display('(I do not think the use of different variables in the anonymous function definitions is necessary.)')
dd=@(z) z^2
pause
display('doublegauss(fxn,0,2,cc,dd,2,3)')
doublegauss(fnx,0,2,cc,dd,2,3)
pause
display('We know the value of the integral is 17/5=3.4.  So, m=2. n=3 is not exact, but m=n=3 is exact.')
display('doublegauss(fxn,0,2,cc,dd,3,3)')
doublegauss(fnx,0,2,cc,dd,3,3)
pause
fxn2=@(u,v) exp(v/u)
display('compute the integral in example 3')
pause
display('vpa(doublegauss(fxn2,.1,.5,@(x) x^3, @(x) x^2,5,5))')
digits(10)
vpa(doublegauss(fxn2,.1,.5,@(s) s^3, @(t) t^2,5,5))
pause
display('Example 1, Section 4.8 uses Simpson Method.  Look at the bottom of page 228.  Ouch!')
display('We do the same example using Gaussian Quadrature.')
pause
display('vpa(doublegauss(@(x,y) log(x+2*y), 1.4, 2.0, @(x) 1.0, @(x) 1.5, 3,3))')
display('Notice that my algorithm requires the y-limits to be functions of x.')
pause
vpa(doublegauss(@(a,b) log(a+2*b),1.4,2.0,@(c) 1.0, @(d) 1.5,3,3))
pause
display('Simpson needed 15 function evaluations for absolute error .0000021.')
display('Gaussian Quad needs 9 function evaluations for absolute error .0000000048.')