Mercurial > octave

diff scripts/ode/private/ode_event_handler.m @ 31263:449ed6f427cb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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 e1788b1a315f
children 88fff8521d76
line wrap: on
line diff
--- 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