Mercurial > octave

changeset 22639:7efa2d0e22c9 stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
More Matlab-compatible implementation of OutputFcn. Clean up odeplot.m. * ode23.m, ode45.m: Call OutputFcn with null inputs when using "done" flag. Rewrite BIST OutputFcn fout to follow expected conventions. * integrate_adaptive.m: Rename variable pltret to stop_solve. Fix issue where only final result of OutputFcn could stop integration. * odeplot.m: Rewrite docstring. Rename output variable to stop_solve. Always return false from function since odeplot never interrupts integration. Set x limits of plot in "init" for better performance.
author Rik <rik@octave.org>
date Tue, 18 Oct 2016 11:31:13 -0700
parents 0128795eeac6
children 58c1c6aeb737
files scripts/ode/ode23.m scripts/ode/ode45.m scripts/ode/odeplot.m scripts/ode/private/integrate_adaptive.m
diffstat 4 files changed, 81 insertions(+), 54 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/ode/ode23.m	Mon Oct 17 14:03:45 2016 -0700
+++ b/scripts/ode/ode23.m	Tue Oct 18 11:31:13 2016 -0700
@@ -243,8 +243,7 @@
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
-    feval (odeopts.OutputFcn, solution.t(end),
-           solution.x(end,:)', "done", odeopts.funarguments{:});
+    feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
     ode_event_handler (odeopts.Events, solution.t(end),
@@ -360,15 +359,21 @@
 %!  mas = sparse ([1, 0; 0, 1]);    # A sparse dummy matrix
 %!endfunction
 %!function out = fout (t, y, flag, varargin)
-%!  if (regexp (char (flag), "init") == 1)
-%!    if (any (size (t) != [2, 1])) error ('"fout" step "init"'); endif
+%!  out = false;
+%!  if (strcmp (flag, "init"))
+%!    if (! isequal (size (t), [2, 1]))
+%!      error ('fout: step "init"');
+%!    endif
 %!  elseif (isempty (flag))
-%!    if (any (size (t) != [1, 1])) error ('"fout" step "calc"'); endif
-%!    out = false;
-%!  elseif (regexp (char (flag), "done") == 1)
-%!    if (any (size (t) != [1, 1])) error ('"fout" step "done"'); endif
+%!    if (! isequal (size (t), [1, 1]))
+%!      error ('fout: step "calc"');
+%!    endif
+%!  elseif (strcmp (flag, "done"))
+%!    if (! isempty (t))
+%!      warning ('fout: step "done"');
+%!    endif
 %!  else
-%!    error ('"fout" invalid flag');
+%!    error ("fout: invalid flag <%s>", flag);
 %!  endif
 %!endfunction
 %!
--- a/scripts/ode/ode45.m	Mon Oct 17 14:03:45 2016 -0700
+++ b/scripts/ode/ode45.m	Tue Oct 18 11:31:13 2016 -0700
@@ -224,8 +224,7 @@
 
   ## Postprocessing, do whatever when terminating integration algorithm
   if (odeopts.haveoutputfunction)  # Cleanup plotter
-    feval (odeopts.OutputFcn, solution.t(end),
-           solution.x(end,:)', "done", odeopts.funarguments{:});
+    feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:});
   endif
   if (! isempty (odeopts.Events))   # Cleanup event function handling
     ode_event_handler (odeopts.Events, solution.t(end),
@@ -341,15 +340,21 @@
 %!  mas = sparse ([1, 0; 0, 1]);   # A sparse dummy matrix
 %!endfunction
 %!function out = fout (t, y, flag, varargin)
-%!  if (regexp (char (flag), "init") == 1)
-%!    if (any (size (t) != [2, 1])) error ('"fout" step "init"'); endif
+%!  out = false;
+%!  if (strcmp (flag, "init"))
+%!    if (! isequal (size (t), [2, 1]))
+%!      error ('fout: step "init"');
+%!    endif
 %!  elseif (isempty (flag))
-%!    if (any (size (t) != [1, 1])) error ('"fout" step "calc"'); endif
-%!    out = false;
-%!  elseif (regexp (char (flag), 'done') == 1)
-%!    if (any (size (t) != [1, 1])) error ('"fout" step "done"'); endif
+%!    if (! isequal (size (t), [1, 1]))
+%!      error ('fout: step "calc"');
+%!    endif
+%!  elseif (strcmp (flag, "done"))
+%!    if (! isempty (t))
+%!      warning ('fout: step "done"');
+%!    endif
 %!  else
-%!    error ('"fout" invalid flag');
+%!    error ("fout: invalid flag <%s>", flag);
 %!  endif
 %!endfunction
 %!
--- a/scripts/ode/odeplot.m	Mon Oct 17 14:03:45 2016 -0700
+++ b/scripts/ode/odeplot.m	Tue Oct 18 11:31:13 2016 -0700
@@ -19,38 +19,34 @@
 ## Author: Thomas Treichl <treichl@users.sourceforge.net>
 
 ## -*- texinfo -*-
-## @deftypefn {} {@var{retval} =} odeplot (@var{t}, @var{y}, @var{flag})
+## @deftypefn {} {@var{stop_solve} =} odeplot (@var{t}, @var{y}, @var{flag})
 ##
-## Open a new figure window and plot input @var{y} over time during the
-## solving of an ode problem.
+## Open a new figure window and plot the solution of an ode problem at each
+## time step during the integration.
 ##
-## The input @var{y} is a column vector.  The types and values of the input
-## parameter @var{t} and output parameter @var{ret} depend on the input
-## @var{flag} that is of type string.  If @var{flag} is
+## The types and values of the input parameters @var{t} and @var{y} depend on
+## the input @var{flag} that is of type string.  Valid values of @var{flag}
+## are:
 ##
 ## @table @option
 ## @item @qcode{"init"}
-## then @var{t} must be a column vector of length 2 with the first and
-## the last time step;  Nothing is returned from this function.
+## The input @var{t} must be a column vector of length 2 with the first and
+## last time step (@code{[@var{tfirst} @var{tlast}]}.  The input @var{y}
+## contains the initial conditions for the ode problem (@var{y0}).
 ##
 ## @item @qcode{""}
-## then @var{t} must be a scalar double specifying the actual time step;
-## The return value is false (resp. value 0) for @qcode{"not stop solving"}.
+## The input @var{t} must be a scalar double specifying the time for which
+## the solution in input @var{y} was calculated.
 ##
 ## @item @qcode{"done"}
-## then @var{t} must be a scalar double specifying the last time step;
-## Nothing is returned from this function.
+## The inputs should be empty, but are ignored if they are present.
 ## @end table
 ##
-## This function is called by an ode solver function if it was specified in
-## an options structure with @code{odeset}.  This function is an internal
-## helper function; It should never be necessary for this function to be
-## directly called by a user.  There is only minimal error detection
-## implemented in order to to achieve the highest performance.
+## @code{odeplot} always returns false, i.e., don't stop the ode solver.
 ##
-## For example, solve an anonymous implementation of the
+## Example: solve an anonymous implementation of the
 ## @nospell{@qcode{"Van der Pol"}} equation and display the results while
-## solving
+## solving.
 ##
 ## @example
 ## @group
@@ -61,19 +57,31 @@
 ## @end group
 ## @end example
 ##
-## @seealso{odeset, odeget}
+## Background Information:
+## This function is called by an ode solver function if it was specified in
+## the @qcode{"OutputFcn"} property of an options structure created with
+## @code{odeset}.  The ode solver will initially call the function with the
+## syntax @code{odeplot ([@var{tfirst}, @var{tlast}], @var{y0}, "init")}.  The
+## function initializes internal variables, creates a new figure window, and
+## sets the x limits of the plot.  Subsequently, at each time step during the
+## integration the ode solver calls @code{odeplot (@var{t}, @var{y}, [])}.
+## At the end of the solution the ode solver calls
+## @code{odeplot ([], [], "end")} so that odeplot can perform any clean-up
+## actions required.
+## @seealso{odeset, odeget, ode23, ode45}
 ## @end deftypefn
 
-function retval = odeplot (t, y, flag)
+function stop_solve = odeplot (t, y, flag)
 
   ## No input argument checking is done for better performance
   persistent hlines num_lines told yold;
   persistent idx = 1;   # Don't remove.  Required for Octave parser.
 
+  ## odeplot never stops the integration
+  stop_solve = false;
+
   if (isempty (flag))
     ## Default case, plot and return a value
-    ## FALSE for "not stopping the integration"
-    ## TRUE  for "stopping the integration"
     idx += 1;
     told(idx,1) = t(1,1);
     yold(:,idx) = y(:,1);
@@ -85,14 +93,14 @@
     retval = false;
 
   elseif (strcmp (flag, "init"))
-    ## Nothing to return
     ## t is either the time slot [tstart tstop] or [t0, t1, ..., tn]
-    ## y is the initial value vector "init"
+    ## y is the initial value vector for the ode solution
     idx = 1;
     told = t(1,1);
     yold = y(:,1);
     figure ();
     hlines = plot (told, yold, "-", "marker", ".", "markersize", 9);
+    xlim ([t(1), t(end)]);  # Fix limits which also speeds up plotting
     num_lines = numel (hlines);
 
   elseif (strcmp (flag, "done"))
--- a/scripts/ode/private/integrate_adaptive.m	Mon Oct 17 14:03:45 2016 -0700
+++ b/scripts/ode/private/integrate_adaptive.m	Tue Oct 18 11:31:13 2016 -0700
@@ -96,8 +96,8 @@
     else
       solution.retout = x;
     endif
-    feval (options.OutputFcn, tspan, solution.retout,
-           "init", options.funarguments{:});
+    feval (options.OutputFcn, tspan, solution.retout, "init",
+           options.funarguments{:});
   endif
 
   ## Initialize the EventFcn
@@ -186,7 +186,7 @@
           endif
 
           ## Call OutputFcn only if a valid result has been found.
-          ## Stop integration if function returns false.
+          ## Stop integration if function returns true.
           if (options.haveoutputfunction)
             cnt = options.Refine + 1;
             approxtime = linspace (t_old, t_new, cnt);
@@ -196,12 +196,16 @@
             if (! isempty (options.OutputSel))
               approxvals = approxvals(options.OutputSel, :);
             endif
+            stop_solve = false;
             for ii = 1:numel (approxtime)
-              pltret = feval (options.OutputFcn, approxtime(ii),
-                              approxvals(:, ii), [],
-                              options.funarguments{:});
+              stop_solve = feval (options.OutputFcn,
+                                  approxtime(ii), approxvals(:, ii), [],
+                                  options.funarguments{:});
+              if (stop_solve)
+                break;  # break from inner loop
+              endif
             endfor
-            if (pltret)  # Leave main loop
+            if (stop_solve)  # Leave main loop
               solution.unhandledtermination = false;
               break;
             endif
@@ -230,7 +234,7 @@
         endif
 
         ## Call OutputFcn only if a valid result has been found.
-        ## Stop integration if function returns false.
+        ## Stop integration if function returns true.
         if (options.haveoutputfunction)
           cnt = options.Refine + 1;
           approxtime = linspace (t_old, t_new, cnt);
@@ -240,11 +244,16 @@
           if (! isempty (options.OutputSel))
             approxvals = approxvals(options.OutputSel, :);
           endif
+          stop_solve = false;
           for ii = 1:numel (approxtime)
-            pltret = feval (options.OutputFcn, approxtime(ii),
-                            approxvals(:, ii), [], options.funarguments{:});
+            stop_solve = feval (options.OutputFcn,
+                                approxtime(ii), approxvals(:, ii), [],
+                                options.funarguments{:});
+            if (stop_solve)
+              break;  # break from inner loop
+            endif
           endfor
-          if (pltret)  # Leave main loop
+          if (stop_solve)  # Leave main loop
             solution.unhandledtermination = false;
             break;
           endif