Mercurial > octave

annotate 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
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
1 ########################################################################
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
2 ##
30564
796f54d4ddbf update Octave Project Developers copyright for the new year
John W. Eaton <jwe@octave.org>
parents: 29358
diff changeset
3 ## Copyright (C) 2006-2022 The Octave Project Developers
27918
b442ec6dda5c use centralized file for copyright info for individual contributors
John W. Eaton <jwe@octave.org>
parents: 26376
diff changeset
4 ##
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
5 ## See the file COPYRIGHT.md in the top-level directory of this
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
6 ## distribution or <https://octave.org/copyright/>.
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
7 ##
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
8 ## This file is part of Octave.
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
9 ##
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23565
diff changeset
10 ## Octave is free software: you can redistribute it and/or modify it
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
11 ## under the terms of the GNU General Public License as published by
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23565
diff changeset
12 ## the Free Software Foundation, either version 3 of the License, or
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22655
diff changeset
13 ## (at your option) any later version.
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
14 ##
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
15 ## Octave is distributed in the hope that it will be useful, but
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
16 ## WITHOUT ANY WARRANTY; without even the implied warranty of
22755
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22655
diff changeset
17 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
3a2b891d0b33 maint: Standardize Copyright formatting.
Rik <rik@octave.org>
parents: 22655
diff changeset
18 ## GNU General Public License for more details.
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
19 ##
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
20 ## You should have received a copy of the GNU General Public License
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
21 ## along with Octave; see the file COPYING. If not, see
24534
194eb4bd202b maint: Update punctuation for GPL v3 license text.
Rik <rik@octave.org>
parents: 23565
diff changeset
22 ## <https://www.gnu.org/licenses/>.
27923
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
23 ##
bd51beb6205e update formatting of copyright notices
John W. Eaton <jwe@octave.org>
parents: 27919
diff changeset
24 ########################################################################
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
25
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
26 ## -*- texinfo -*-
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
27 ## @deftypefn {} {@var{retval} =} ode_event_handler (@var{@@evtfcn}, @var{t}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{})
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20536
diff changeset
28 ##
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
29 ## Return the solution of the event function (@var{@@evtfcn}) which is
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
30 ## specified in the form of a function handle.
23565
3a730821e4a2 doc: Peridoc grammarcheck of documentation.
Rik <rik@octave.org>
parents: 23220
diff changeset
31 ##
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
32 ## The second input argument @var{t} is a scalar double and specifies the time
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
33 ## of the event evaluation.
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
34 ##
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
35 ## The third input argument @var{y} may be a column vector of type double
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
36 ## (for ODEs and DAEs) which specifies the solutions. Alternatives, @var{y}
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
37 ## may be a cell array (for IDEs and DDEs) which specifies the solutions and
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
38 ## derivatives.
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
39 ##
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
40 ## The fourth input argument @var{flag} is of type string. Valid values are:
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20536
diff changeset
41 ##
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
42 ## @table @option
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20536
diff changeset
43 ## @item @qcode{"init"}
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
44 ## Initialize internal persistent variables of the function
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
45 ## @code{ode_event_handler} and return an empty cell array of size 4.
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20536
diff changeset
46 ##
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20536
diff changeset
47 ## @item @qcode{"calc"}
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
48 ## Evaluate the event function and return the solution @var{retval} as a cell
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
49 ## array of size 4.
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20536
diff changeset
50 ##
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20536
diff changeset
51 ## @item @qcode{"done"}
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
52 ## Clean up internal variables of the function @code{ode_event_handler} and
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
53 ## return an empty cell array of size 4.
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
54 ## @end table
20548
25623ef2ff4f doc: Rewrite docstrings for ode* family of functions.
Rik <rik@octave.org>
parents: 20536
diff changeset
55 ##
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
56 ## If additional input arguments @var{par1}, @var{par2}, @dots{} are given
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
57 ## these parameters are passed directly to the event function.
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
58 ##
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
59 ## This function is an ODE internal helper function and it should never be
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
60 ## necessary to call it directly.
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
61 ## @end deftypefn
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
62
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
63 function retval = ode_event_handler (evtfcn, t, y, flag = "", varargin)
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
64
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
65 ## No error handling has been implemented in this function to achieve
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
66 ## the highest performance possible.
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
67
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
68 ## retval{1} is true (to terminate) or false (to continue)
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
69 ## retval{2} is the index information for which an event occurred
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
70 ## retval{3} is the time information column vector
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
71 ## retval{4} is the line by line result information matrix
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
72
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
73 ## These persistent variables store the results and time value from the
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
74 ## processing in the previous time stamp.
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
75 ## evtold the results from the event function
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
76 ## told the time stamp
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
77 ## yold the ODE result
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
78 ## retcell the return values cell array
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
79 ## evtcnt the counter for how often this function has been called
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
80 persistent evtold told yold retcell;
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
81 persistent evtcnt = 1; # Don't remove. Required for Octave parser.
23079
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
82 persistent firstrun = true;
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
83
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
84 if (isempty (flag))
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
85 ## Process the event, i.e.,
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
86 ## find the zero crossings for either a rising or falling edge
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
87 if (! iscell (y))
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
88 inpargs = {evtfcn, t, y};
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
89 else
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
90 inpargs = {evtfcn, t, y{1}, y{2}};
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
91 y = y{1}; # Delete cell element 2
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
92 endif
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
93 if (nargin > 4)
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
94 inpargs = {inpargs{:}, varargin{:}};
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
95 endif
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
96 [evt, term, dir] = feval (inpargs{:});
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
97
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
98 ## We require that all return values be row vectors
22655
6b134d294d61 ode solvers: use ordinary transpose instead of Hermitian conjugate (bug #49410).
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22648
diff changeset
99 evt = evt(:).'; term = term(:).'; dir = dir(:).';
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
100
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
101 ## Check if one or more signs of the event has changed
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
102 signum = (sign (evtold) != sign (evt));
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
103 if (any (signum)) # One or more values have changed
29193
b2d5ee958d7f ode_event_handler.m: Fix mishandling of event edge types and multiple events (bug #59709).
Rik <rik@octave.org>
parents: 27923
diff changeset
104 ## Find events where either rising and falling edges are counted (dir==0)
b2d5ee958d7f ode_event_handler.m: Fix mishandling of event edge types and multiple events (bug #59709).
Rik <rik@octave.org>
parents: 27923
diff changeset
105 ## or where the specified edge type matches the event edge type.
b2d5ee958d7f ode_event_handler.m: Fix mishandling of event edge types and multiple events (bug #59709).
Rik <rik@octave.org>
parents: 27923
diff changeset
106 idx = signum & (dir == 0 | dir == sign (evt));
b2d5ee958d7f ode_event_handler.m: Fix mishandling of event edge types and multiple events (bug #59709).
Rik <rik@octave.org>
parents: 27923
diff changeset
107 idx = find (idx); # convert logical to numeric index or []
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
108
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
109 ## Create new output values if a valid index has been found
20634
80e630b37ba1 maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
Rik <rik@octave.org>
parents: 20631
diff changeset
110 if (! isempty (idx))
20536
6256f6e366ac Fix copyright text in private ode functions
Carlo de Falco <carlo.defalco@polimi.it>
parents: 20533
diff changeset
111 ## Change the persistent result cell array
23079
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
112 if (firstrun)
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
113 ## Matlab compatibility requires ignoring condition on first run.
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
114 retcell{1} = false;
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
115 else
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
116 retcell{1} = any (term(idx)); # Stop integration or not
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
117 endif
31263
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
118 evtcntnew = 1;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
119 ## Add all events this step to the output.
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
120 for idx2 = idx # Loop through all values of idx
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
121 ## Calculate the time stamp when the event function returned 0 and
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
122 ## calculate new values for the integration results, we do both by
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
123 ## a linear interpolation.
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
124 tnew = t - evt(idx2) * (t - told) / (evt(idx2) - evtold(idx2));
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
125 ynew = (y - (t - tnew) * (y - yold) / (t - told)).';
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
126 tnews(evtcntnew, 1) = tnew;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
127 ynews(evtcntnew, :) = ynew;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
128 terms(evtcntnew, 1) = term(idx2);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
129 evtcntnew += 1;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
130 endfor
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
131 ## Sort by time of event
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
132 if length (idx) > 1
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
133 [tnews, idx_sort] = sort (tnews, "ascend");
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
134 idxs = idx(idx_sort);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
135 ynews = ynews(idx_sort,:);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
136 terms = terms(idx_sort);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
137 else
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
138 idxs = idx;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
139 endif
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
140 ## Check for terminal events and remove any events after terminal.
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
141 ## Any events at same time as first terminal event will be retained.
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
142 idx3 = find (terms, 1); # Find first terminal event by time
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
143 if ! isempty (idx3)
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
144 t_cutoff = tnews(idx3);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
145 ## Last index to return
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
146 evtcntnew = find (tnews == t_cutoff, 1, "last");
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
147 else
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
148 evtcntnew = length (terms); # Return all indices if no terminal
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
149 endif
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
150 idxs = idxs(1:evtcntnew);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
151 tnews = tnews(1:evtcntnew);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
152
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
153 ## Populate return values with sorted, clipped values
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
154 evtcntrange = evtcnt - 1 + (1:evtcntnew);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
155 evtcnt += evtcntnew;
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
156 retcell{2}(evtcntrange, 1) = idxs(:);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
157 retcell{3}(evtcntrange, 1) = tnews(:);
449ed6f427cb ode45/23/23s: Implement Events, OutputFcn, & Refine options (bug #49408 and #63063)
Ken Marek <marek_ka@mercer.edu>
parents: 30893
diff changeset
158 retcell{4}(evtcntrange, :) = ynews(1:evtcntnew,:);
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
159 endif
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
160
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
161 endif
23079
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
162
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
163 firstrun = false;
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
164 evtold = evt; told = t; yold = y;
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
165 retval = retcell;
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
166
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
167 elseif (strcmp (flag, "init"))
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
168 ## Call the event function if an event function has been defined to
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
169 ## initialize the internal variables of the event function and to get
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
170 ## a value for evtold.
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
171
23079
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
172 firstrun = true;
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
173
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
174 if (! iscell (y))
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
175 inpargs = {evtfcn, t, y};
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
176 else
30893
e1788b1a315f maint: Use "fcn" as preferred abbreviation for "function" in m-files.
Rik <rik@octave.org>
parents: 30564
diff changeset
177 inpargs = {evtfcn, t, y{1}, y{2}};
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
178 y = y{1}; # Delete cell element 2
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
179 endif
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
180 if (nargin > 4)
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
181 inpargs = {inpargs{:}, varargin{:}};
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
182 endif
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
183 [evtold, ~, ~] = feval (inpargs{:});
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
184
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
185 ## We require that all return values be row vectors
22655
6b134d294d61 ode solvers: use ordinary transpose instead of Hermitian conjugate (bug #49410).
Carlo de Falco <carlo.defalco@polimi.it>
parents: 22648
diff changeset
186 evtold = evtold(:).'; told = t; yold = y;
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
187 evtcnt = 1;
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
188 retval = retcell = cell (1,4);
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
189
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
190 elseif (strcmp (flag, "done"))
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
191 ## Clear this event handling function
23079
17dc6c7ef427 Ignore ODE event function at t==0 for Matlab compatibility (bug #49919).
Rik <rik@octave.org>
parents: 22755
diff changeset
192 firstrun = true;
22648
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
193 evtold = told = yold = evtcnt = [];
9bc03a3f7a34 ode_event_handler.m: Overhaul function.
Rik <rik@octave.org>
parents: 22626
diff changeset
194 retval = retcell = cell (1,4);
20533
fcb792acab9b Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
jcorno <jacopo.corno@gmail.com>
parents:
diff changeset
195
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
196 endif
20621
b92f8e148936 maint: Continued clean-up of functions in ode/private dir.
Rik <rik@octave.org>
parents: 20554
diff changeset
197
20552
eb9e2d187ed2 maint: Use Octave coding conventions in scripts/ode/private dir.
Rik <rik@octave.org>
parents: 20548
diff changeset
198 endfunction