Mercurial > forge

comparison extra/integration/quad2dg.m @ 0:6b33357c7561 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 = quad2dg(fun,xlow,xhigh,ylow,yhigh,tol)
2 %
3 %usage: int = quad2dg('Fun',xlow,xhigh,ylow,yhigh)
4 %or
5 % int = quad2dg('Fun',xlow,xhigh,ylow,yhigh,tol)
6 %
7 %This function is similar to QUAD or QUAD8 for 2-dimensional integration,
8 %but it uses a Gaussian quadrature integration scheme.
9 % int -- value of the integral
10 % Fun -- Fun(x,y) (function to be integrated)
11 % xlow -- lower x limit of integration
12 % xhigh -- upper x limit of integration
13 % ylow -- lower y limit of integration
14 % yhigh -- upper y limit of integration
15 % tol -- tolerance parameter (optional)
16 %Note that if there are discontinuities the region of integration
17 %should be broken up into separate pieces. And if there are singularities,
18 %a more appropriate integration quadrature should be used
19 %(such as the Gauss-Chebyshev for a specific type of singularity).
20
21 %This routine could be optimized.
22
23 if ( exist('tol') != 1 )
24 tol=1e-3;
25 elseif ( tol==[])
26 tol=1e-3;
27 endif
28
29 n=length(xlow);
30 nquad=2*ones(n,1);
31 [bpx,bpy,wfxy] = grule2d(2,2);
32 int_old=gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy);
33
34 converge=0;
35 for i=1:7
36 lim = 2^(i+1);
37 [bpx,bpy,wfxy] = grule2d(lim,lim);
38 int=gquad2d(fun,xlow,xhigh,ylow,yhigh,bpx,bpy,wfxy);
39
40 if ( abs(int_old-int) < abs(tol*int) )
41 converge=1;
42 break;
43 endif
44 int_old=int;
45 endfor
46
47 if ( converge==0 )
48 disp('Integral did not converge--singularity likely')
49 endif
50
51 endfunction