# HG changeset patch
# User David Bateman <dbateman@free.fr>
# Date 1251753847 -7200
# Node ID b18a50c56144884e069a08d35fcc478976d554d4
# Parent  8e42bb4ad34d0270b9c330329fd36dea589512e3
More care with contour integrals in quadgk

diff -r 8e42bb4ad34d -r b18a50c56144 scripts/ChangeLog
--- a/scripts/ChangeLog	Sun Aug 30 21:56:18 2009 +0200
+++ b/scripts/ChangeLog	Mon Aug 31 23:24:07 2009 +0200
@@ -1,3 +1,8 @@
+2009-08-31  David Bateman  <dbateman@free.fr>
+
+	* general/quadgk.m: More care with the interval length and
+	convergence of contour integrals.
+
 2009-08-29  John W. Eaton  <jwe@octave.org>
 
 	* time/datestr.m: Add missing semicolon.
diff -r 8e42bb4ad34d -r b18a50c56144 scripts/general/quadgk.m
--- a/scripts/general/quadgk.m	Sun Aug 30 21:56:18 2009 +0200
+++ b/scripts/general/quadgk.m	Mon Aug 31 23:24:07 2009 +0200
@@ -184,7 +184,8 @@
     if (iscomplex (a) || iscomplex (b) || iscomplex(waypoints))
       ## contour integral, no transform
       subs = [a; waypoints; b];
-      h = b - a;
+      h = sum (abs (diff (subs)));
+      h0 = h;
       trans = @(t) t;
     elseif (isinf (a) && isinf(b))
       ## Standard Infinite to finite integral transformation.
@@ -201,6 +202,7 @@
 	subs = linspace (-1, 1, 11)'; 
       endif
       h = 2;
+      h0 = b - a;
       trans = @(t) t ./ (1 - t.^2);
       f = @(t) f (t ./ (1 - t .^ 2)) .* (1 + t .^ 2) ./ ((1 - t .^ 2) .^ 2);
     elseif (isinf(a))
@@ -226,6 +228,7 @@
 	subs = linspace (0, 1, 11)'; 
       endif
       h = 1;
+      h0 = b - a;
       trans = @(t) b - (t ./ (1 + t)).^2;
       f = @(s) - 2 * s .* f (b -  (s ./ (1 + s)) .^ 2) ./ ((1 + s) .^ 3);
     elseif (isinf(b))
@@ -250,6 +253,7 @@
 	subs = linspace (0, 1, 11)'; 
       endif
       h = 1;
+      h0 = b - a;
       trans = @(t) a + (t ./ (1 - t)).^2;
       f = @(s) 2 * s .* f (a +  (s ./ (1 - s)) .^ 2) ./ ((1 - s) .^ 3);
     else
@@ -271,6 +275,7 @@
 	subs = linspace(-1, 1, 11)'; 
       endif
       h = 2;
+      h0 = b - a;
       trans = @(t) ((b - a) ./ 4) * t .* (3 - t.^2) + (b + a) ./ 2;
       f = @(t) f((b - a) ./ 4 .* t .* (3 - t.^2) + (b + a) ./ 2) .* ...
            3 .* (b - a) ./ 4 .* (1 - t.^2);
@@ -305,7 +310,7 @@
 	## Check for sub-intervals that are too small. Test must be
 	## performed in untransformed sub-intervals. What is a good
 	## value for this test. Shampine suggests 100*eps
-	if (any (diff (trans (subs), [], 2) / (b - a) < 100 * myeps))
+	if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps))
 	  q = q0;
 	  err = err0;
 	  break;
@@ -329,7 +334,7 @@
 	endif
 
 	## Accept the sub-intervals that meet the convergence criteria
-	idx = find (abs (q_errs) < tol .* diff (subs, [], 2) ./ h);
+	idx = find (abs (q_errs) < tol .* abs(diff (subs, [], 2)) ./ h);
 	if (first)
 	  q = sum (q_subs (idx));
 	  err = sum (q_errs(idx));