octave_packages/odepkg-0.8.2/odepkg_event_handle.m

   1 %# Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net>
   2 %# OdePkg - A package for solving ordinary differential equations and more
   3 %#
   4 %# This program is free software; you can redistribute it and/or modify
   5 %# it under the terms of the GNU General Public License as published by
   6 %# the Free Software Foundation; either version 2 of the License, or
   7 %# (at your option) any later version.
   8 %#
   9 %# This program is distributed in the hope that it will be useful,
  10 %# but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 %# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 %# GNU General Public License for more details.
  13 %#
  14 %# You should have received a copy of the GNU General Public License
  15 %# along with this program; If not, see <http://www.gnu.org/licenses/>.
  16
  17 %# -*- texinfo -*-
  18 %# @deftypefn {Function File} {[@var{sol}] =} odepkg_event_handle (@var{@@fun}, @var{time}, @var{y}, @var{flag}, [@var{par1}, @var{par2}, @dots{}])
  19 %#
  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
  21 %# @table @option
  22 %# @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,
  24 %# @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,
  26 %# @item  @code{"done"}
  27 %# then cleanup internal variables of the function @command{odepkg_event_handle} and return an empty cell array of size 4.
  28 %# @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.
  30 %#
  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.
  32 %# @end deftypefn
  33 %#
  34 %# @seealso{odepkg}
  35
  36 function [vretval] = odepkg_event_handle (vevefun, vt, vy, vflag, varargin)
  37
  38   %# No error handling has been implemented in this function to achieve
  39   %# the highest performance available.
  40
  41   %# vretval{1} is true or false; either to terminate or to continue
  42   %# vretval{2} is the index information for which event occured
  43   %# vretval{3} is the time information column vector
  44   %# vretval{4} is the line by line result information matrix
  45
  46   %# These persistent variables are needed to store the results and the
  47   %# time value from the processing in the time stamp before, veveold
  48   %# are the results from the event function, vtold the time stamp,
  49   %# vretcell the return values cell array, vyold the result of the ode
  50   %# and vevecnt the counter for how often this event handling
  51   %# has been called
  52   persistent veveold;  persistent vtold;
  53   persistent vretcell; persistent vyold;
  54   persistent vevecnt;
  55
  56   %# Call the event function if an event function has been defined to
  57   %# initialize the internal variables of the event function an to get
  58   %# a value for veveold
  59   if (strcmp (vflag, 'init'))
  60
  61     if (~iscell (vy))
  62       vinpargs = {vevefun, vt, vy};
  63     else
  64       vinpargs = {vevefun, vt, vy{1}, vy{2}};
  65       vy = vy{1}; %# Delete cell element 2
  66     end
  67     if (nargin > 4)
  68       vinpargs = {vinpargs{:}, varargin{:}};
  69     end
  70     [veveold, vterm, vdir] = feval (vinpargs{:});
  71
  72     %# We assume that all return values must be column vectors
  73     veveold = veveold(:)'; vterm = vterm(:)'; vdir = vdir(:)';
  74     vtold = vt; vyold = vy; vevecnt = 1; vretcell = cell (1,4);
  75
  76   %# Process the event, find the zero crossings either for a rising
  77   %# or for a falling edge
  78   elseif (isempty (vflag))
  79
  80     if (~iscell (vy))
  81       vinpargs = {vevefun, vt, vy};
  82     else
  83       vinpargs = {vevefun, vt, vy{1}, vy{2}};
  84       vy = vy{1}; %# Delete cell element 2
  85     end
  86     if (nargin > 4)
  87       vinpargs = {vinpargs{:}, varargin{:}};
  88     end
  89     [veve, vterm, vdir] = feval (vinpargs{:});
  90
  91     %# We assume that all return values must be column vectors
  92     veve = veve(:)'; vterm = vterm(:)'; vdir = vdir(:)';
  93
  94     %# Check if one or more signs of the event has changed
  95     vsignum = (sign (veveold) ~= sign (veve));
  96     if (any (vsignum))         %# One or more values have changed
  97       vindex = find (vsignum); %# Get the index of the changed values
  98
  99       if (any (vdir(vindex) == 0))
 100         %# Rising or falling (both are possible)
 101         %# Don't change anything, keep the index
 102       elseif (any (vdir(vindex) == sign (veve(vindex))))
 103         %# Detected rising or falling, need a new index
 104         vindex = find (vdir == sign (veve));
 105       else
 106         %# Found a zero crossing but must not be notified
 107         vindex = [];
 108       end
 109
 110       %# Create new output values if a valid index has been found
 111       if (~isempty (vindex))
 112         %# Change the persistent result cell array
 113         vretcell{1} = any (vterm(vindex));    %# Stop integration or not
 114         vretcell{2}(vevecnt,1) = vindex(1,1); %# Take first event found
 115         %# Calculate the time stamp when the event function returned 0 and
 116         %# calculate new values for the integration results, we do both by
 117         %# a linear interpolation
 118         vtnew = vt - veve(1,vindex) * (vt - vtold) / (veve(1,vindex) - veveold(1,vindex));
 119         vynew = (vy - (vt - vtnew) * (vy - vyold) / (vt - vtold))';
 120         vretcell{3}(vevecnt,1) = vtnew;
 121         vretcell{4}(vevecnt,:) = vynew;
 122         vevecnt = vevecnt + 1;
 123       end %# if (~isempty (vindex))
 124
 125     end %# Check for one or more signs ...
 126     veveold = veve; vtold = vt; vretval = vretcell; vyold = vy;
 127
 128   elseif (strcmp (vflag, 'done')) %# Clear this event handling function
 129     clear ('veveold', 'vtold', 'vretcell', 'vyold', 'vevecnt');
 130     vretcell = cell (1,4);
 131
 132   end
 133
 134 %!function [veve, vterm, vdir] = feveode (vt, vy, varargin)
 135 %!  veve  = vy(1); %# Which event component should be treaded
 136 %!  vterm =     1; %# Terminate if an event is found
 137 %!  vdir  =    -1; %# In which direction, -1 for falling
 138 %!function [veve, vterm, vdir] = feveide (vt, vy, vyd, varargin)
 139 %!  veve  = vy(1); %# Which event component should be treaded
 140 %!  vterm =     1; %# Terminate if an event is found
 141 %!  vdir  =    -1; %# In which direction, -1 for falling
 142 %!
 143 %!test %# First call to initialize the odepkg_event_handle function
 144 %!  odepkg_event_handle (@feveode, 0.0, [0 1 2 3], 'init', 123, 456);
 145 %!test %# Two calls to find the event that may occur with ODE/DAE syntax
 146 %!  A = odepkg_event_handle (@feveode, 2.0, [ 0 0 3 2], '', 123, 456);
 147 %!  B = odepkg_event_handle (@feveode, 3.0, [-1 0 3 2], '', 123, 456);
 148 %!  assert (A, {[], [], [], []});
 149 %!  assert (B{4}, [0 0 3 2]);
 150 %!test %# Last call to cleanup the odepkg_event_handle function
 151 %!  odepkg_event_handle (@feveode, 4.0, [0 1 2 3], 'done', 123, 456);
 152 %!test %# First call to initialize the odepkg_event_handle function
 153 %!  odepkg_event_handle (@feveide, 0.0, {[0 1 2 3], [0 1 2 3]}, 'init', 123, 456);
 154 %!test %# Two calls to find the event that may occur with IDE/DDE syntax
 155 %!  A = odepkg_event_handle (@feveide, 2.0, {[0 0 3 2], [0 0 3 2]}, '', 123, 456);
 156 %!  B = odepkg_event_handle (@feveide, 3.0, {[-1 0 3 2], [0 0 3 2]}, '', 123, 456);
 157 %!  assert (A, {[], [], [], []});
 158 %!  assert (B{4}, [0 0 3 2]);
 159 %!test %# Last call to cleanup the odepkg_event_handle function
 160 %!  odepkg_event_handle (@feveide, 4.0, {[0 1 2 3], [0 1 2 3]}, 'done', 123, 456);
 161
 162 %# Local Variables: ***
 163 %# mode: octave ***
 164 %# End: ***