Mercurial > forge
diff extra/integration/quad2dggen.m @ 0:6b33357c7561 octave-forge
Initial revision
author | pkienzle |
---|---|
date | Wed, 10 Oct 2001 19:54:49 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extra/integration/quad2dggen.m Wed Oct 10 19:54:49 2001 +0000 @@ -0,0 +1,55 @@ +function int = quad2dggen(fun,xlow,xhigh,ylow,yhigh,tol) +% +%usage: int = quad2dggen('Fun','funxlow','funxhigh',ylow,yhigh) +%or +% int = quad2dggen('Fun','funxlow','funxhigh',ylow,yhigh,tol) +% +%This function is similar to QUAD or QUAD8 for 2-dimensional integration +%over a general 2-dimensional region, but it uses a Gaussian quadrature +%integration scheme. +%The integral is like: +% yhigh funxhigh(y) +% int = Int Int Fun(x,y) dx dy +% ylow funxlow(y) +% +% int -- value of the integral +% Fun -- Fun(x,y) (function to be integrated) +% funxlow -- funxlow(y) +% funxhigh-- funxhigh(y) +% ylow -- lower y limit of integration +% yhigh -- upper y limit of integration +% tol -- tolerance parameter (optional) +%Note that if there are discontinuities the region of integration +%should be broken up into separate pieces. And if there are singularities, +%a more appropriate integration quadrature should be used +%(such as the Gauss-Chebyshev for a specific type of singularity). + +%This routine could be optimized. + +if ( exist('tol') != 1) + tol=1e-3; +elseif ( tol==[] ) + tol=1e-3; +endif + +oldint=gquad2dgen(fun,xlow,xhigh,ylow,yhigh,2,2); + +converge=0; +for i=1:7 + lim = 2^(i+1); + int=gquad2dgen(fun,xlow,xhigh,ylow,yhigh,lim,lim); + + diff = oldint - int; + limit = abs(tol*int); + if ( abs(diff) < limit ) + converge=1; + break; + endif + oldint=int; +endfor + +if ( converge==0 ) + disp('Integral did not converge--singularity likely') +endif + +endfunction