Mercurial > octave
changeset 31263:449ed6f427cb
ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
* scripts/ode/ode23.m: Remove disabling of Refine option with struct output.
Modify solution struct to output two sets of solution variables: output_t,
output_x and ode_t and ode_x, and transpose struct output variables for
improved Matlab compatibility. Update BISTs and perform minor code formatting.
* scripts/ode/ode23s.m: Make same changes as ode23.m.
* scripts/ode/ode45.m: Make same changes as ode23.m. Remove comment indicating
that Refine is not implemented.
* scripts/ode/private/integrate_adaptive.m: Update internal handling of
variables t and x, separating them into ode_t & ode_x for internal integration
and output_t & output_x for function output or calls to OutputFcn. Replace
prior attempt at Refine option with new implementation. Specify time output or
Refine != 0 are both interpolated from internal variables (ode_t, ode_x) for
output of non-struct variables and/or for use with OutputFcn. Improve event
handling when multiple Events (including at least one terminal Event) are
detected in a single simulation step so that all Events up to and including the
first terminal one are reported, and final data point is set to that of
terminal Event. Send multiple data points in a single call to OutputFcn if they
are all interpolated from a single integration step. Remove warning for
termination when term signal is received from Events or OutputFcn. Return both
internal variables (ode_t, ode_x) and interpolated variables (output_t,
output_x) to allow calling function to correctly return either struct or
separate variables.
* scripts/ode/private/ode_event_handler.m: Sort multiple Events in ascending
time order when they are encountered in one integration step. Remove any events
after the time of a terminal Event.
* scripts/ode/odeset.m: Update docstring to remove indication that Refine is
not implemented
* scripts/ode/odeplot.m: Update docstring to indicate that input t can be a
scalar or vector. Add file test.
* etc/NEWS.8.md: Add descriptions of changes under General improvements and
Matlab compatibility.
| author | Ken Marek <marek_ka@mercer.edu> |
|---|---|
| date | Wed, 05 Oct 2022 16:53:01 -0400 |
| parents | b8f4ec18e728 |
| children | af3752d9a59f |
| files | etc/NEWS.8.md scripts/ode/ode23.m scripts/ode/ode23s.m scripts/ode/ode45.m scripts/ode/odeplot.m scripts/ode/odeset.m scripts/ode/private/integrate_adaptive.m scripts/ode/private/ode_event_handler.m |
| diffstat | 8 files changed, 261 insertions(+), 225 deletions(-) [+] |
line wrap: on
line diff
--- a/etc/NEWS.8.md Wed Oct 05 12:03:52 2022 -0700 +++ b/etc/NEWS.8.md Wed Oct 05 16:53:01 2022 -0400 @@ -29,6 +29,9 @@ outward normal vectors. Input type checking has also been added for improved error handling. +- `Refine` option is now implemented in functions `ode45`, `ode23`, + and `ode23s`. + ### Graphical User Interface @@ -75,6 +78,10 @@ Object | Property | Default State ------------|------------------|------------ `figure` | `"dockcontrols"` | `"on"` + +- `ode45`, `ode23`, and `ode23s` have improved results for options `Events`, + `OutputFcn`, and `Refine`, along with corrected orientation of struct + outputs. ### Alphabetical list of new functions added in Octave 8
--- a/scripts/ode/ode23.m Wed Oct 05 12:03:52 2022 -0700 +++ b/scripts/ode/ode23.m Wed Oct 05 16:53:01 2022 -0400 @@ -196,8 +196,8 @@ odeopts.InitialStep = odeopts.direction * ... starting_stepsize (order, fcn, trange(1), init, odeopts.AbsTol, odeopts.RelTol, - strcmpi (odeopts.NormControl, "on"), - odeopts.funarguments); + strcmpi (odeopts.NormControl, + "on"), odeopts.funarguments); endif if (! isempty (odeopts.Mass)) @@ -229,12 +229,7 @@ endif endif - if (nargout == 1) - ## Single output requires auto-selected intermediate times, - ## which is obtained by NOT specifying specific solution times. - trange = [trange(1); trange(end)]; - odeopts.Refine = []; # disable Refine when single output requested - elseif (numel (trange) > 2) + if (numel (trange) > 2) odeopts.Refine = []; # disable Refine when specific times requested endif @@ -246,8 +241,9 @@ feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:}); endif if (! isempty (odeopts.Events)) # Cleanup event function handling - ode_event_handler (odeopts.Events, solution.t(end), - solution.x(end,:).', "done", odeopts.funarguments{:}); + ode_event_handler (odeopts.Events, solution.ode_t(end), ... + solution.ode_x(end,:).', "done", ... + odeopts.funarguments{:}); endif ## Print additional information if option Stats is set @@ -265,16 +261,16 @@ endif if (nargout == 2) - varargout{1} = solution.t; # Time stamps are first output argument - varargout{2} = solution.x; # Results are second output argument + varargout{1} = solution.output_t; # Time stamps are first output argument + varargout{2} = solution.output_x; # Results are second output argument elseif (nargout == 1) - varargout{1}.x = solution.t.'; # Time stamps saved in field x (row vector) - varargout{1}.y = solution.x.'; # Results are saved in field y (row vector) - varargout{1}.solver = solver; # Solver name is saved in field solver + varargout{1}.x = solution.ode_t.'; #Time stamps saved in field x (row vect.) + varargout{1}.y = solution.ode_x.'; #Results are saved in field y (row vect.) + varargout{1}.solver = solver; # Solver name is saved in field solver if (! isempty (odeopts.Events)) - varargout{1}.xe = solution.event{3}; # Time info when an event occurred - varargout{1}.ye = solution.event{4}; # Results when an event occurred - varargout{1}.ie = solution.event{2}; # Index info which event occurred + varargout{1}.xe = solution.event{3}.'; # Time info when an event occurred + varargout{1}.ye = solution.event{4}.'; # Results when an event occurred + varargout{1}.ie = solution.event{2}.'; # Index info which event occurred endif if (strcmpi (odeopts.Stats, "on")) varargout{1}.stats = struct (); @@ -287,8 +283,8 @@ endif elseif (nargout > 2) varargout = cell (1,5); - varargout{1} = solution.t; - varargout{2} = solution.x; + varargout{1} = solution.output_t; + varargout{2} = solution.output_x; if (! isempty (odeopts.Events)) varargout{3} = solution.event{3}; # Time info when an event occurred varargout{4} = solution.event{4}; # Results when an event occurred @@ -355,7 +351,8 @@ %! error ('fout: step "init"'); %! endif %! elseif (isempty (flag)) -%! if (! isequal (size (t), [1, 1])) +%! # Multiple steps can be sent in one function call +%! if (! isequal ( size (t), size (y))) %! error ('fout: step "calc"'); %! endif %! elseif (strcmp (flag, "done")) @@ -370,6 +367,14 @@ %!test # two output arguments %! [t, y] = ode23 (@fpol, [0 2], [2 0]); %! assert ([t(end), y(end,:)], [2, fref], 1e-3); +%!test # correct number of steps with Refine +%! [t1, y1] = ode23 (@fpol, [0 2], [2 0], odeset ("Refine", 1)); +%! [t2, y2] = ode23 (@fpol, [0 2], [2 0], odeset ("Refine", 4)); +%! [t3, y3] = ode23 (@fpol, [0 2], [2 0]); #default Refine=1 +%! s = ode23 (@fpol, [0 2], [2 0], odeset ("Refine", 4)); +%! assert (length (t1) == length (t3)); +%! assert (length (t2) == 4*length (t1) - 3); +%! assert (length (s.x) == length (t1)); %!test # anonymous function instead of real function %! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; %! [t, y] = ode23 (fvdp, [0 2], [2 0]); @@ -435,8 +440,8 @@ %! opt = odeset ("NonNegative", 2); %! sol = ode23 (@fpol, [0 2], [2 0], opt); %! assert ([sol.x(end); sol.y(:,end)], [2; 2; 0], 1e-1); -%!test # Details of OutputSel and Refine can't be tested -%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); +%!test # Details of OutputSel can't be tested +%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1); %! sol = ode23 (@fpol, [0 2], [2 0], opt); %!test # Stats must add further elements in sol %! opt = odeset ("Stats", "on"); @@ -452,13 +457,11 @@ %! assert (isfield (sol, "xe")); %! assert (isfield (sol, "ye")); %!test # Events option, now stop integration -%! warning ("off", "integrate_adaptive:unexpected_termination", "local"); %! opt = odeset ("Events", @fevn, "NormControl", "on"); %! sol = ode23 (@fpol, [0 10], [2 0], opt); -%! assert ([sol.ie, sol.xe, sol.ye], +%! assert ([sol.ie, sol.xe, sol.ye.'], %! [2.0, 2.496110, -0.830550, -2.677589], .5e-1); %!test # Events option, five output arguments -%! warning ("off", "integrate_adaptive:unexpected_termination", "local"); %! opt = odeset ("Events", @fevn, "NormControl", "on"); %! [t, y, vxe, ye, vie] = ode23 (@fpol, [0 10], [2 0], opt); %! assert ([vie, vxe, ye], [2.0, 2.496110, -0.830550, -2.677589], 1e-1); @@ -501,6 +504,12 @@ %! [x, y] = ode23 (@(x,y) 1, [0 1], 1i); %! assert (imag (y), ones (size (y))); +## FIXME: convert to demo or a visible=off test with failable assert/error +## statemments +##%!test # Make sure odeplot works (default OutputFcn when no return value) +##%! ode23 (@fpol, [0 2], [2 0]); +##%! close all + ## Test input validation %!error <Invalid call> ode23 () %!error <Invalid call> ode23 (1)
--- a/scripts/ode/ode23s.m Wed Oct 05 12:03:52 2022 -0700 +++ b/scripts/ode/ode23s.m Wed Oct 05 16:53:01 2022 -0400 @@ -205,12 +205,7 @@ ## Starting the initialization of the core solver ode23s - if (nargout == 1) - ## Single output requires auto-selected intermediate times, - ## which is obtained by NOT specifying specific solution times. - trange = [trange(1); trange(end)]; - odeopts.Refine = []; # disable Refine when single output requested - elseif (numel (trange) > 2) + if (numel (trange) > 2) odeopts.Refine = []; # disable Refine when specific times requested endif @@ -222,8 +217,9 @@ feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:}); endif if (! isempty (odeopts.Events)) # Cleanup event function handling - ode_event_handler (odeopts.Events, solution.t(end), ... - solution.x(end,:).', "done", odeopts.funarguments{:}); + ode_event_handler (odeopts.Events, solution.ode_t(end), ... + solution.ode_x(end,:).', "done", ... + odeopts.funarguments{:}); endif ## Print additional information if option Stats is set @@ -241,16 +237,16 @@ endif if (nargout == 2) - varargout{1} = solution.t; # Time stamps are first output argument - varargout{2} = solution.x; # Results are second output argument + varargout{1} = solution.output_t; # Time stamps are first output argument + varargout{2} = solution.output_x; # Results are second output argument elseif (nargout == 1) - varargout{1}.x = solution.t.'; # Time stamps saved in field x (row vector) - varargout{1}.y = solution.x.'; # Results are saved in field y (row vector) + varargout{1}.x = solution.ode_t.'; #Time stamps saved in field x (row vect.) + varargout{1}.y = solution.ode_x.'; #Results are saved in field y (row vect.) varargout{1}.solver = solver; # Solver name is saved in field solver if (! isempty (odeopts.Events)) - varargout{1}.xe = solution.event{3}; # Time info when an event occurred - varargout{1}.ye = solution.event{4}; # Results when an event occurred - varargout{1}.ie = solution.event{2}; # Index info which event occurred + varargout{1}.xe = solution.event{3}.'; # Time info when an event occurred + varargout{1}.ye = solution.event{4}.'; # Results when an event occurred + varargout{1}.ie = solution.event{2}.'; # Index info which event occurred endif if (strcmpi (odeopts.Stats, "on")) varargout{1}.stats = struct (); @@ -263,8 +259,8 @@ endif elseif (nargout == 5) varargout = cell (1,5); - varargout{1} = solution.t; - varargout{2} = solution.x; + varargout{1} = solution.output_t; + varargout{2} = solution.output_x; if (! isempty (odeopts.Events)) varargout{3} = solution.event{3}; # Time info when an event occurred varargout{4} = solution.event{4}; # Results when an event occurred @@ -389,7 +385,8 @@ %! error ('fout: step "init"'); %! endif %! elseif (isempty (flag)) -%! if (! isequal (size (t), [1, 1])) +%! # Multiple steps can be sent in one function call +%! if (! isequal ( size (t), size (y))) %! error ('fout: step "calc"'); %! endif %! elseif (strcmp (flag, "done")) @@ -404,6 +401,14 @@ %!test # two output arguments %! [t, y] = ode23s (@fpol, [0 2], [2 0]); %! assert ([t(end), y(end,:)], [2, fref], 1e-3); +%!test # correct number of steps with Refine +%! [t1, y1] = ode23s (@fpol, [0 2], [2 0], odeset ("Refine", 1)); +%! [t2, y2] = ode23s (@fpol, [0 2], [2 0], odeset ("Refine", 4)); +%! [t3, y3] = ode23s (@fpol, [0 2], [2 0]); #default Refine=1 +%! s = ode23s (@fpol, [0 2], [2 0], odeset ("Refine", 4)); +%! assert (length (t1) == length (t3)); +%! assert (length (t2) == 4*length (t1) - 3); +%! assert (length (s.x) == length (t1)); %!test # anonymous function instead of real function %! fvdp = @(t,y) [y(2); 10 * (1 - y(1)^2) * y(2) - y(1)]; %! [t, y] = ode23s (fvdp, [0 2], [2 0]); @@ -435,8 +440,8 @@ %! opt = odeset ("RelTol", 1e-8, "NormControl", "on"); %! sol = ode23s (@fpol, [0 2], [2 0], opt); %! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-4); -%!test # Details of OutputSel and Refine can't be tested -%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); +%!test # Details of OutputSel can't be tested +%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1); %! sol = ode23s (@fpol, [0 2], [2 0], opt); %!test # Stats must add further elements in sol %! opt = odeset ("Stats", "on"); @@ -499,6 +504,12 @@ %! [x, y] = ode23s (@(x,y) 1, [0 1], 1i); %! assert (imag (y), ones (size (y))); +## FIXME: convert to demo or a visible=off test with failable assert/error +## statemments +##%!test # Make sure odeplot works (default OutputFcn when no return value) +##%! ode23s (@fpol, [0 2], [2 0]); +##%! close all + ## Test input validation %!error <Invalid call> ode23s () %!error <Invalid call> ode23s (1)
--- a/scripts/ode/ode45.m Wed Oct 05 12:03:52 2022 -0700 +++ b/scripts/ode/ode45.m Wed Oct 05 16:53:01 2022 -0400 @@ -156,7 +156,6 @@ [defaults, classes, attributes] = odedefaults (numel (init), trange(1), trange(end)); - ## FIXME: Refine is not correctly implemented yet defaults = odeset (defaults, "Refine", 4); persistent ode45_ignore_options = ... @@ -196,8 +195,8 @@ odeopts.InitialStep = odeopts.direction * ... starting_stepsize (order, fcn, trange(1), init, odeopts.AbsTol, odeopts.RelTol, - strcmpi (odeopts.NormControl, "on"), - odeopts.funarguments); + strcmpi (odeopts.NormControl, + "on"), odeopts.funarguments); endif if (! isempty (odeopts.Mass)) @@ -229,12 +228,7 @@ endif endif - if (nargout == 1) - ## Single output requires auto-selected intermediate times, - ## which is obtained by NOT specifying specific solution times. - trange = [trange(1); trange(end)]; - odeopts.Refine = []; # disable Refine when single output requested - elseif (numel (trange) > 2) + if (numel (trange) > 2) odeopts.Refine = []; # disable Refine when specific times requested endif @@ -246,8 +240,9 @@ feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:}); endif if (! isempty (odeopts.Events)) # Cleanup event function handling - ode_event_handler (odeopts.Events, solution.t(end), - solution.x(end,:).', "done", odeopts.funarguments{:}); + ode_event_handler (odeopts.Events, solution.ode_t(end), ... + solution.ode_x(end,:).', "done", ... + odeopts.funarguments{:}); endif ## Print additional information if option Stats is set @@ -265,16 +260,16 @@ endif if (nargout == 2) - varargout{1} = solution.t; # Time stamps are first output argument - varargout{2} = solution.x; # Results are second output argument + varargout{1} = solution.output_t; # Time stamps are first output argument + varargout{2} = solution.output_x; # Results are second output argument elseif (nargout == 1) - varargout{1}.x = solution.t.'; # Time stamps saved in field x (row vector) - varargout{1}.y = solution.x.'; # Results are saved in field y (row vector) - varargout{1}.solver = solver; # Solver name is saved in field solver + varargout{1}.x = solution.ode_t.'; #Time stamps saved in field x (row vect.) + varargout{1}.y = solution.ode_x.'; #Results are saved in field y (row vect.) + varargout{1}.solver = solver; # Solver name is saved in field solver if (! isempty (odeopts.Events)) - varargout{1}.xe = solution.event{3}; # Time info when an event occurred - varargout{1}.ye = solution.event{4}; # Results when an event occurred - varargout{1}.ie = solution.event{2}; # Index info which event occurred + varargout{1}.xe = solution.event{3}.'; # Time info when an event occurred + varargout{1}.ye = solution.event{4}.'; # Results when an event occurred + varargout{1}.ie = solution.event{2}.'; # Index info which event occurred endif if (strcmpi (odeopts.Stats, "on")) varargout{1}.stats = struct (); @@ -287,8 +282,8 @@ endif elseif (nargout > 2) varargout = cell (1,5); - varargout{1} = solution.t; - varargout{2} = solution.x; + varargout{1} = solution.output_t; + varargout{2} = solution.output_x; if (! isempty (odeopts.Events)) varargout{3} = solution.event{3}; # Time info when an event occurred varargout{4} = solution.event{4}; # Results when an event occurred @@ -355,7 +350,8 @@ %! error ('fout: step "init"'); %! endif %! elseif (isempty (flag)) -%! if (! isequal (size (t), [1, 1])) +%! # Multiple steps can be sent in one function call +%! if (! isequal ( size (t), size (y))) %! error ('fout: step "calc"'); %! endif %! elseif (strcmp (flag, "done")) @@ -371,8 +367,16 @@ %! [t, y] = ode45 (@fpol, [0 2], [2 0]); %! assert ([t(end), y(end,:)], [2, fref], 1e-2); %!test # not too many steps -%! [t, y] = ode45 (@fpol, [0 2], [2 0]); +%! [t, y] = ode45 (@fpol, [0 2], [2 0], odeset("Refine",1)); %! assert (size (t) < 20); +%!test # correct number of steps with Refine +%! [t1, y1] = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 1)); +%! [t2, y2] = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 4)); +%! [t3, y3] = ode45 (@fpol, [0 2], [2 0]); #default Refine=4 +%! s = ode45 (@fpol, [0 2], [2 0], odeset ("Refine", 4)); +%! assert (length (t2) == length (t3)); +%! assert (length (t2) == 4*length (t1) - 3); +%! assert (length (s.x) == length (t1)); %!test # anonymous function instead of real function %! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; %! [t, y] = ode45 (fvdp, [0 2], [2 0]); @@ -392,7 +396,7 @@ %! [t, y] = ode45 (@(t,y) y, [-2 0], 2); %! assert ([t(end), y(end,:)], vref, 1e-1); %!test # InitialStep option -%! opt = odeset ("InitialStep", 1e-8); +%! opt = odeset ("InitialStep", 1e-8, "Refine", 1); %! [t, y] = ode45 (@fpol, [0 0.2], [2 0], opt); %! assert ([t(2)-t(1)], [1e-8], 1e-9); %!test # MaxStep option @@ -447,8 +451,8 @@ %! opt = odeset ("NonNegative", 2); %! sol = ode45 (@fpol, [0 2], [2 0], opt); %! assert ([sol.x(end); sol.y(:,end)], [2; 2; 0], 0.5); -%!test # Details of OutputSel and Refine can't be tested -%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); +%!test # Details of OutputSel can't be tested +%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1); %! sol = ode45 (@fpol, [0 2], [2 0], opt); %!test # Stats must add further elements in sol %! opt = odeset ("Stats", "on"); @@ -464,13 +468,11 @@ %! assert (isfield (sol, "xe")); %! assert (isfield (sol, "ye")); %!test # Events option, now stop integration -%! warning ("off", "integrate_adaptive:unexpected_termination", "local"); %! opt = odeset ("Events", @fevn, "NormControl", "on"); %! sol = ode45 (@fpol, [0 10], [2 0], opt); -%! assert ([sol.ie, sol.xe, sol.ye], +%! assert ([sol.ie, sol.xe, sol.ye.'], %! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); %!test # Events option, five output arguments -%! warning ("off", "integrate_adaptive:unexpected_termination", "local"); %! opt = odeset ("Events", @fevn, "NormControl", "on"); %! [t, y, vxe, ye, vie] = ode45 (@fpol, [0 10], [2 0], opt); %! assert ([vie, vxe, ye], @@ -511,9 +513,15 @@ %!test # Check that imaginary part of solution does not get inverted %! sol = ode45 (@(x,y) 1, [0 1], 1i); %! assert (imag (sol.y), ones (size (sol.y))); -%! [x, y] = ode45 (@(x,y) 1, [0 1], 1i); +%! [x, y] = ode45 (@(x,y) 1, [0 1], 1i, odeset ("Refine", 1)); %! assert (imag (y), ones (size (y))); +## FIXME: convert to demo or a visible=off test with failable assert/error +## statemments +##%!test # Make sure odeplot works (default OutputFcn when no return value) +##%! ode45 (@fpol, [0 2], [2 0]); +##%! close all + %!error <Invalid call> ode45 () %!error <Invalid call> ode45 (1) %!error <Invalid call> ode45 (1,2)
--- a/scripts/ode/odeplot.m Wed Oct 05 12:03:52 2022 -0700 +++ b/scripts/ode/odeplot.m Wed Oct 05 16:53:01 2022 -0400 @@ -40,8 +40,8 @@ ## contains the initial conditions for the ode problem (@var{y0}). ## ## @item @qcode{""} -## The input @var{t} must be a scalar double specifying the time for which -## the solution in input @var{y} was calculated. +## The input @var{t} must be a scalar double or vector specifying the time(s) +## for which the solution in input @var{y} was calculated. ## ## @item @qcode{"done"} ## The inputs should be empty, but are ignored if they are present. @@ -120,3 +120,15 @@ %! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; %! opt = odeset ("OutputFcn", @odeplot, "RelTol", 1e-6); %! sol = ode45 (fvdp, [0 20], [2 0], opt); + +## FIXME: convert to demo or a visible=off test with failable assert/error +## statemments +## Test that the function works for expected ode OutputFcn calls +## %!test +## %! t = linspace(0,2,10); +## %! y = [0.2*t; -0.1*t.^2-1; sin(t)]; +## %! odeplot ([0 2], y(:,1), "init"); +## %! odeplot (t(2), y(:,2), []); +## %! odeplot (t(3:end), y(:, 3:end), []); +## %! odeplot ([], [], "done"); +## %! close all
--- a/scripts/ode/odeset.m Wed Oct 05 12:03:52 2022 -0700 +++ b/scripts/ode/odeset.m Wed Oct 05 16:53:01 2022 -0400 @@ -130,7 +130,6 @@ ## time step or also at intermediate time instances. The value should be ## a scalar indicating the number of equally spaced time points to use ## within each timestep at which to return output. -## @emph{Note}: This option is not yet implemented. ## ## @item @code{RelTol}: positive scalar ## Relative error tolerance.
--- a/scripts/ode/private/integrate_adaptive.m Wed Oct 05 12:03:52 2022 -0700 +++ b/scripts/ode/private/integrate_adaptive.m Wed Oct 05 16:53:01 2022 -0400 @@ -73,13 +73,13 @@ fixed_times = numel (tspan) > 2; - t_new = t_old = t = tspan(1); - x_new = x_old = x = x0(:); + t_new = t_old = ode_t = output_t = tspan(1); + x_new = x_old = ode_x = output_x = x0(:); ## Get first initial timestep dt = options.InitialStep; if (isempty (dt)) - dt = starting_stepsize (order, fcn, t, x, + dt = starting_stepsize (order, fcn, ode_t, ode_x, options.AbsTol, options.RelTol, strcmp (options.NormControl, "on"), options.funarguments); @@ -95,12 +95,23 @@ facmax = 1.5; fac = 0.38^(1/(order+1)); # formula taken from Hairer + ## Initialize Refine value + refine = options.Refine; + if isempty (refine) + refine = 1; + elseif ((refine != round (refine)) || (refine < 1)) + refine = 1; + warning ("integrate_adaptive:invalid_refine", + ["Invalid value of Refine. Refine must be a positive " ... + "integer. Setting Refine = 1."] ); + endif + ## Initialize the OutputFcn if (options.haveoutputfunction) if (! isempty (options.OutputSel)) - solution.retout = x(options.OutputSel,end); + solution.retout = output_x(options.OutputSel, end); else - solution.retout = x; + solution.retout = output_x; endif feval (options.OutputFcn, tspan, solution.retout, "init", options.funarguments{:}); @@ -110,7 +121,7 @@ have_EventFcn = false; if (! isempty (options.Events)) have_EventFcn = true; - ode_event_handler (options.Events, tspan(1), x, + ode_event_handler (options.Events, tspan(1), ode_x, "init", options.funarguments{:}); endif @@ -150,130 +161,84 @@ solution.cntloop += 1; ireject = 0; # Clear reject counter - - ## if output time steps are fixed - if (fixed_times) - - t_caught = find ((dir * tspan(iout:end) > dir * t_old) - & (dir * tspan(iout:end) <= dir * t_new)); - t_caught = t_caught + iout - 1; + terminal_event = false; + terminal_output = false; - if (! isempty (t_caught)) - t(t_caught) = tspan(t_caught); - iout = max (t_caught); - x(:, t_caught) = ... - runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], - t(t_caught), new_k_vals, dt, fcn, - options.funarguments); - - istep += 1; + istep++; + ode_t(istep) = t_new; + ode_x(:, istep) = x_new; + iadd = 0; # Number of output points added this iteration - ## Call Events function only if a valid result has been found. - ## Stop integration if eventbreak is true. - if (have_EventFcn) - break_loop = false; - for idenseout = 1:numel (t_caught) - id = t_caught(idenseout); - td = t(id); - solution.event = ... - ode_event_handler (options.Events, t(id), x(:, id), [], - options.funarguments{:}); - if (! isempty (solution.event{1}) && solution.event{1} == 1) - t(id) = solution.event{3}(end); - t = t(1:id); - x(:, id) = solution.event{4}(end, :).'; - x = x(:,1:id); - solution.unhandledtermination = false; - break_loop = true; - break; - endif - endfor - if (break_loop) - break; - endif - endif + ## Check for Events + if (have_EventFcn) + solution.event = ode_event_handler (options.Events, t_new, x_new, ... + [], options.funarguments{:}); + ## Check for terminal Event + if (! isempty (solution.event{1}) && solution.event{1} == 1) + ode_t(istep) = solution.event{3}(end); + ode_x(:, istep) = solution.event{4}(end, :).'; + solution.unhandledtermination = false; + terminal_event = true; + endif + endif - ## Call OutputFcn only if a valid result has been found. - ## Stop integration if function returns true. - if (options.haveoutputfunction) - cnt = options.Refine + 1; - approxtime = linspace (t_old, t_new, cnt); - approxvals = interp1 ([t_old, t(t_caught), t_new], - [x_old, x(:, t_caught), x_new] .', - approxtime, "linear") .'; - if (isvector (approxvals) && ! isscalar (approxtime)) - approxvals = approxvals.'; - endif - if (! isempty (options.OutputSel)) - approxvals = approxvals(options.OutputSel, :); - endif - stop_solve = false; - for ii = 1:numel (approxtime) - stop_solve = feval (options.OutputFcn, - approxtime(ii), approxvals(:, ii), [], - options.funarguments{:}); - if (stop_solve) - break; # break from inner loop - endif - endfor - if (stop_solve) # Leave main loop - solution.unhandledtermination = false; - break; - endif - endif - + ## Interpolate to specified or Refined points + if (fixed_times) + t_caught = find ((dir * tspan(iout:end) > dir * t_old) ... + & (dir * tspan(iout:end) <= dir * ode_t(istep))); + t_caught = t_caught + iout - 1; + iadd = length (t_caught); + if (! isempty (t_caught)) + output_t(t_caught) = tspan(t_caught); + iout = max (t_caught); + output_x(:, t_caught) = ... + runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], ... + output_t(t_caught), new_k_vals, dt, ... + fcn, options.funarguments); endif - - else # not fixed times - - t(++istep) = t_new; - x(:, istep) = x_new; - iout = istep; + ## Add a possible additional output value if we found a terminal Event + if ((terminal_event == true) && ... + (dir * ode_t(istep) > dir * output_t(iout))) + iadd += 1; + iout += 1; + output_x(:, iout) = ode_x(:, istep); + output_t(iout) = ode_t(istep); + endif + elseif (refine > 1) + iadd = refine; + tadd = linspace (t_old, ode_t(istep), refine + 1); + tadd = tadd(2:end); + output_x(:, iout + (1:iadd)) = ... + runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], ... + tadd, new_k_vals, dt, fcn, ... + options.funarguments); + output_t(iout + (1:iadd)) = tadd; + iout = length (output_t); + else # refine = 1 + iadd = 1; + iout += iadd; + output_x(:, iout) = ode_x(:, istep); + output_t(iout) = ode_t(istep); + end - ## Call Events function only if a valid result has been found. - ## Stop integration if eventbreak is true. - if (have_EventFcn) - solution.event = ... - ode_event_handler (options.Events, t(istep), x(:, istep), [], - options.funarguments{:}); - if (! isempty (solution.event{1}) && solution.event{1} == 1) - t(istep) = solution.event{3}(end); - x(:, istep) = solution.event{4}(end, :).'; - solution.unhandledtermination = false; - break; - endif + ## Call OutputFcn + if ((options.haveoutputfunction) && (iadd > 0)) + xadd = output_x(:, (iout-iadd+1):end); + tadd = output_t((iout-iadd+1):end); + if (! isempty (options.OutputSel)) + xadd = xadd(options.OutputSel, :); endif + stop_solve = feval (options.OutputFcn, tadd, xadd, ... + [], options.funarguments{:}); + if (stop_solve) + solution.unhandledtermination = false; + terminal_output = true; + endif + endif + if (terminal_event || terminal_output) + break; # break from main loop + endif - ## Call OutputFcn only if a valid result has been found. - ## Stop integration if function returns true. - if (options.haveoutputfunction) - cnt = options.Refine + 1; - approxtime = linspace (t_old, t_new, cnt); - approxvals = interp1 ([t_old, t_new], - [x_old, x_new] .', - approxtime, "linear") .'; - if (isvector (approxvals) && ! isscalar (approxtime)) - approxvals = approxvals.'; - endif - if (! isempty (options.OutputSel)) - approxvals = approxvals(options.OutputSel, :); - endif - stop_solve = false; - for ii = 1:numel (approxtime) - stop_solve = feval (options.OutputFcn, - approxtime(ii), approxvals(:, ii), [], - options.funarguments{:}); - if (stop_solve) - break; # break from inner loop - endif - endfor - if (stop_solve) # Leave main loop - solution.unhandledtermination = false; - break; - endif - endif - - endif ## move to next time-step t_old = t_new; @@ -303,7 +268,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 (ode_t(end)))) break; endif @@ -313,7 +278,7 @@ endwhile ## Check if integration of the ode has been successful - if (dir * t(end) < dir * tspan(end)) + if (dir * ode_t(end) < dir * tspan(end)) if (solution.unhandledtermination == true) warning ("integrate_adaptive:unexpected_termination", [" Solving was not successful. ", ... @@ -322,20 +287,14 @@ " This may happen if the stepsize becomes too small. ", ... " Try to reduce the value of 'InitialStep'", ... " and/or 'MaxStep' with the command 'odeset'."], - t(end), tspan(end)); - else - warning ("integrate_adaptive:unexpected_termination", - ["Solver was stopped by a call of 'break'", ... - " in the main iteration loop at time", ... - " t = %f before the endpoint at tend = %f was reached. ", ... - " This may happen because the @odeplot function", ... - " returned 'true' or the @event function returned 'true'."], - t(end), tspan(end)); + ode_t(end), tspan(end)); endif endif ## Set up return structure - solution.t = t(:); - solution.x = x.'; + solution.ode_t = ode_t(:); + solution.ode_x = ode_x.'; + solution.output_t = output_t(:); + solution.output_x = output_x.'; endfunction
--- a/scripts/ode/private/ode_event_handler.m Wed Oct 05 12:03:52 2022 -0700 +++ b/scripts/ode/private/ode_event_handler.m Wed Oct 05 16:53:01 2022 -0400 @@ -115,16 +115,47 @@ else retcell{1} = any (term(idx)); # Stop integration or not endif - idx = idx(1); # Use first event found if there are multiple. - retcell{2}(evtcnt,1) = idx; - ## Calculate the time stamp when the event function returned 0 and - ## calculate new values for the integration results, we do both by - ## a linear interpolation. - tnew = t - evt(idx) * (t - told) / (evt(idx) - evtold(idx)); - ynew = (y - (t - tnew) * (y - yold) / (t - told)).'; - retcell{3}(evtcnt,1) = tnew; - retcell{4}(evtcnt,:) = ynew; - evtcnt += 1; + evtcntnew = 1; + ## Add all events this step to the output. + for idx2 = idx # Loop through all values of idx + ## Calculate the time stamp when the event function returned 0 and + ## calculate new values for the integration results, we do both by + ## a linear interpolation. + tnew = t - evt(idx2) * (t - told) / (evt(idx2) - evtold(idx2)); + ynew = (y - (t - tnew) * (y - yold) / (t - told)).'; + tnews(evtcntnew, 1) = tnew; + ynews(evtcntnew, :) = ynew; + terms(evtcntnew, 1) = term(idx2); + evtcntnew += 1; + endfor + ## Sort by time of event + if length (idx) > 1 + [tnews, idx_sort] = sort (tnews, "ascend"); + idxs = idx(idx_sort); + ynews = ynews(idx_sort,:); + terms = terms(idx_sort); + else + idxs = idx; + endif + ## Check for terminal events and remove any events after terminal. + ## Any events at same time as first terminal event will be retained. + idx3 = find (terms, 1); # Find first terminal event by time + if ! isempty (idx3) + t_cutoff = tnews(idx3); + ## Last index to return + evtcntnew = find (tnews == t_cutoff, 1, "last"); + else + evtcntnew = length (terms); # Return all indices if no terminal + endif + idxs = idxs(1:evtcntnew); + tnews = tnews(1:evtcntnew); + + ## Populate return values with sorted, clipped values + evtcntrange = evtcnt - 1 + (1:evtcntnew); + evtcnt += evtcntnew; + retcell{2}(evtcntrange, 1) = idxs(:); + retcell{3}(evtcntrange, 1) = tnews(:); + retcell{4}(evtcntrange, :) = ynews(1:evtcntnew,:); endif endif