Fwd: quadl test request for Matlab

[Dmitri A. Sergatskov]

Sorry, forgot to Cc to the list.


---------- Forwarded message ----------
From: Dmitri A. Sergatskov <address@hidden>
Date: Wed, Aug 24, 2011 at 8:50 PM
Subject: Re: quadl test request for Matlab
To: Rik <address@hidden>


On Wed, Aug 24, 2011 at 2:30 PM, Rik <address@hidden> wrote:
> On 08/24/2011 11:24 AM, Judd Storrs wrote:
>> On Wed, Aug 24, 2011 at 2:21 PM, Rik <address@hidden> wrote:
>>> In order to verify that could someone with ML access check what the
>>> following code produces?
>>>> quadl (@(x) sin (2 + 3*x).^2, 0, 10, 0.3, [])
>> ans =
>>
>>     4.8597
>>
> Well then.  It looks like it is the Octave implementation of the algorithm
> which fails for large tolerances.  I suppose we could make this an xtest
> and ignore the results for now.  I don't think anyone is going to leap to
> fix this and it does make people believe something is wrong with their
> installed version of Octave.
>

With tolerance above 0.25 the recursion never happens because in the
(is+(i1-i2) == is) condition "is" = 1/eps and condition is true.
One can force the recursion to happen at least once by doing some dirty trick
like that:


--- /usr/share/octave/3.4.0/m/general/quadl.m   2011-04-08
11:54:40.000000000 -0500
+++ quadl.m     2011-08-24 20:43:08.626997473 -0500
@@ -57,6 +57,8 @@

 function Q = quadl (f, a, b, tol, trace, varargin)
  need_warning (1);
+  global iter;
+  iter = 0;
  if (nargin < 4)
    tol = [];
  endif
@@ -142,6 +144,7 @@
 ##   Walter Gautschi, 08/03/98

 function Q = adaptlobstp (f, a, b, fa, fb, is, trace, varargin)
+  global iter;
  h = (b-a)/2;
  m = (a+b)/2;
  alpha = sqrt(2/3);
@@ -159,7 +162,7 @@
  fmrr = y(5);
  i2 = (h/6)*(fa + fb + 5*(fml+fmr));
  i1 = (h/1470)*(77*(fa+fb) + 432*(fmll+fmrr) + 625*(fml+fmr) + 672*fm);
-  if (is+(i1-i2) == is || mll <= a || b <= mrr)
+  if ((iter != 0)&&(is+(i1-i2) == is || mll <= a || b <= mrr))
    if ((m <= a || b <= m) && need_warning ())
      warning ("quadl: interval contains no more machine number");
      warning ("quadl: required tolerance may not be met");
@@ -170,6 +173,7 @@
      disp ([a, b-a, Q]);
    endif
  else
+    iter++;
    Q = (adaptlobstp (f, a, mll, fa, fmll, is, trace, varargin{:})
         + adaptlobstp (f, mll, ml, fmll, fml, is, trace, varargin{:})
         + adaptlobstp (f, ml, m, fml, fm, is, trace, varargin{:})


I do not claim it to be a proper fix...



> --Rik
>

Regards,
Dmitri.
--