Documentation of inteul

Global Index (all files) (short | long) | Local contents | Local Index (files in subdir) (short | long)

Function Synopsis

[Timesim, Result] = inteul(SIM_F, Times, XINIT, Options, Control, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10)

Help text

 Integration of ODE's using EULER's formula

 This function implements the vectorized EULER integrator.

 Syntax:  [Timevec, Result] = inteul(SIM_F, Times, XINIT, Options, Control, P1-Px)

 Input parameters:
    SIM_F     - function to simulate
    Times     - vector containing TStart (start time - optional) and TEnd (end time)
    XINIT     - initial values for the simulation of the individuals
                every row contains the initial values for 1 individual
    Options   - vector containing [ERRMAX: max. rel. error, STEPMIN: min. stepsize, STEPMAX: max. stepsize, ...
                                   (optional) NCONTR: number of control variables]
    Control   - control vector/matrix, every row is one control vector
                format:   1. ControlVariable  2. ControlVariable ...  NCONTR. ControlVariable
                         --------------------------------------------------------------------
                  Ind1: | 1:ControlSteps      1:ControlSteps     ...  1:ControlSteps
                  Ind2: | 1:ControlSteps      1:ControlSteps     ...  1:ControlSteps
    P1 - Px   - (optional) parameters passed through to simulation function

 Output parameters:
    Timesim   - vector containing time of simulation steps
    Result    - results, state vectors
                format:   1. ControlVariable  2. ControlVariable ...  NCONTR. ControlVariable
                         --------------------------------------------------------------------
                  Ind1: | 1:ControlSteps      1:ControlSteps     ...  1:ControlSteps
                  Ind2: | 1:ControlSteps      1:ControlSteps     ...  1:ControlSteps

 See also: objdopi, intrk4

Cross-Reference Information

Listing of function inteul

% Author:     Hartmut Pohlheim
% History:    24.10.94     file created
%             26.02.95     aditional parameters added
%             16.05.95     same format as rk23/rk45, here vectorized
%             17.05.95     final definition of parameterformat (quite complicated)
%             18.05.95     renamed to inteul.m 
%             29.04.96     definition of simevalstr changed, made problem with error
%             21.01.98     small adjustement for Matlab5 (TIMEVEC(:), Part(:))


function [Timesim, Result] = inteul(SIM_F, Times, XINIT, Options, Control, P1, P2, P3, P4, P5, P6, P7, P8, P9, P10)

% Time construction
   if nargin < 2, Times = []; end
   if isempty(Times), error('At least the end time must be defined for the simulation.'); end
   [ts, te] = size(Times);
   if (ts*te) == 1, Times = [ 0, Times(1)]; end
   TStart = Times(1); TEnd = Times(2);
   
% check input parameter XINIT
   if nargin < 3, XINIT = []; end
   if isempty(XINIT), error('The initial value(s) must be defined for the simulation.'); end

% check input parameter Options
   if nargin < 4, Options = []; end
   if isempty(Options), Options = 1e-3; end
   if length(Options) < 2, Options(2) = (TEnd-TStart)/2000; end
   if length(Options) < 3, Options(3) = (TEnd-TStart)/50; end
   if length(Options) < 4, Options(4) = 0; end
   ERRMAX = Options(1); STEPMIN = Options(2); STEPMAX = Options(3); NCONTR = Options(4);

% check input parameter Control
   if nargin < 5, Control = []; end
   if isempty(Control),
      if NCONTR ~= 0, error('No control vector/matrix defined for the simulation.'); end,
   end

% Compute population parameters
   [NIND, Nvar] = size(XINIT);
   [NI, NCont] = size(Control);

% wenn Control noch nicht kopiert oder konstant fuer alle Schritte
   if NCONTR > 0,
      if rem(NCont, NCONTR) ~= 0, error('NCONTR and size of Control matrix disagree.'); end
      % compute number of control steps
      ControlSteps = NCont/NCONTR;
      % if Control is konstant for all steps, we need only one Control for all steps;
      % however, 2 steps are needed for table lookup
      if ControlSteps == 1,
         Control = expandm(Control, [1, 2]);
         NCont = 2 * NCONTR; ControlSteps = 2;
      end
      % compute time vector according to number of control steps
      TIMEVEC = linspace(TStart, TEnd, ControlSteps);
      TIMEVEC = TIMEVEC(:);
      % reshape Control for table lookup, every row contains all variables for one time step
      % format:                  1. Ind      2. Ind
      %         ControlStep 1    1:NCONTR    1:NCONTR
      %         ControlStep 2    1:NCONTR    1:NCONTR
      ControlTable = [];
      for i = 1:ControlSteps, 
         % Part(:) = Control(:, (i-1)*NCONTR+1:i*NCONTR)';
         Part = Control(:, i:ControlSteps:i+(NCONTR-1)*ControlSteps)';
         Part = Part(:);
         ControlTable = [ControlTable; Part'];
      end
   end

% Build string of simulation function
   simevalstr = ['xdot = ' SIM_F]; flag = 1;
   simevalstr = [simevalstr, '(T, X, U, flag'];
   for i = 1:nargin-5
      simevalstr = [simevalstr, ',P', int2str(i)];
   end
   simevalstr = [simevalstr, ');'];

% perform simulation steps
   T = TStart; Timesim = [T]; DoSim = 1; X = XINIT'; Result = [XINIT];
   SimStep = STEPMIN;
   while DoSim == 1,
      % compute size of next step
      if (TEnd-TStart) > 0,
         if T+SimStep > TEnd, SimStep = TEnd - T; end,
      else 
         if T+SimStep < TEnd, SimStep = TEnd - T; end,
      end
      % check, if simulation reached end of time, depends on direction of simulation
      if (TEnd-TStart) > 0, if T+eps >= TEnd, DoSim = 0; end,
      else                  if T-eps <= TEnd, DoSim = 0; end, end

      if DoSim == 1,
         if NCONTR ~= 0,
            % lookup the contol values for the actual step
            U = table1([TIMEVEC, ControlTable], T+SimStep/2); 
            U = reshape(U, NCONTR, NIND);      % one column per individium
         else
            U = [];
         end

         % compute 1 simulation step
         eval(simevalstr);
         X = X + (xdot * SimStep);
         T = T + SimStep;

         % X(:,6) = X(:,6) .* (X(:,6) > 0.0 & X(:,6) < 25.0) + 25.0 .* (X(:,6) >= 25.0);

         % save Result in 1 row per individual,
         % all individuals of 1 time step at once in rows underneath,
         % then the individuls of the next time step 
         Result = [Result; X'];
         Timesim = [Timesim; T];

      end
   end


% End of function