Sat 22 Sep 2012 12:06:31 PM UTC, comment #2:
Hi Rik,
Sure, here's an example of integrating a function, in this case sin(x), that has random occurrences of 'NaN' in it:
function y = mysin ( x )
y = sin(x);
ind = ceil( rand(1) * size(y) );
y(ind) = NaN;
end
Integrating the non-faulty function from -pi to pi gives, as expected, zero:
octave:2> [ q , err , npoints ] = quadcc( @sin , -pi , pi )
q = -2.8189e-17
err = 7.3915e-16
npoints = 95
Integrating the faulty function in the same interval does not converge at all:
octave:5> [ q , err , npoints ] = quadcc( @mysin , -pi , pi )
q = -0.026062
err = 2.1504
npoints = 167589
However, with the code fixes (fixed and re-compiled), the integrator gets to within 13 digits of the correct result:
octave:3> [ q , err , npoints ] = quadcc( @mysin , -pi , pi )
q = 3.7813e-14
err = 6.1902e-12
npoints = 198919
Incidentally, this is a nice example of the robustness of this code, at least in the corrected version, for an integrand at which all other codes fail :)
|
Thu 20 Sep 2012 04:31:33 PM UTC, original submission:
I just got an e-mail from someone who found a bug in the c-language implementation of quadcc.
In lines 1718, 1847, 1992 and 2087 of quadcc.cc, the following
statements
iv->fx[i] = octave_NaN;
ivl->fx[i] = octave_NaN;
ivr->fx[i] = octave_NaN;
should be changed to
iv->fx[nans[i]] = octave_NaN;
ivl->fx[nans[i]] = octave_NaN;
ivr->fx[nans[i]] = octave_NaN;
respectively. I would have sent a diff-file, but I'm having trouble accessing the repository.
|