Mercurial > octave

changeset 31063:451fb63a10a0

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update integral to call quadgk for 'ArrayValued' integrations (bug #62468) * integral.m: Modify integrator selection so that calls with ‘ArrayValued’ go to quadgk instead of quadv. Remove error checks for previously incompatible parameter combinations. Update docstring to remove mention of quadv, point 'ArrayValued' reference to quadgk, change returned error parameter description to match current behavior, and remove parameter incompatibility note. Add BIST to verify combined parameter functionality, and change BISTs checking quadv err parameter. * quadgk.m: Correct parameter name in docstring. * NEWS.8.md: Under General Improvements add note about quadgk now accepting 'ArrayValued' parameter and update integral improvement description of optional returned error parameter. Under Matlab Compatibility add note about integral now accepting all parameter combinations.
author Nicholas R. Jankowski <jankowski.nicholas@gmail.com>
date Thu, 02 Jun 2022 19:56:15 -0400
parents ebba770cd852
children 74d97efb7573
files etc/NEWS.8.md scripts/general/integral.m scripts/general/quadgk.m
diffstat 3 files changed, 63 insertions(+), 83 deletions(-) [+]
line wrap: on
line diff
--- a/etc/NEWS.8.md	Thu Jun 02 17:08:32 2022 -0400
+++ b/etc/NEWS.8.md	Thu Jun 02 19:56:15 2022 -0400
@@ -9,13 +9,16 @@
 (as in previous versions).
 
 - `integral` can now output a second argument passing the error
-parameter from the underlying integrator.
+measurement used by the underlying integrator.
 
 - `perms` now accepts a second argument "unique" to return only unique
 permutations for inputs with repeated elements.  It is faster and takes
 less memory to call `perms ('aaaabbbbcccc', "unique")` than to call
 `unique (perms ('aaaabbbbcccc'), "rows")`.
 
+- `quadgk` can now accept the `ArrayValued` input parameter to handle
+array-valued input functions.
+
 ### Graphical User Interface
 
 
@@ -50,6 +53,9 @@
 - `var` and `std` now optionally output a second argument containing the mean
   or weighted mean.
 
+- `integral` can now accept the 'ArrayValued' option in combination with
+  'RelTol' and 'WayPoints'.
+
 - The default state for certain graphics properties has been made
   consistent with Matlab.
 
--- a/scripts/general/integral.m	Thu Jun 02 17:08:32 2022 -0400
+++ b/scripts/general/integral.m	Thu Jun 02 19:56:15 2022 -0400
@@ -32,11 +32,10 @@
 ## adaptive quadrature.
 ##
 ## @code{integral} is a wrapper for @code{quadcc} (general real-valued, scalar
-## integrands and limits), @code{quadgk} (integrals with specified integration
-## paths), and @code{quadv} (array-valued integrands) that is intended to
-## provide @sc{matlab} compatibility.  More control of the numerical
-## integration may be achievable by calling the various quadrature functions
-## directly.
+## integrands and limits), and @code{quadgk} (integrals with specified
+## integration paths and array-valued integrands) that is intended to provide
+## @sc{matlab} compatibility.  More control of the numerical integration may be
+## achievable by calling the various quadrature functions directly.
 ##
 ## @var{f} is a function handle, inline function, or string containing the name
 ## of the function to evaluate.  The function @var{f} must be vectorized and
@@ -63,7 +62,7 @@
 ## @var{arrayvalued} is specified as true.  This option will cause
 ## @code{integral} to perform the integration over the entire array and return
 ## @var{q} with the same dimensions as returned by @var{f}.  For more detail
-## @pxref{XREFquadv,,@code{quadv}}.
+## @pxref{XREFquadgk,,@code{quadgk}}.
 ##
 ## @item AbsTol
 ## Define the absolute error tolerance for the quadrature.  The default
@@ -74,10 +73,8 @@
 ## relative tolerance is 1e-6 (1e-4 for single).
 ## @end table
 ##
-## The optional output @var{err} contains an integration quality measure from
-## the called integrator.  This is an absolute error estimate from @code{quadcc}
-## and @code{quadgk}, and the number of function evaluations for array-valued
-## functions passed to @code{quadv}.
+## The optional output @var{err} contains the absolute error estimate used by
+## the called integrator.
 ##
 ## Adaptive quadrature is used to minimize the estimate of error until the
 ## following is satisfied:
@@ -104,12 +101,6 @@
 ## above.  If tighter tolerances are desired they must be specified.
 ## @sc{matlab} leaves the tighter tolerances appropriate for @code{double}
 ## inputs in place regardless of the class of the integration limits.
-##
-## @item
-## As a consequence of using @code{quadcc}, @code{quadgk}, and @code{quadv},
-## certain option combinations are not supported.  Currently,
-## @qcode{"ArrayValued"} cannot be combined with @qcode{"RelTol"} or
-## @qcode{"Waypoints"}.
 ## @end enumerate
 ##
 ## @seealso{integral2, integral3, quad, quadgk, quadv, quadl, quadcc, trapz,
@@ -124,7 +115,7 @@
 
   error_flag = (nargout == 2);
 
-  ## quadcc can't handle complex limits or integrands, but quadgk & quadv can.
+  ## quadcc can't handle complex limits or integrands, but quadgk can.
   ## Check for simple cases of complex limits and integrand.
   f_is_complex = false;
   if (iscomplex (a) || iscomplex (b))
@@ -145,7 +136,7 @@
         endif
       catch quaderror
         if (strcmp (quaderror.message,
-                    "quadcc: integrand F must return a single, real-valued vector"))
+                "quadcc: integrand F must return a single, real-valued vector"))
           if (error_flag)
             [q, err] = quadgk (f, a, b);
           else
@@ -171,6 +162,7 @@
     reltol = [];
     waypoints = [];
     arrayvalued = false;
+    use_quadgk = false;
 
     idx = 1;
     while (idx < nargin - 3)
@@ -186,8 +178,10 @@
           abstol = varargin{idx++};
         case "waypoints"
           waypoints = varargin{idx++}(:);
+          use_quadgk = true;
         case "arrayvalued"
           arrayvalued = varargin{idx++};
+          use_quadgk = true;
         otherwise
           error ("integral: unknown property '%s'", prop);
       endswitch
@@ -196,73 +190,51 @@
     issingle = (isa (a, "single") || isa (b, "single")
                 || isa (waypoints, "single"));
 
-    if (arrayvalued)
-      ## Pass vector-valued function to quadv, checking for conflicting params
-
-      ## FIXME: Replace warning when have array compatible call with waypoints
-      if (! isempty (waypoints))
-        warning (["integral: array-valued quadrature routine currently ", ...
-                 "unable to handle WayPoints.  WayPoints are ignored."]);
-      endif
+    if (isempty (abstol))
+      abstol = ifelse (issingle, 1e-5, 1e-10);
+    endif
+    if (isempty (reltol))
+      reltol = ifelse (issingle, 1e-4, 1e-6);
+    endif
 
-      ## FIXME: Remove warning once we have reltol compatible arrayval'd quadfn
-      if (! isempty (reltol))
-        warning (["integral: array-valued quadrature only accepts AbsTol.", ...
-                 "  RelTol ignored."]);
-      endif
-      if (isempty (abstol))
-        abstol = ifelse (issingle, 1e-5, 1e-10);
-      endif
+    if (use_quadgk)
+      ## Array valued functions or waypoint definitions require quadgk
+      ## no need to test for complex components
       if (error_flag)
-        [q, err] = quadv (f, a, b, abstol);
+        [q, err] = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol,
+                            "WayPoints", waypoints, "ArrayValued", arrayvalued);
       else
-        q = quadv (f, a, b, abstol);
+        q = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol,
+                            "WayPoints", waypoints, "ArrayValued", arrayvalued);
       endif
 
     else
-      if (isempty (abstol))
-        abstol = ifelse (issingle, 1e-5, 1e-10);
-      endif
-      if (isempty (reltol))
-        reltol = ifelse (issingle, 1e-4, 1e-6);
-      endif
-
-      if (! isempty (waypoints))
-        if (error_flag)
-          [q, err] = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol,
-                             "WayPoints", waypoints);
-        else
-          q = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol,
-                             "WayPoints", waypoints);
-        endif
-
-      else
-        if (! f_is_complex)
-          try
+      ## otherwise try quadcc first, switch to quadgk if complex test fails
+      if (! f_is_complex)
+        try
+          if (error_flag)
+            [q, err] = quadcc (f, a, b, [abstol, reltol]);
+          else
+            q = quadcc (f, a, b, [abstol, reltol]);
+          endif
+        catch quaderror
+          if (strcmp (quaderror.message,
+                "quadcc: integrand F must return a single, real-valued vector"))
             if (error_flag)
-              [q, err] = quadcc (f, a, b, [abstol, reltol]);
+              [q, err] = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol);
             else
-              q = quadcc (f, a, b, [abstol, reltol]);
+              q = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol);
             endif
-          catch quaderror
-            if (strcmp (quaderror.message,
-                        "quadcc: integrand F must return a single, real-valued vector"))
-              if (error_flag)
-                [q, err] = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol);
-              else
-                q = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol);
-              endif
-            else
-              error (quaderror.message);
-            endif
-          end_try_catch
+          else
+            error (quaderror.message);
+          endif
+        end_try_catch
+      else
+        ## Complex-valued integral
+        if (error_flag)
+          [q, err] = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol);
         else
-          ## Complex-valued integral
-          if (error_flag)
-            [q, err] = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol);
-          else
-            q = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol);
-          endif
+          q = quadgk (f, a, b, "AbsTol", abstol, "RelTol", reltol);
         endif
       endif
     endif
@@ -337,15 +309,17 @@
 %!assert (integral (@(z) log (z),1+1i,1+1i, "WayPoints", [1-1i, -1,-1i, -1+1i]),
 %!        complex (0, pi), 1e-10)
 
-## tests from quadv
 ## Test vector-valued functions
 %!assert (integral (@(x) [(sin (x)), (sin (2*x))], 0, pi, "ArrayValued", 1),
 %!        [2, 0], 1e-10)
 
 ## Test matrix-valued functions
-%!test
-%! assert (integral (@(x) [x,x,x; x,exp(x),x; x,x,x], 0, 1, "ArrayValued", 1),
-%!         [0.5,0.5,0.5; 0.5,(exp (1) - 1),0.5; 0.5,0.5,0.5], 1e-10);
+%!assert (integral (@(x) [x,x,x; x,exp(x),x; x,x,x], 0, 1, "ArrayValued", 1),
+%!                    [0.5,0.5,0.5; 0.5,(exp (1) - 1),0.5; 0.5,0.5,0.5], 1e-10);
+
+## Test combined parameters
+%!assert (integral (@(x) [sin(x), cos(x)], 0, pi, "ArrayValued", 1,
+%!                   "Waypoints", [0.5]), [2, 0], eps);
 
 ##test 2nd output
 %!test <*62412>
@@ -353,8 +327,8 @@
 %! assert (err, 0, 5*eps);  # err ~3e-16
 %! [~, err] = integral (@(x) ones (size (x)), 0, 1, "waypoints", 1);  # quadgk
 %! assert (err, 0, 1000*eps);  # err ~7e-14
-%! [~, err] = integral (@(x) ones (size (x)), 0, 1, "arrayvalued", true); # quadv
-%! assert (err, 0, 20);  # nfev ~13
+%! [~, err] = integral (@(x) ones (size (x)), 0, 1, "arrayvalued", true);  # quadgk
+%! assert (err, 0, 1000*eps);  # err ~7e-14
 
 ## Test input validation
 %!error integral (@sin)
--- a/scripts/general/quadgk.m	Thu Jun 02 17:08:32 2022 -0400
+++ b/scripts/general/quadgk.m	Thu Jun 02 19:56:15 2022 -0400
@@ -52,7 +52,7 @@
 ## operator @code{./} and all user functions to @code{quadgk} should do the
 ## same.
 ##
-## The optional argument @var{tol} defines the absolute tolerance used to stop
+## The optional argument @var{abstol} defines the absolute tolerance used to stop
 ## the integration procedure.  The default value is 1e-10 (1e-5 for single).
 ##
 ## The algorithm used by @code{quadgk} involves subdividing the integration