[Octave-cvsupdate] octave-forge/main/odepkg/inst odepkg_equations_roessler.m, NONE, 1.1 odepkg_event_handle.m, NONE, 1.1

[Thomas T. <tr...@us...>]

Update of /cvsroot/octave/octave-forge/main/odepkg/inst
In directory sc8-pr-cvs3.sourceforge.net:/tmp/cvs-serv22374

Added Files:
	odepkg_equations_roessler.m odepkg_event_handle.m 
Log Message:
Further improved functionality. See main/odepkg/doc/ChangeLog for further details.


--- NEW FILE: odepkg_equations_roessler.m ---
%# Copyright (C) 2006, Thomas Treichl <tr...@us...>
%# OdePkg - Package for solving ordinary differential equations with octave
%#
%# This program is free software; you can redistribute it and/or modify
%# it under the terms of the GNU General Public License as published by
%# the Free Software Foundation; either version 2 of the License, or
%# (at your option) any later version.
%#
%# This program is distributed in the hope that it will be useful,
%# but WITHOUT ANY WARRANTY; without even the implied warranty of
%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%# GNU General Public License for more details.
%#
%# You should have received a copy of the GNU General Public License
%# along with this program; if not, write to the Free Software
%# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

%# -*- texinfo -*-
%# @deftypefn {Function} {@var{[y]} =} odepkg_equations_lorenz (@var{t}, @var{x})
%# TODO
%# @end deftypefn
%# 
%# @seealso{odepkg}

function y = odepkg_equations_roessler (t, x)
  y = [- ( x(2) + x(3) );
       x(1) + 0.2 * x(2);
       0.2 + x(1) * x(3) - 5.7 * x(3)];

%!test A = odeset ('MaxStep', 1e-1);
%!     [t, y] = ode78 (@odepkg_equations_roessler, [0 70], [0.1 0.3 0.1], A);

%!demo
%!
%! A = odeset ('MaxStep', 1e-1);
%! [t, y] = ode78 (@odepkg_equations_roessler, [0 70], [0.1 0.3 0.1], A);
%!
%! subplot (2, 2, 1); grid ('on'); plot (t, y(:,1), '-b;f_x(t);', ...
%!    t, y(:,2), '-g;f_y(t);', t, y(:,3), '-r;f_z(t);');
%! subplot (2, 2, 2); grid ('on'); plot (y(:,1), y(:,2), '-b;f_{xyz}(x, y);');
%! subplot (2, 2, 3); grid ('on'); plot (y(:,2), y(:,3), '-b;f_{xyz}(y, z);');
%! subplot (2, 2, 4); grid ('on'); plot3 (y(:,1), y(:,2), y(:,3), ...
%!    '-b;f_{xyz}(x, y, z);');
%!
%! % -------------------------------------------------------------------------
%! % TODO

%# Local Variables: ***
%# mode: octave ***
%# End: ***
--- NEW FILE: odepkg_event_handle.m ---
%# Copyright (C) 2006, Thomas Treichl <tr...@us...>
%# OdePkg - Package for solving ordinary differential equations with octave
%#
%# This program is free software; you can redistribute it and/or modify
%# it under the terms of the GNU General Public License as published by
%# the Free Software Foundation; either version 2 of the License, or
%# (at your option) any later version.
%#
%# This program is distributed in the hope that it will be useful,
%# but WITHOUT ANY WARRANTY; without even the implied warranty of
%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%# GNU General Public License for more details.
%#
%# You should have received a copy of the GNU General Public License
%# along with this program; if not, write to the Free Software
%# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

%# -*- texinfo -*-
%# @deftypefn {Function} {@var{[sol]} =} odepkg_event_handle (@var{@@fun, time, y, flag, [P1, P2, @dots{}]})
%# Evaluates the event function of an ode solver call, if before the event option was set with the command odeset. This function is an internal helper function, therefore this function should never be directly called by a user. No error handling has been implemented in this function.
%# @end deftypefn
%#
%# @seealso{odepkg}

%# Maintainer: Thomas Treichl
%# Created: 20060928
%# ChangeLog:

function [vretval] = odepkg_event_handle (vevefun, vt, vy, vflag, varargin)
  %# vretval{1} is true or false; either to terminate or to continue
  %# vretval{2} is the index information for which event occured
  %# vretval{3} is the time information column vector
  %# vretval{4} is the line by line result information matrix

  %# These persistent variables are needed to store the results and the
  %# time value from the processing in the time stamp before, veveold
  %# are the results from the event function, vtold the time stamp,
  %# vretcell the return values cell array, vyold the result of the ode
  %# and vevecnt the counter for how often this event handling
  %# has been called
  persistent veveold;  persistent vtold;
  persistent vretcell; persistent vyold;
  persistent vevecnt;

  %# Call the event function if an event function has been defined to
  %# initialize the internal variables of the event function an to get
  %# a value for veveold
  if (strcmp (vflag, 'init') == true)

    if (isempty (varargin) == true)
      [veveold, vterm, vdir] = feval (vevefun, vt, vy);
    else
      [veveold, vterm, vdir] = feval (vevefun, vt, vy, varargin);
    end

    %# My first implementation was to return row vectors from the event
    %# functions but this is faulty if we use LabMat events. So this has
    %# changed 20060928 and by now these return vectors must be column
    %# vectors. The not so nice implementation then is:
    veveold = veveold(:)'; vterm = vterm(:)'; vdir = vdir(:)';

    vtold = vt; %# veveold has been set in the code lines before
    for vcnt = 1:4, vretcell{vcnt} = []; end
    vyold = vy; vevecnt = 1;

  %# Process the event, find the zero crossings either for a rising
  %# or for a falling edge
  elseif (isempty (vflag) == true)

    %# Call the event function if an event function has been defined to
    %# to get the value(s) veve
    if (isempty (varargin) == true)
      [veve, vterm, vdir] = feval (vevefun, vt, vy);
    else
      [veve, vterm, vdir] = feval (vevefun, vt, vy, varargin);
    end
    veve = veve(:)'; vterm = vterm(:)'; vdir = vdir(:)';

    %# Check if one or more signs of the event function results changed
    vsignum = (sign (veveold) ~= sign (veve));
    if (any (vsignum) == true)     %# One or more values have changed
      vindex = find (vsignum);     %# Get the index of the changed values

      if (any (vdir(vindex) == 0)) %# Rising or falling (both are possible)
        %# Don't change anything, keep the index
      elseif (any (vdir(vindex) == sign (veve(vindex))) == true)
        %# Detected rising or falling, need a new index
        vindex = find (vdir == sign (veve));
      else
        %# Found a zero crossing that will not be notified
        vindex = [];
      end

      %# We create a new output values if we found a valid index 
      if (isempty (vindex) == false)
        %# Change the persistent result cell array
        vretcell{1} = any (vterm(vindex) == true);  %# Stop integration or not
        vretcell{2}(vevecnt,1) = vindex(1,1); %# Take first event that occurs
        %# 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
        vtnew = vt - veve(1,vindex) * (vt - vtold) / (veve(1,vindex) - veveold(1,vindex));
        vynew = (vy - (vt - vtnew) * (vy - vyold) / (vt - vtold));
        vretcell{3}(vevecnt,1) = vtnew;
        vretcell{4}(vevecnt,:) = vynew';
        vevecnt = vevecnt + 1;
      end

    end %# Check if one or more signs ...
    veveold = veve; vtold = vt;
    vretval = vretcell; vyold = vy;

  %# Reset this event handling function
  elseif (strcmp (vflag, 'done') == true)
    clear ('veveold', 'vtold', 'vretcell', 'vyold', 'vevecnt');
    for vcnt = 1:4, vretcell{vcnt} = []; end

  end %# if (strcmp (vflag, ...) == true)

%# Local Variables: ***
%# mode: octave ***
%# End: ***