Mercurial > octave-nkf

comparison scripts/ode/private/odepkg_event_handle.m @ 20570:6256f6e366ac

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fix copyright text in private ode functions * script/ode/private/*.m: Fox the copyright text.
author Carlo de Falco <carlo.defalco@polimi.it>
date Fri, 02 Oct 2015 05:50:43 +0200
parents fcb792acab9b
children 25623ef2ff4f
comparison
equal deleted inserted replaced
20569:b70cc4bd8109 20570:6256f6e366ac
1 %# Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> 1 ## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net>
2 %# OdePkg - A package for solving ordinary differential equations and more 2 ##
3 %# 3 ## This file is part of Octave.
4 %# This program is free software; you can redistribute it and/or modify 4 ##
5 %# it under the terms of the GNU General Public License as published by 5 ## Octave is free software; you can redistribute it and/or modify it
6 %# the Free Software Foundation; either version 2 of the License, or 6 ## under the terms of the GNU General Public License as published by
7 %# (at your option) any later version. 7 ## the Free Software Foundation; either version 3 of the License, or (at
8 %# 8 ## your option) any later version.
9 %# This program is distributed in the hope that it will be useful, 9 ##
10 %# but WITHOUT ANY WARRANTY; without even the implied warranty of 10 ## Octave is distributed in the hope that it will be useful, but
11 %# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 %# GNU General Public License for more details. 12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 %# 13 ## General Public License for more details.
14 %# You should have received a copy of the GNU General Public License 14 ##
15 %# along with this program; If not, see <http://www.gnu.org/licenses/>. 15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, see
17 ## <http://www.gnu.org/licenses/>.
16 18
17 %# -*- texinfo -*- 19 ## -*- texinfo -*-
18 %# @deftypefn {Function File} {[@var{sol}] =} odepkg_event_handle (@var{@@fun}, @var{time}, @var{y}, @var{flag}, [@var{par1}, @var{par2}, @dots{}]) 20 ## @deftypefn {Function File} {[@var{sol}] =} odepkg_event_handle (@var{@@fun}, @var{time}, @var{y}, @var{flag}, [@var{par1}, @var{par2}, @dots{}])
19 %# 21 ##
20 %# Return the solution of the event function that is specified as the first input argument @var{@@fun} in form of a function handle. The second input argument @var{time} 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 22 ## Return the solution of the event function that is specified as the first input argument @var{@@fun} in form of a function handle. The second input argument @var{time} 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
21 %# @table @option 23 ## @table @option
22 %# @item @code{"init"} 24 ## @item @code{"init"}
23 %# then initialize internal persistent variables of the function @command{odepkg_event_handle} and return an empty cell array of size 4, 25 ## then initialize internal persistent variables of the function @command{odepkg_event_handle} and return an empty cell array of size 4,
24 %# @item @code{"calc"} 26 ## @item @code{"calc"}
25 %# then do the evaluation of the event function and return the solution @var{sol} as type cell array of size 4, 27 ## then do the evaluation of the event function and return the solution @var{sol} as type cell array of size 4,
26 %# @item @code{"done"} 28 ## @item @code{"done"}
27 %# then cleanup internal variables of the function @command{odepkg_event_handle} and return an empty cell array of size 4. 29 ## then cleanup internal variables of the function @command{odepkg_event_handle} and return an empty cell array of size 4.
28 %# @end table 30 ## @end table
29 %# Optionally if further input arguments @var{par1}, @var{par2}, @dots{} of any type are given then pass these parameters through @command{odepkg_event_handle} to the event function. 31 ## Optionally if further input arguments @var{par1}, @var{par2}, @dots{} of any type are given then pass these parameters through @command{odepkg_event_handle} to the event function.
30 %# 32 ##
31 %# This function is an OdePkg 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. 33 ## This function is an OdePkg 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.
32 %# @end deftypefn 34 ## @end deftypefn
33 %# 35 ##
34 %# @seealso{odepkg} 36 ## @seealso{odepkg}
35 37
36 function [vretval] = odepkg_event_handle (vevefun, vt, vy, vflag, varargin) 38 function [vretval] = odepkg_event_handle (vevefun, vt, vy, vflag, varargin)
37 39
38 %# No error handling has been implemented in this function to achieve 40 ## No error handling has been implemented in this function to achieve
39 %# the highest performance available. 41 ## the highest performance available.
40 42
41 %# vretval{1} is true or false; either to terminate or to continue 43 ## vretval{1} is true or false; either to terminate or to continue
42 %# vretval{2} is the index information for which event occured 44 ## vretval{2} is the index information for which event occured
43 %# vretval{3} is the time information column vector 45 ## vretval{3} is the time information column vector
44 %# vretval{4} is the line by line result information matrix 46 ## vretval{4} is the line by line result information matrix
45 47
46 %# These persistent variables are needed to store the results and the 48 ## These persistent variables are needed to store the results and the
47 %# time value from the processing in the time stamp before, veveold 49 ## time value from the processing in the time stamp before, veveold
48 %# are the results from the event function, vtold the time stamp, 50 ## are the results from the event function, vtold the time stamp,
49 %# vretcell the return values cell array, vyold the result of the ode 51 ## vretcell the return values cell array, vyold the result of the ode
50 %# and vevecnt the counter for how often this event handling 52 ## and vevecnt the counter for how often this event handling
51 %# has been called 53 ## has been called
52 persistent veveold; persistent vtold; 54 persistent veveold; persistent vtold;
53 persistent vretcell; persistent vyold; 55 persistent vretcell; persistent vyold;
54 persistent vevecnt; 56 persistent vevecnt;
55 57
56 %# Call the event function if an event function has been defined to 58 ## Call the event function if an event function has been defined to
57 %# initialize the internal variables of the event function an to get 59 ## initialize the internal variables of the event function an to get
58 %# a value for veveold 60 ## a value for veveold
59 if (strcmp (vflag, 'init')) 61 if (strcmp (vflag, 'init'))
60 62
61 if (~iscell (vy)) 63 if (~iscell (vy))
62 vinpargs = {vevefun, vt, vy}; 64 vinpargs = {vevefun, vt, vy};
63 else 65 else
64 vinpargs = {vevefun, vt, vy{1}, vy{2}}; 66 vinpargs = {vevefun, vt, vy{1}, vy{2}};
65 vy = vy{1}; %# Delete cell element 2 67 vy = vy{1}; ## Delete cell element 2
66 end 68 end
67 if (nargin > 4) 69 if (nargin > 4)
68 vinpargs = {vinpargs{:}, varargin{:}}; 70 vinpargs = {vinpargs{:}, varargin{:}};
69 end 71 end
70 [veveold, vterm, vdir] = feval (vinpargs{:}); 72 [veveold, vterm, vdir] = feval (vinpargs{:});
71 73
72 %# We assume that all return values must be column vectors 74 ## We assume that all return values must be column vectors
73 veveold = veveold(:)'; vterm = vterm(:)'; vdir = vdir(:)'; 75 veveold = veveold(:)'; vterm = vterm(:)'; vdir = vdir(:)';
74 vtold = vt; vyold = vy; vevecnt = 1; vretcell = cell (1,4); 76 vtold = vt; vyold = vy; vevecnt = 1; vretcell = cell (1,4);
75 77
76 %# Process the event, find the zero crossings either for a rising 78 ## Process the event, find the zero crossings either for a rising
77 %# or for a falling edge 79 ## or for a falling edge
78 elseif (isempty (vflag)) 80 elseif (isempty (vflag))
79 81
80 if (~iscell (vy)) 82 if (~iscell (vy))
81 vinpargs = {vevefun, vt, vy}; 83 vinpargs = {vevefun, vt, vy};
82 else 84 else
83 vinpargs = {vevefun, vt, vy{1}, vy{2}}; 85 vinpargs = {vevefun, vt, vy{1}, vy{2}};
84 vy = vy{1}; %# Delete cell element 2 86 vy = vy{1}; ## Delete cell element 2
85 end 87 end
86 if (nargin > 4) 88 if (nargin > 4)
87 vinpargs = {vinpargs{:}, varargin{:}}; 89 vinpargs = {vinpargs{:}, varargin{:}};
88 end 90 end
89 [veve, vterm, vdir] = feval (vinpargs{:}); 91 [veve, vterm, vdir] = feval (vinpargs{:});
90 92
91 %# We assume that all return values must be column vectors 93 ## We assume that all return values must be column vectors
92 veve = veve(:)'; vterm = vterm(:)'; vdir = vdir(:)'; 94 veve = veve(:)'; vterm = vterm(:)'; vdir = vdir(:)';
93 95
94 %# Check if one or more signs of the event has changed 96 ## Check if one or more signs of the event has changed
95 vsignum = (sign (veveold) ~= sign (veve)); 97 vsignum = (sign (veveold) ~= sign (veve));
96 if (any (vsignum)) %# One or more values have changed 98 if (any (vsignum)) ## One or more values have changed
97 vindex = find (vsignum); %# Get the index of the changed values 99 vindex = find (vsignum); ## Get the index of the changed values
98 100
99 if (any (vdir(vindex) == 0)) 101 if (any (vdir(vindex) == 0))
100 %# Rising or falling (both are possible) 102 ## Rising or falling (both are possible)
101 %# Don't change anything, keep the index 103 ## Don't change anything, keep the index
102 elseif (any (vdir(vindex) == sign (veve(vindex)))) 104 elseif (any (vdir(vindex) == sign (veve(vindex))))
103 %# Detected rising or falling, need a new index 105 ## Detected rising or falling, need a new index
104 vindex = find (vdir == sign (veve)); 106 vindex = find (vdir == sign (veve));
105 else 107 else
106 %# Found a zero crossing but must not be notified 108 ## Found a zero crossing but must not be notified
107 vindex = []; 109 vindex = [];
108 end 110 end
109 111
110 %# Create new output values if a valid index has been found 112 ## Create new output values if a valid index has been found
111 if (~isempty (vindex)) 113 if (~isempty (vindex))
112 %# Change the persistent result cell array 114 ## Change the persistent result cell array
113 vretcell{1} = any (vterm(vindex)); %# Stop integration or not 115 vretcell{1} = any (vterm(vindex)); ## Stop integration or not
114 vretcell{2}(vevecnt,1) = vindex(1,1); %# Take first event found 116 vretcell{2}(vevecnt,1) = vindex(1,1); ## Take first event found
115 %# Calculate the time stamp when the event function returned 0 and 117 ## Calculate the time stamp when the event function returned 0 and
116 %# calculate new values for the integration results, we do both by 118 ## calculate new values for the integration results, we do both by
117 %# a linear interpolation 119 ## a linear interpolation
118 vtnew = vt - veve(1,vindex) * (vt - vtold) / (veve(1,vindex) - veveold(1,vindex)); 120 vtnew = vt - veve(1,vindex) * (vt - vtold) / (veve(1,vindex) - veveold(1,vindex));
119 vynew = (vy - (vt - vtnew) * (vy - vyold) / (vt - vtold))'; 121 vynew = (vy - (vt - vtnew) * (vy - vyold) / (vt - vtold))';
120 vretcell{3}(vevecnt,1) = vtnew; 122 vretcell{3}(vevecnt,1) = vtnew;
121 vretcell{4}(vevecnt,:) = vynew; 123 vretcell{4}(vevecnt,:) = vynew;
122 vevecnt = vevecnt + 1; 124 vevecnt = vevecnt + 1;
123 end %# if (~isempty (vindex)) 125 end ## if (~isempty (vindex))
124 126
125 end %# Check for one or more signs ... 127 end ## Check for one or more signs ...
126 veveold = veve; vtold = vt; vretval = vretcell; vyold = vy; 128 veveold = veve; vtold = vt; vretval = vretcell; vyold = vy;
127 129
128 elseif (strcmp (vflag, 'done')) %# Clear this event handling function 130 elseif (strcmp (vflag, 'done')) ## Clear this event handling function
129 clear ('veveold', 'vtold', 'vretcell', 'vyold', 'vevecnt'); 131 clear ('veveold', 'vtold', 'vretcell', 'vyold', 'vevecnt');
130 vretcell = cell (1,4); 132 vretcell = cell (1,4);
131 133
132 end 134 end
133
134 %# Local Variables: ***
135 %# mode: octave ***
136 %# End: ***