Quadndg has trouble with functions whose integrals evaluate to zero.
Quadndg will run though all its iterations and claim the function does
not converge if the final result is close to zero.
Steps to reproduce:
GNU Octave, version 3.0.1
octave:1> function z=f(x); z=0;endfunction
octave:2> quadndg("f",[0 0 0],[4 4 4],1e-1)
The function will take a considerable length of time to complete and
will claim that the integral does not converge.
The offending line appears to be the condition in quadndg.m, Line 39:
if ( abs(int_old-int) < abs(tol*int) )
converge=1;
break;
endif
If the value of int is zero, the rhs of the inequality is zero and the
condition will never be true. Thus the function continues to iterate
until it complete all 7 feasible iterations.
I would recommend changing the condition to
if ( norm(int_old-int)/max(norm(int),1) < abs(tol) )
This allows for integrals which approach zero and also for vector, and
even matrix valued functions. I've attached the modified file in case
anyone is interested. I haven't had much time to test it though.
Hope this helps.
--
May The Maths Be With You
function int = quadndg(fun,xlow,xhigh,tol)
%
%usage: int = quadndg('Fun',xlow,xhigh)
%or
% int = quadndg('Fun',xlow,xhigh,tol)
%
%This function is similar to QUAD or QUAD8 for n-dimensional integration,
%but it uses a Gaussian quadrature integration scheme.
% int -- value of the integral
% Fun -- Fun(x) (function to be integrated) in this case treat
% all the different values of x as different variables
% as opposed to different instances of the same variable
% x -- n length vector of coordinates
% xlow -- n length vector of lower limits of integration
% xhigh -- n length vector of upper limits 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).
%Example usage:
% function z=f(x); z=x(1)*x(2)*x(3); endfunction;
% quadndg("f",[0;0;0],[4;4;4]);
% Should evaluate to 512
%This routine could be optimized.
if ( exist('tol') != 1 )
tol=1e-3;
elseif ( tol==[] )
tol=1e-3;
endif
n=length(xlow);
nquad=2*ones(n,1);
int_old=gquadnd(fun,xlow,xhigh,nquad);
converge=0;
for i=1:7
nquad=(2^(i+1))*ones(n,1);
int=gquadnd(fun,xlow,xhigh,nquad)
if ( norm(int_old-int)/max(1,norm(int)) < abs(tol) )
converge=1;
break;
endif
int_old=int;
endfor
if ( converge==0 )
disp('Integral did not converge--singularity likely')
endif
endfunction
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
