Mercurial > octave
changeset 31409:88fff8521d76
Improve "Events" location accuracy in ode45, ode23, ode23s (bug #63162)
* ode_event_handler.m: Find root using polynomial interpolation function.
Remove unused y cell input code. Update documentation. Update function input
and persistent variables. Add helper function "evtfcn_val" for use with root
finder.
* ode23.m, ode23s.m, ode45m: Update self tests, add FIXME comments, modify
ode_event_handler function call to match above changes.
* runge_kutta_interpolate.m: Remove unneeded second and higher-order
interpolation code, add explanatory comments, change bahavior for invalid
input.
* integrate_adaptive.m: Modify function calls. Fix function calls to
runge_kutta_interpolate to match above changes.
| author | Nicholas R. Jankowski <jankowski.nicholas@gmail.com> |
|---|---|
| date | Thu, 10 Nov 2022 19:17:35 -0500 |
| parents | 870036573716 |
| children | 4114ec68f10a |
| files | scripts/ode/ode23.m scripts/ode/ode23s.m scripts/ode/ode45.m scripts/ode/private/integrate_adaptive.m scripts/ode/private/ode_event_handler.m scripts/ode/private/runge_kutta_interpolate.m |
| diffstat | 6 files changed, 88 insertions(+), 73 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ode/ode23.m Mon Nov 07 16:50:39 2022 -0500 +++ b/scripts/ode/ode23.m Thu Nov 10 19:17:35 2022 -0500 @@ -154,6 +154,10 @@ "ode23: FCN must be a valid function handle"); endif + ## FIXME: Warn user if ! isempty (funarguments) + ## Not a documented behavior and may be deprecated + + ## Start preprocessing, have a look which options are set in odeopts, ## check if an invalid or unused option is set. [defaults, classes, attributes] = odedefaults (numel (init), @@ -241,9 +245,7 @@ feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:}); endif if (! isempty (odeopts.Events)) # Cleanup event function handling - ode_event_handler (odeopts.Events, solution.ode_t(end), ... - solution.ode_x(end,:).', "done", ... - odeopts.funarguments{:}); + ode_event_handler ([], [], [], [], [], "done"); endif ## Print additional information if option Stats is set @@ -460,11 +462,11 @@ %! opt = odeset ("Events", @fevn, "NormControl", "on"); %! sol = ode23 (@fpol, [0 10], [2 0], opt); %! assert ([sol.ie, sol.xe, sol.ye.'], -%! [2.0, 2.496110, -0.830550, -2.677589], .5e-1); +%! [2.0, 2.496110, -0.830550, -2.677589], 1e-3); %!test # Events option, five output arguments %! 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); +%! assert ([vie, vxe, ye], [2.0, 2.496110, -0.830550, -2.677589], 1e-3); %!test # Mass option as function %! opt = odeset ("Mass", @fmas); %! sol = ode23 (@fpol, [0 2], [2 0], opt);
--- a/scripts/ode/ode23s.m Mon Nov 07 16:50:39 2022 -0500 +++ b/scripts/ode/ode23s.m Thu Nov 10 19:17:35 2022 -0500 @@ -155,6 +155,10 @@ "ode23s: FCN must be a valid function handle"); endif + ## FIXME: Warn user if ! isempty (funarguments) + ## Not a documented behavior and may be deprecated + + ## Start preprocessing, have a look which options are set in odeopts, ## check if an invalid or unused option is set. [defaults, classes, attributes] = odedefaults (numel (init), @@ -217,9 +221,7 @@ feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:}); endif if (! isempty (odeopts.Events)) # Cleanup event function handling - ode_event_handler (odeopts.Events, solution.ode_t(end), ... - solution.ode_x(end,:).', "done", ... - odeopts.funarguments{:}); + ode_event_handler ([], [], [], [], [], "done"); endif ## Print additional information if option Stats is set @@ -456,6 +458,15 @@ %! assert (sol.ie(1), 2); %! assert (isfield (sol, "xe")); %! assert (isfield (sol, "ye")); +%!test # Events option, now stop integration +%! opt = odeset ("Events", @fevn, "NormControl", "on"); +%! sol = ode23s (@fpol, [0 10], [2 0], opt); +%! assert ([sol.ie, sol.xe, sol.ye.'], +%! [2.0, 9.094439, -0.996480, -14.180147], 2e-2); +%!test # Events option, five output arguments +%! opt = odeset ("Events", @fevn, "NormControl", "on"); +%! [t, y, vxe, ye, vie] = ode23s (@fpol, [0 10], [2 0], opt); +%! assert ([vie, vxe, ye], [2.0, 9.094439, -0.996480, -14.180147], 2e-2); %!test # Mass option as function %! opt = odeset ("Mass", @fmas); %! sol = ode23s (@fpol, [0 2], [2 0], opt);
--- a/scripts/ode/ode45.m Mon Nov 07 16:50:39 2022 -0500 +++ b/scripts/ode/ode45.m Thu Nov 10 19:17:35 2022 -0500 @@ -151,6 +151,10 @@ "ode45: FCN must be a valid function handle"); endif + ## FIXME: Warn user if ! isempty (funarguments) + ## Not a documented behavior and may be deprecated + + ## Start preprocessing, have a look which options are set in odeopts, ## check if an invalid or unused option is set [defaults, classes, attributes] = odedefaults (numel (init), @@ -240,9 +244,7 @@ feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:}); endif if (! isempty (odeopts.Events)) # Cleanup event function handling - ode_event_handler (odeopts.Events, solution.ode_t(end), ... - solution.ode_x(end,:).', "done", ... - odeopts.funarguments{:}); + ode_event_handler ([], [], [], [], [], "done"); endif ## Print additional information if option Stats is set @@ -471,12 +473,12 @@ %! opt = odeset ("Events", @fevn, "NormControl", "on"); %! sol = ode45 (@fpol, [0 10], [2 0], opt); %! assert ([sol.ie, sol.xe, sol.ye.'], -%! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); +%! [2.0, 2.496110, -0.830550, -2.677589], 2e-3); %!test # Events option, five output arguments %! opt = odeset ("Events", @fevn, "NormControl", "on"); %! [t, y, vxe, ye, vie] = ode45 (@fpol, [0 10], [2 0], opt); %! assert ([vie, vxe, ye], -%! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); +%! [2.0, 2.496110, -0.830550, -2.677589], 2e-3); %!test # Mass option as function %! opt = odeset ("Mass", @fmas); %! sol = ode45 (@fpol, [0 2], [2 0], opt);
--- a/scripts/ode/private/integrate_adaptive.m Mon Nov 07 16:50:39 2022 -0500 +++ b/scripts/ode/private/integrate_adaptive.m Thu Nov 10 19:17:35 2022 -0500 @@ -121,8 +121,9 @@ have_EventFcn = false; if (! isempty (options.Events)) have_EventFcn = true; - ode_event_handler (options.Events, tspan(1), ode_x, - "init", options.funarguments{:}); + options.Events = @(t,y) options.Events (t, y, options.funarguments{:}); + ode_event_handler (options.Events, tspan(1), ode_x, ... + [], order, "init"); endif if (options.havenonnegative) @@ -171,8 +172,8 @@ ## Check for Events if (have_EventFcn) - solution.event = ode_event_handler (options.Events, t_new, x_new, ... - [], options.funarguments{:}); + solution.event = ode_event_handler ([], t_new, x_new, ... + new_k_vals, [], []); ## Check for terminal Event if (! isempty (solution.event{1}) && solution.event{1} == 1) ode_t(istep) = solution.event{3}(end); @@ -193,8 +194,7 @@ 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); + output_t(t_caught), new_k_vals); endif ## Add a possible additional output value if we found a terminal Event if ((terminal_event == true) && ... @@ -210,8 +210,7 @@ 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); + tadd, new_k_vals); output_t(iout + (1:iadd)) = tadd; iout = length (output_t); else # refine = 1
--- a/scripts/ode/private/ode_event_handler.m Mon Nov 07 16:50:39 2022 -0500 +++ b/scripts/ode/private/ode_event_handler.m Thu Nov 10 19:17:35 2022 -0500 @@ -24,43 +24,45 @@ ######################################################################## ## -*- texinfo -*- -## @deftypefn {} {@var{retval} =} ode_event_handler (@var{@@evtfcn}, @var{t}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{}) +## @deftypefn {} {@var{retval} =} ode_event_handler (@var{@@evt_fcn}, @var{t}, @var{y}, @var{k_vals}, @var{ord}, @var{flag}) ## -## Return the solution of the event function (@var{@@evtfcn}) which is +## Return the solution of the event function (@var{@@evt_fcn}) 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 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 third input argument @var{y} is a scalar or a column vector of type +## double which specifies the solution(s) at time @var{t}. ## -## The fourth input argument @var{flag} is of type string. Valid values are: +## The fourth input argument @var{k_vals} is a vector or matrix with the +## k values obtained from the most recent integration step. +## +## The fifth input argument @var{ord} is the order of the integration technique. +## +## The sixth input argument @var{flag} is of type string. Valid values are: ## ## @table @option ## @item @qcode{"init"} ## Initialize internal persistent variables of the function ## @code{ode_event_handler} and return an empty cell array of size 4. ## -## @item @qcode{"calc"} -## Evaluate the event function and return the solution @var{retval} as a cell -## array of size 4. +## @item @qcode{""} +## (default) Evaluate the event function and return the solution @var{retval} +## as a cell array of size 4. Inputs @var{ord} and @var{@@evt_fcn} are ignored, +## since they are set in the @qcode{"init"} step. ## ## @item @qcode{"done"} ## Clean up internal variables of the function @code{ode_event_handler} and -## return an empty cell array of size 4. +## return an empty cell array of size 4. All other inputs are ignored with +## this flag. ## @end table ## -## 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 and it should never be ## necessary to call it directly. ## @end deftypefn -function retval = ode_event_handler (evtfcn, t, y, flag = "", varargin) +function retval = ode_event_handler (evt_fcn, t, y, k_vals, ord, flag = "") ## No error handling has been implemented in this function to achieve ## the highest performance possible. @@ -77,23 +79,14 @@ ## 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; + persistent evtold told yold retcell order evtfcn; persistent evtcnt = 1; # Don't remove. Required for Octave parser. persistent firstrun = true; if (isempty (flag)) ## Process the event, i.e., ## find the zero crossings for either a rising or falling edge - if (! iscell (y)) - inpargs = {evtfcn, t, y}; - else - inpargs = {evtfcn, t, y{1}, y{2}}; - y = y{1}; # Delete cell element 2 - endif - if (nargin > 4) - inpargs = {inpargs{:}, varargin{:}}; - endif - [evt, term, dir] = feval (inpargs{:}); + [evt, term, dir] = evtfcn (t, y); ## We require that all return values be row vectors evt = evt(:).'; term = term(:).'; dir = dir(:).'; @@ -119,10 +112,21 @@ ## 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)).'; + ## calculate new values for the integration results. We do both by + ## root solving, calling the Event function with y values from + ## the RK interpolation polynomial. Set tolerance to zero (actually + ## uses machine tolerance based criterion in this case) since we + ## don't have a way for the user to specify an acceptable tolerance. + ## For simple Event functions, this means we're basically root + ## finding on a small interval for a polynomial, so we expect + ## pretty quick convergence. + tvals = [told t]; + yvals = [yold y]; + tnew = fzero (@(t2) evtfcn_val (evtfcn, t2, ... + runge_kutta_interpolate (order, tvals, yvals, ... + t2, k_vals), idx2), tvals, optimset ("TolX", 0)); + ynew = runge_kutta_interpolate (order, tvals, yvals, tnew, k_vals); + tnews(evtcntnew, 1) = tnew; ynews(evtcntnew, :) = ynew; terms(evtcntnew, 1) = term(idx2); @@ -170,17 +174,10 @@ ## a value for evtold. firstrun = true; + order = ord; + evtfcn = evt_fcn; - if (! iscell (y)) - inpargs = {evtfcn, t, y}; - else - inpargs = {evtfcn, t, y{1}, y{2}}; - y = y{1}; # Delete cell element 2 - endif - if (nargin > 4) - inpargs = {inpargs{:}, varargin{:}}; - endif - [evtold, ~, ~] = feval (inpargs{:}); + [evtold, ~, ~] = evtfcn (t, y); ## We require that all return values be row vectors evtold = evtold(:).'; told = t; yold = y; @@ -190,9 +187,15 @@ elseif (strcmp (flag, "done")) ## Clear this event handling function firstrun = true; - evtold = told = yold = evtcnt = []; + evtold = told = yold = evtcnt = order = evtfcn = []; retval = retcell = cell (1,4); endif endfunction + +function val = evtfcn_val (evtfcn, t, y, ind) + [evt, ~, ~] = evtfcn (t, y); + val = evt(ind); +endfunction +
--- a/scripts/ode/private/runge_kutta_interpolate.m Mon Nov 07 16:50:39 2022 -0500 +++ b/scripts/ode/private/runge_kutta_interpolate.m Thu Nov 10 19:17:35 2022 -0500 @@ -23,22 +23,21 @@ ## ######################################################################## -function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, fcn, args) +function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals) switch (order) case 1 + ## Unused u_interp = interp1 (z, u.', t, "linear"); case 2 - if (! isempty (k_vals)) - der = k_vals(:,1); - else - der = feval (fcn, z(1) , u(:,1), args); - endif + ## ode23s with Rosenbrock scheme: + der = k_vals(:,1); u_interp = quadratic_interpolation (z, u, der, t); case 3 + ## ode23 with Bogacki-Shampine scheme: u_interp = hermite_cubic_interpolation (z, u, k_vals, t); case 5 @@ -46,11 +45,10 @@ u_interp = hermite_quartic_interpolation (z, u, k_vals, t); otherwise - warning (["High order interpolation not yet implemented: ", ... - "using cubic interpolation instead"]); - der(:,1) = feval (fcn, z(1), u(:,1), args); - der(:,2) = feval (fcn, z(2), u(:,2), args); - u_interp = hermite_cubic_interpolation (z, u, der, t); + ## This should never happen + warning (["runge_kutta_interpolate: Invalid/unimplemented ", ... + "interpolation order: using linear interpolation instead"]); + u_interp = interp1 (z, u.', t, "linear"); endswitch