Mercurial > forge
comparison extra/integration/quad2dggen.m @ 0:6b33357c7561 octave-forge
Initial revision
| author | pkienzle |
|---|---|
| date | Wed, 10 Oct 2001 19:54:49 +0000 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| -1:000000000000 | 0:6b33357c7561 |
|---|---|
| 1 function int = quad2dggen(fun,xlow,xhigh,ylow,yhigh,tol) | |
| 2 % | |
| 3 %usage: int = quad2dggen('Fun','funxlow','funxhigh',ylow,yhigh) | |
| 4 %or | |
| 5 % int = quad2dggen('Fun','funxlow','funxhigh',ylow,yhigh,tol) | |
| 6 % | |
| 7 %This function is similar to QUAD or QUAD8 for 2-dimensional integration | |
| 8 %over a general 2-dimensional region, but it uses a Gaussian quadrature | |
| 9 %integration scheme. | |
| 10 %The integral is like: | |
| 11 % yhigh funxhigh(y) | |
| 12 % int = Int Int Fun(x,y) dx dy | |
| 13 % ylow funxlow(y) | |
| 14 % | |
| 15 % int -- value of the integral | |
| 16 % Fun -- Fun(x,y) (function to be integrated) | |
| 17 % funxlow -- funxlow(y) | |
| 18 % funxhigh-- funxhigh(y) | |
| 19 % ylow -- lower y limit of integration | |
| 20 % yhigh -- upper y limit of integration | |
| 21 % tol -- tolerance parameter (optional) | |
| 22 %Note that if there are discontinuities the region of integration | |
| 23 %should be broken up into separate pieces. And if there are singularities, | |
| 24 %a more appropriate integration quadrature should be used | |
| 25 %(such as the Gauss-Chebyshev for a specific type of singularity). | |
| 26 | |
| 27 %This routine could be optimized. | |
| 28 | |
| 29 if ( exist('tol') != 1) | |
| 30 tol=1e-3; | |
| 31 elseif ( tol==[] ) | |
| 32 tol=1e-3; | |
| 33 endif | |
| 34 | |
| 35 oldint=gquad2dgen(fun,xlow,xhigh,ylow,yhigh,2,2); | |
| 36 | |
| 37 converge=0; | |
| 38 for i=1:7 | |
| 39 lim = 2^(i+1); | |
| 40 int=gquad2dgen(fun,xlow,xhigh,ylow,yhigh,lim,lim); | |
| 41 | |
| 42 diff = oldint - int; | |
| 43 limit = abs(tol*int); | |
| 44 if ( abs(diff) < limit ) | |
| 45 converge=1; | |
| 46 break; | |
| 47 endif | |
| 48 oldint=int; | |
| 49 endfor | |
| 50 | |
| 51 if ( converge==0 ) | |
| 52 disp('Integral did not converge--singularity likely') | |
| 53 endif | |
| 54 | |
| 55 endfunction |