# HG changeset patch # User Rik # Date 1476915633 25200 # Node ID 2d9393be9a18ef1eabb18bbd18ae8c39a56e0b97 # Parent 859d48b5648db5af297cc2eada1b4829a50e3455# Parent 9bc03a3f7a345c5b0551f6a93a19e950aad5a96c maint: Periodic merge of stable to default. diff -r 859d48b5648d -r 2d9393be9a18 scripts/ode/ode23.m --- a/scripts/ode/ode23.m Mon Oct 17 16:26:10 2016 -0400 +++ b/scripts/ode/ode23.m Wed Oct 19 15:20:33 2016 -0700 @@ -89,15 +89,6 @@ ## @seealso{odeset, odeget, ode45} ## @end deftypefn -## FIXME: We store ChangeLog information in Mercurial. Can this be deleted? -## ChangeLog: -## 20010703 the function file "ode23.m" was written by Marc Compere -## under the GPL for the use with this software. This function has been -## taken as a base for the following implementation. -## 20060810, Thomas Treichl -## This function was adapted to the new syntax that is used by the -## new OdePkg for Octave and is compatible to Matlab's ode23. - function varargout = ode23 (fun, trange, init, varargin) if (nargin < 3) @@ -161,12 +152,6 @@ ## Start preprocessing, have a look which options are set in odeopts, ## check if an invalid or unused option is set. - ## FIXME: Why persistent? Won't these have different values for every - ## run of ode23? - persistent defaults = []; - persistent classes = []; - persistent attributes = []; - [defaults, classes, attributes] = odedefaults (numel (init), trange(1), trange(end)); @@ -243,8 +228,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), @@ -312,9 +296,6 @@ %! err(i) = norm (y .* exp (t) - 1, Inf); %! endfor %! -%! ## Estimate order numerically -%! p = diff (log (err)) ./ diff (log (h)) -%! %! ## Estimate order visually %! loglog (h, tol, "-ob", %! h, err, "-b", @@ -326,6 +307,9 @@ %! title ("Convergence plot for ode23"); %! legend ("imposed tolerance", "ode23 (relative) error", %! "order 2", "order 3", "location", "northwest"); +%! +%! ## Estimate order numerically +%! p = diff (log (err)) ./ diff (log (h)) ## We are using the "Van der Pol" implementation for all tests that are done ## for this function. @@ -360,15 +344,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 %! @@ -376,8 +366,8 @@ %! [t, y] = ode23 (@fpol, [0 2], [2 0]); %! assert ([t(end), y(end,:)], [2, fref], 1e-3); %!test # anonymous function instead of real function -%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; -%! [t, y] = ode23 (fvdb, [0 2], [2 0]); +%! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%! [t, y] = ode23 (fvdp, [0 2], [2 0]); %! assert ([t(end), y(end,:)], [2, fref], 1e-3); %!test # extra input arguments passed through %! [t, y] = ode23 (@fpol, [0 2], [2 0], 12, 13, "KL"); @@ -484,9 +474,9 @@ %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); ## FIXME: Missing tests. -## test for MvPattern option is missing ## test for InitialSlope option is missing ## test for MaxOrder option is missing +## test for MvPattern option is missing ## Test input validation %!error ode23 () diff -r 859d48b5648d -r 2d9393be9a18 scripts/ode/ode45.m --- a/scripts/ode/ode45.m Mon Oct 17 16:26:10 2016 -0400 +++ b/scripts/ode/ode45.m Wed Oct 19 15:20:33 2016 -0700 @@ -143,11 +143,6 @@ ## Start preprocessing, have a look which options are set in odeopts, ## check if an invalid or unused option is set - ## FIXME: Why persistent when it is changed with every run of ode45? - persistent defaults = []; - persistent classes = []; - persistent attributes = []; - [defaults, classes, attributes] = odedefaults (numel (init), trange(1), trange(end)); @@ -224,8 +219,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), @@ -293,9 +287,6 @@ %! err(i) = norm (y .* exp (t) - 1, Inf); %! endfor %! -%! ## Estimate order numerically -%! p = diff (log (err)) ./ diff (log (h)) -%! %! ## Estimate order visually %! loglog (h, tol, "-ob", %! h, err, "-b", @@ -307,6 +298,9 @@ %! title ("Convergence plot for ode45"); %! legend ("imposed tolerance", "ode45 (relative) error", %! "order 4", "order 5", "location", "northwest"); +%! +%! ## Estimate order numerically +%! p = diff (log (err)) ./ diff (log (h)) ## We are using the Van der Pol equation for all tests that are done ## for this function. @@ -341,15 +335,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 %! @@ -360,8 +360,8 @@ %! [t, y] = ode45 (@fpol, [0 2], [2 0]); %! assert (size (t) < 20); %!test # anonymous function instead of real function -%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; -%! [t, y] = ode45 (fvdb, [0 2], [2 0]); +%! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%! [t, y] = ode45 (fvdp, [0 2], [2 0]); %! assert ([t(end), y(end,:)], [2, fref], 1e-2); %!test # string instead of function %! [t, y] = ode45 ("fpol", [0 2], [2 0]); @@ -482,9 +482,9 @@ %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); ## FIXME: Missing tests. -## test for MvPattern option is missing ## test for InitialSlope option is missing ## test for MaxOrder option is missing +## test for MvPattern option is missing %!error ode45 () %!error ode45 (1) diff -r 859d48b5648d -r 2d9393be9a18 scripts/ode/odeget.m --- a/scripts/ode/odeget.m Mon Oct 17 16:26:10 2016 -0400 +++ b/scripts/ode/odeget.m Wed Oct 19 15:20:33 2016 -0700 @@ -37,9 +37,7 @@ ## @seealso{odeset} ## @end deftypefn -## FIXME: 4th input argument "opt" is undocumented. - -function val = odeget (ode_opt, field, default = [], opt = "") +function val = odeget (ode_opt, field, default = []) validateattributes (ode_opt, {"struct"}, {"nonempty"}); validateattributes (field, {"char"}, {"nonempty"}); @@ -58,9 +56,9 @@ %!demo -%! # Return the manually changed value RelTol of the ODE options -%! # structure A. If RelTol wouldn't have been changed then an -%! # empty matrix value would have been returned. +%! ## Return the manually changed value RelTol of the ODE options +%! ## structure A. If RelTol wouldn't have been changed then an +%! ## empty matrix value would have been returned. %! %! A = odeset ("RelTol", 1e-1, "AbsTol", 1e-2); %! odeget (A, "RelTol", []) @@ -73,7 +71,9 @@ %!assert (odeget (odeset (), "Mass"), []) %!assert (odeget (odeset (), "AbsTol", 1e-9), 1e-9) %!assert (odeget (odeset ("AbsTol", 1e-9), "AbsTol", []), 1e-9) -%!assert (odeget (odeset ("foo", 42), "foo"), 42) +%!test +%! warning ("off", "Octave:invalid-input-arg", "local"); +%! assert (odeget (odeset ("foo", 42), "foo"), 42); %!error odeget () %!error odeget (1) diff -r 859d48b5648d -r 2d9393be9a18 scripts/ode/odeplot.m --- a/scripts/ode/odeplot.m Mon Oct 17 16:26:10 2016 -0400 +++ b/scripts/ode/odeplot.m Wed Oct 19 15:20:33 2016 -0700 @@ -19,80 +19,93 @@ ## Author: Thomas Treichl ## -*- texinfo -*- -## @deftypefn {} {@var{ret} =} 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 -## fvdb = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; +## fvdp = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; ## ## opt = odeset ("OutputFcn", @@odeplot, "RelTol", 1e-6); -## sol = ode45 (fvdb, [0 20], [2 0], opt); +## sol = ode45 (fvdp, [0 20], [2 0], opt); ## @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 ret = odeplot (t, y, flag, varargin) +function stop_solve = odeplot (t, y, flag) - ## No input argument check is done for a higher processing speed - persistent fig told yold counter; + ## 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 (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" - counter = 1; - fig = figure (); + if (isempty (flag)) + ## Default case, plot and return a value + idx += 1; + told(idx,1) = t(1,1); + yold(:,idx) = y(:,1); + for i = 1:num_lines + set (hlines(i), "xdata", told, "ydata", yold(i,:)); + endfor + drawnow; + + retval = false; + + elseif (strcmp (flag, "init")) + ## t is either the time slot [tstart tstop] or [t0, t1, ..., tn] + ## y is the initial value vector for the ode solution + idx = 1; told = t(1,1); yold = y(:,1); - - elseif (isempty (flag)) - ## Return something, either false for "not stopping - ## the integration" or true for "stopping the integration" - counter += 1; - figure (fig); - told(counter,1) = t(1,1); - yold(:,counter) = y(:,1); - ## FIXME: Why not use '.' rather than 'o' and skip the markersize? - ## FIXME: Why not just update the xdata, ydata properties? - ## Calling plot involves a lot of overhead. - plot (told, yold, "-o", "markersize", 1); drawnow; - ret = false; + 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")) ## Cleanup after ode solver has finished. - clear ("figure", "told", "yold", "counter"); + hlines = num_lines = told = yold = idx = []; endif @@ -100,9 +113,9 @@ %!demo -%! # Solve an anonymous implementation of the Van der Pol equation -%! # and display the results while solving -%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%! ## Solve an anonymous implementation of the Van der Pol equation +%! ## and display the results while solving +%! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; %! opt = odeset ("OutputFcn", @odeplot, "RelTol", 1e-6); -%! sol = ode45 (fvdb, [0 20], [2 0], opt); +%! sol = ode45 (fvdp, [0 20], [2 0], opt); diff -r 859d48b5648d -r 2d9393be9a18 scripts/ode/odeset.m --- a/scripts/ode/odeset.m Mon Oct 17 16:26:10 2016 -0400 +++ b/scripts/ode/odeset.m Wed Oct 19 15:20:33 2016 -0700 @@ -87,10 +87,6 @@ ## Mass matrix, specified as a constant matrix or a function of ## time and state. ## -## @item MassConstant -## Specify whether the mass matrix is a constant matrix or depends on -## the state. -## ## @item MassSingular ## Specify whether the mass matrix is singular. Accepted values include ## @qcode{"yes"}, @qcode{"no"}, @qcode{"maybe"}. @@ -165,8 +161,6 @@ p.addParameter ("JConstant", []); p.addParameter ("JPattern", []); p.addParameter ("Mass", []); - ## FIXME: MassConstant does not appear in Matlab documentation for odeset - p.addParameter ("MassConstant", []); p.addParameter ("MassSingular", []); p.addParameter ("MaxOrder", []); p.addParameter ("MaxStep", []); @@ -190,22 +184,17 @@ odestruct = p.Results; odestruct_extra = p.Unmatched; - ## FIXME: For speed, shouldn't this merge of structures only occur - ## when there is something in odestruct_extra? - ## FIXME: Should alphabetical order of fieldnames be maintained - ## by using sort? - s1 = cellfun (@(x) ifelse (iscell (x), {x}, x), - struct2cell (odestruct), - "UniformOutput", false); + xtra_fields = fieldnames (odestruct_extra); + if (! isempty (xtra_fields)) + ## Merge extra fields into existing odestruct + for fldname = sort (xtra_fields.') + fldname = fldname{1}; + warning ("Octave:invalid-input-arg", + "odeset: unknown option \"%s\"\n", fldname); + odestruct.(fldname) = odestruct_extra.(fldname); + endfor + endif - s2 = cellfun (@(x) ifelse (iscell (x), {x}, x), - struct2cell (odestruct_extra), - "UniformOutput", false); - - C = [fieldnames(odestruct) s1; - fieldnames(odestruct_extra) s2]; - - odestruct = struct (C'{:}); endif endfunction @@ -225,7 +214,6 @@ disp (' JConstant: binary, {["off"], "on"}'); disp (' JPattern: sparse matrix, []'); disp (' Mass: matrix or function_handle, []'); - disp (' MassConstant: binary, {["off"], "on"}'); disp (' MassSingular: switch, {["maybe"], "no", "yes"}'); disp (' MaxOrder: switch, {[5], 1, 2, 3, 4, }'); disp (' MaxStep: scalar, >0, []'); @@ -244,19 +232,19 @@ %!demo -%! # A new ODE options structure with default values is created. +%! ## A new ODE options structure with default values is created. %! %! odeoptA = odeset (); %!demo -%! # A new ODE options structure with manually set options -%! # for "AbsTol" and "RelTol" is created. +%! ## A new ODE options structure with manually set options +%! ## for "AbsTol" and "RelTol" is created. %! %! odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1); %!demo -%! # A new ODE options structure is created from odeoptB with -%! # a modified value for option "NormControl". +%! ## A new ODE options structure is created from odeoptB with +%! ## a modified value for option "NormControl". %! %! odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1); %! odeoptC = odeset (odeoptB, "NormControl", "on"); @@ -264,7 +252,7 @@ %!test %! odeoptA = odeset (); %! assert (isstruct (odeoptA)); -%! assert (numfields (odeoptA), 23); +%! assert (numfields (odeoptA), 22); %! assert (all (structfun ("isempty", odeoptA))); %!shared odeoptB, odeoptC @@ -284,13 +272,9 @@ ## Test custom user-defined option %!test -%! wstate = warning ("off", "Octave:invalid-input-arg"); -%! unwind_protect -%! odeopt = odeset ("NewtonTol", 3); -%! assert (odeopt.NewtonTol, 3); -%! unwind_protect_cleanup -%! warning (wstate); -%! end_unwind_protect +%! warning ("off", "Octave:invalid-input-arg", "local"); +%! odeopt = odeset ("NewtonTol", 3); +%! assert (odeopt.NewtonTol, 3); ## FIXME: Add an inexact match option once it is available in inputParser. ## See bug #49364. @@ -299,12 +283,7 @@ ## Test input validation %!error odeset ("opt1") -%!error odeset (1, 1) +%!error odeset (1, 1) %!error odeset (odeset (), "opt1") -%!error odeset (odeset (), 1, 1) +%!error odeset (odeset (), 1, 1) -##FIXME: Add not exact match option -## %!warning odeset ("Rel", 1); -## %!error odeset ("Initial", 1) - - diff -r 859d48b5648d -r 2d9393be9a18 scripts/ode/private/integrate_adaptive.m --- a/scripts/ode/private/integrate_adaptive.m Mon Oct 17 16:26:10 2016 -0400 +++ b/scripts/ode/private/integrate_adaptive.m Wed Oct 19 15:20:33 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 @@ -116,6 +116,7 @@ solution.unhandledtermination = true; ireject = 0; + NormControl = strcmp (options.NormControl, "on"); k_vals = []; iout = istep = 1; @@ -133,9 +134,8 @@ x_est(nn, end) = abs (x_est(nn, end)); endif - ## FIXME: Take strcmp out of while loop and calculate just once err = AbsRel_norm (x_new, x_old, options.AbsTol, options.RelTol, - strcmp (options.NormControl, "on"), x_est); + NormControl, x_est); ## Accept solution only if err <= 1.0 if (err <= 1) @@ -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 @@ -280,7 +289,7 @@ err += eps; # avoid divisions by zero dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1)))); dt = dir * min (abs (dt), options.MaxStep); - if (! (abs (dt) > eps (t (end)))) + if (! (abs (dt) > eps (t(end)))) break; endif diff -r 859d48b5648d -r 2d9393be9a18 scripts/ode/private/ode_event_handler.m --- a/scripts/ode/private/ode_event_handler.m Mon Oct 17 16:26:10 2016 -0400 +++ b/scripts/ode/private/ode_event_handler.m Wed Oct 19 15:20:33 2016 -0700 @@ -19,44 +19,44 @@ ## -*- texinfo -*- ## @deftypefn {} {@var{retval} =} ode_event_handler (@var{@@evtfun}, @var{t}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{}) ## -## Return the solution of the event function that is specified as the first -## input argument @var{@@evtfun} in the form of a function handle. +## Return the solution of the event function (@var{@@evtfun}) which is +## specified in the form of a function handle. +# +## The second input argument @var{t} is a scalar double and specifies the time +## of the event evaluation. ## -## The second input argument @var{t} is of type double scalar and -## specifies the time of the event evaluation, the third input argument -## @var{y} either is of type double column vector (for ODEs and DAEs) and -## specifies the solutions or is of type cell array (for IDEs and DDEs) and -## specifies the derivatives or the history values, the third input argument -## @var{flag} is of type string and can be of the form +## The third input argument @var{y} may be a column vector of type double +## (for ODEs and DAEs) which specifies the solutions. Alternatives, @var{y} +## may be a cell array (for IDEs and DDEs) which specifies the solutions and +## derivatives. +## +## The fourth input argument @var{flag} is of type string. Valid values are: ## ## @table @option ## @item @qcode{"init"} -## then initialize internal persistent variables of the function -## @code{ode_event_handler} and return an empty cell array of size 4, +## Initialize internal persistent variables of the function +## @code{ode_event_handler} and return an empty cell array of size 4. ## ## @item @qcode{"calc"} -## then do the evaluation of the event function and return the solution -## @var{retval} as type cell array of size 4, +## Evaluate the event function and return the solution @var{retval} as a cell +## array of size 4. ## ## @item @qcode{"done"} -## then cleanup internal variables of the function -## @code{ode_event_handler} and return an empty cell array of size 4. +## Clean up internal variables of the function @code{ode_event_handler} and +## return an empty cell array of size 4. ## @end table ## -## Optionally if further input arguments @var{par1}, @var{par2}, @dots{} of -## any type are given then pass these parameters through -## @code{ode_event_handler} to the event function. +## If additional input arguments @var{par1}, @var{par2}, @dots{} are given +## these parameters are passed directly to the event function. ## -## This function is an ODE internal helper function therefore it should -## never be necessary that this function is called directly by a user. There -## is only little error detection implemented in this function file to -## achieve the highest performance. +## This function is an ODE internal helper function and it should never be +## necessary to call it directly. ## @end deftypefn function retval = ode_event_handler (evtfun, t, y, flag = "", varargin) ## No error handling has been implemented in this function to achieve - ## the highest performance available. + ## the highest performance possible. ## retval{1} is true (to terminate) or false (to continue) ## retval{2} is the index information for which an event occurred @@ -70,33 +70,12 @@ ## yold the ODE result ## retcell the return values cell array ## evtcnt the counter for how often this function has been called - persistent evtold told yold retcell evtcnt; - - ## Call the event function if an event function has been defined to - ## initialize the internal variables of the event function and to get - ## a value for evtold. - if (strcmp (flag, "init")) + persistent evtold told yold retcell; + persistent evtcnt = 1; # Don't remove. Required for Octave parser. - if (! iscell (y)) - inpargs = {evtfun, t, y}; - else - inpargs = {evtfun, t, y{1}, y{2}}; - y = y{1}; # Delete cell element 2 - endif - if (nargin > 4) - inpargs = {inpargs{:}, varargin{:}}; - endif - [evtold, term, dir] = feval (inpargs{:}); - - ## FIXME: This actually seems to assume that everything must be row vectors - ## We assume that all return values must be column vectors - evtold = evtold(:)'; term = term(:)'; dir = dir(:)'; - told = t; yold = y; evtcnt = 1; retcell = cell (1,4); - - ## Process the event, i.e., - ## find the zero crossings for either a rising or falling edge - elseif (isempty (flag)) - + if (isempty (flag)) + ## Process the event, i.e., + ## find the zero crossings for either a rising or falling edge if (! iscell (y)) inpargs = {evtfun, t, y}; else @@ -108,8 +87,7 @@ endif [evt, term, dir] = feval (inpargs{:}); - ## FIXME: This actually seems to assume that everything must be row vectors - ## We assume that all return values must be column vectors + ## We require that all return values be row vectors evt = evt(:)'; term = term(:)'; dir = dir(:)'; ## Check if one or more signs of the event has changed @@ -143,12 +121,35 @@ evtcnt += 1; endif - endif # Check for one or more signs ... - evtold = evt; told = t; retval = retcell; yold = y; + endif + evtold = evt; told = t; yold = y; + retval = retcell; + + elseif (strcmp (flag, "init")) + ## Call the event function if an event function has been defined to + ## initialize the internal variables of the event function and to get + ## a value for evtold. - elseif (strcmp (flag, "done")) # Clear this event handling function - clear ("evtold", "told", "retcell", "yold", "evtcnt"); - retcell = cell (1,4); + if (! iscell (y)) + inpargs = {evtfun, t, y}; + else + inpargs = {evtfun, t, y{1}, y{2}}; + y = y{1}; # Delete cell element 2 + endif + if (nargin > 4) + inpargs = {inpargs{:}, varargin{:}}; + endif + [evtold, ~, ~] = feval (inpargs{:}); + + ## We require that all return values be row vectors + evtold = evtold(:)'; told = t; yold = y; + evtcnt = 1; + retval = retcell = cell (1,4); + + elseif (strcmp (flag, "done")) + ## Clear this event handling function + evtold = told = yold = evtcnt = []; + retval = retcell = cell (1,4); endif diff -r 859d48b5648d -r 2d9393be9a18 scripts/ode/private/odedefaults.m --- a/scripts/ode/private/odedefaults.m Mon Oct 17 16:26:10 2016 -0400 +++ b/scripts/ode/private/odedefaults.m Wed Oct 19 15:20:33 2016 -0700 @@ -33,10 +33,9 @@ "JConstant", "off", "JPattern", [], "Mass", [], - "MassConstant", "off", "MassSingular", "maybe", "MaxOrder", 5, - "MaxStep", 0.1 * abs (t0-tf), + "MaxStep", 0.1 * abs (tf - t0), "MStateDependence", "weak", "MvPattern", [], "NonNegative", [], @@ -48,7 +47,7 @@ "Stats", "off", "Vectorized", "off"); - defaults.MaxStep = (0.1 * abs (t0-tf)); + defaults.MaxStep = 0.1 * abs (tf -t0); persistent classes = struct ("AbsTol", {{"float"}}, "BDF", "char", @@ -59,7 +58,6 @@ "JConstant", "char", "JPattern", {{"float"}}, "Mass", {{"float", "function_handle"}}, - "MassConstant", "char", "MassSingular", "char", "MaxOrder", {{"float"}}, "MaxStep", {{"float"}}, @@ -83,7 +81,6 @@ "JConstant", {{"on", "off"}}, "JPattern", {{"vector"}}, "Mass", {{}}, - "MassConstant", {{"on", "off"}}, "MassSingular", {{"no", "maybe", "yes"}}, "MaxOrder", {{">=", 0, "<=", 5, "integer"}}, "MaxStep", {{"positive", "scalar", "real"}},