octave_packages/odepkg-0.8.2/ode78d.m

   1 %# Copyright (C) 2008-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{}] =} ode78d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
  19 %# @deftypefnx {Command} {[@var{sol}] =} ode78d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
  20 %# @deftypefnx {Command} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =} ode78d (@var{@@fun}, @var{slot}, @var{init}, @var{lags}, @var{hist}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
  21 %#
  22 %# This function file can be used to solve a set of non--stiff delay differential equations (non--stiff DDEs) with a modified version of the well known explicit Runge--Kutta method of order (7,8).
  23 %#
  24 %# If this function is called with no return argument then plot the solution over time in a figure window while solving the set of DDEs that are defined in a function and specified by the function handle @var{@@fun}. The second input argument @var{slot} is a double vector that defines the time slot, @var{init} is a double vector that defines the initial values of the states, @var{lags} is a double vector that describes the lags of time, @var{hist} is a double matrix and describes the history of the DDEs, @var{opt} can optionally be a structure array that keeps the options created with the command @command{odeset} and @var{par1}, @var{par2}, @dots{} can optionally be other input arguments of any type that have to be passed to the function defined by @var{@@fun}.
  25 %#
  26 %# In other words, this function will solve a problem of the form
  27 %# @example
  28 %# dy/dt = fun (t, y(t), y(t-lags(1), y(t-lags(2), @dots{})))
  29 %# y(slot(1)) = init
  30 %# y(slot(1)-lags(1)) = hist(1), y(slot(1)-lags(2)) = hist(2), @dots{}
  31 %# @end example
  32 %#
  33 %# If this function is called with one return argument then return the solution @var{sol} of type structure array after solving the set of DDEs. The solution @var{sol} has the fields @var{x} of type double column vector for the steps chosen by the solver, @var{y} of type double column vector for the solutions at each time step of @var{x}, @var{solver} of type string for the solver name and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector that keep the informations of the event function if an event function handle is set in the option argument @var{opt}.
  34 %#
  35 %# If this function is called with more than one return argument then return the time stamps @var{t}, the solution values @var{y} and optionally the extended time stamp information @var{xe}, the extended solution information @var{ye} and the extended index information @var{ie} all of type double column vector.
  36 %#
  37 %# For example:
  38 %# @itemize @minus
  39 %# @item
  40 %# the following code solves an anonymous implementation of a chaotic behavior
  41 %#
  42 %# @example
  43 %# fcao = @@(vt, vy, vz) [2 * vz / (1 + vz^9.65) - vy];
  44 %#
  45 %# vopt = odeset ("NormControl", "on", "RelTol", 1e-3);
  46 %# vsol = ode78d (fcao, [0, 100], 0.5, 2, 0.5, vopt);
  47 %#
  48 %# vlag = interp1 (vsol.x, vsol.y, vsol.x - 2);
  49 %# plot (vsol.y, vlag); legend ("fcao (t,y,z)");
  50 %# @end example
  51 %#
  52 %# @item
  53 %# to solve the following problem with two delayed state variables
  54 %#
  55 %# @example
  56 %# d y1(t)/dt = -y1(t)
  57 %# d y2(t)/dt = -y2(t) + y1(t-5)
  58 %# d y3(t)/dt = -y3(t) + y2(t-10)*y1(t-10)
  59 %# @end example
  60 %#
  61 %# one might do the following
  62 %#
  63 %# @example
  64 %# function f = fun (t, y, yd)
  65 %# f(1) = -y(1);                   %% y1' = -y1(t)
  66 %# f(2) = -y(2) + yd(1,1);         %% y2' = -y2(t) + y1(t-lags(1))
  67 %# f(3) = -y(3) + yd(2,2)*yd(1,2); %% y3' = -y3(t) + y2(t-lags(2))*y1(t-lags(2))
  68 %# endfunction
  69 %# T = [0,20]
  70 %# res = ode78d (@@fun, T, [1;1;1], [5, 10], ones (3,2));
  71 %# @end example
  72 %#
  73 %# @end itemize
  74 %# @end deftypefn
  75 %#
  76 %# @seealso{odepkg}
  77
  78 function [varargout] = ode78d (vfun, vslot, vinit, vlags, vhist, varargin)
  79
  80   if (nargin == 0) %# Check number and types of all input arguments
  81     help ('ode78d');
  82     error ('OdePkg:InvalidArgument', ...
  83       'Number of input arguments must be greater than zero');
  84
  85   elseif (nargin < 5)
  86     print_usage;
  87
  88   elseif (~isa (vfun, 'function_handle'))
  89     error ('OdePkg:InvalidArgument', ...
  90       'First input argument must be a valid function handle');
  91
  92   elseif (~isvector (vslot) || length (vslot) < 2)
  93     error ('OdePkg:InvalidArgument', ...
  94       'Second input argument must be a valid vector');
  95
  96   elseif (~isvector (vinit) || ~isnumeric (vinit))
  97     error ('OdePkg:InvalidArgument', ...
  98       'Third input argument must be a valid numerical value');
  99
 100   elseif (~isvector (vlags) || ~isnumeric (vlags))
 101     error ('OdePkg:InvalidArgument', ...
 102       'Fourth input argument must be a valid numerical value');
 103
 104   elseif ~(isnumeric (vhist) || isa (vhist, 'function_handle'))
 105     error ('OdePkg:InvalidArgument', ...
 106       'Fifth input argument must either be numeric or a function handle');
 107
 108   elseif (nargin >= 6)
 109
 110     if (~isstruct (varargin{1}))
 111       %# varargin{1:len} are parameters for vfun
 112       vodeoptions = odeset;
 113       vfunarguments = varargin;
 114
 115     elseif (length (varargin) > 1)
 116       %# varargin{1} is an OdePkg options structure vopt
 117       vodeoptions = odepkg_structure_check (varargin{1}, 'ode78d');
 118       vfunarguments = {varargin{2:length(varargin)}};
 119
 120     else %# if (isstruct (varargin{1}))
 121       vodeoptions = odepkg_structure_check (varargin{1}, 'ode78d');
 122       vfunarguments = {};
 123
 124     end
 125
 126   else %# if (nargin == 5)
 127     vodeoptions = odeset;
 128     vfunarguments = {};
 129   end
 130
 131   %# Start preprocessing, have a look which options have been set in
 132   %# vodeoptions. Check if an invalid or unused option has been set and
 133   %# print warnings.
 134   vslot = vslot(:)'; %# Create a row vector
 135   vinit = vinit(:)'; %# Create a row vector
 136   vlags = vlags(:)'; %# Create a row vector
 137
 138   %# Check if the user has given fixed points of time
 139   if (length (vslot) > 2), vstepsizegiven = true; %# Step size checking
 140   else vstepsizegiven = false; end
 141
 142   %# Get the default options that can be set with 'odeset' temporarily
 143   vodetemp = odeset;
 144
 145   %# Implementation of the option RelTol has been finished. This option
 146   %# can be set by the user to another value than default value.
 147   if (isempty (vodeoptions.RelTol) && ~vstepsizegiven)
 148     vodeoptions.RelTol = 1e-6;
 149     warning ('OdePkg:InvalidOption', ...
 150       'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
 151   elseif (~isempty (vodeoptions.RelTol) && vstepsizegiven)
 152     warning ('OdePkg:InvalidOption', ...
 153       'Option "RelTol" will be ignored if fixed time stamps are given');
 154   %# This implementation has been added to odepkg_structure_check.m
 155   %# elseif (~isscalar (vodeoptions.RelTol) && ~vstepsizegiven)
 156   %# error ('OdePkg:InvalidOption', ...
 157   %#   'Option "RelTol" must be set to a scalar value for this solver');
 158   end
 159
 160   %# Implementation of the option AbsTol has been finished. This option
 161   %# can be set by the user to another value than default value.
 162   if (isempty (vodeoptions.AbsTol) && ~vstepsizegiven)
 163     vodeoptions.AbsTol = 1e-6;
 164     warning ('OdePkg:InvalidOption', ...
 165       'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
 166   elseif (~isempty (vodeoptions.AbsTol) && vstepsizegiven)
 167     warning ('OdePkg:InvalidOption', ...
 168       'Option "AbsTol" will be ignored if fixed time stamps are given');
 169   else %# create column vector
 170     vodeoptions.AbsTol = vodeoptions.AbsTol(:);
 171   end
 172
 173   %# Implementation of the option NormControl has been finished. This
 174   %# option can be set by the user to another value than default value.
 175   if (strcmp (vodeoptions.NormControl, 'on')), vnormcontrol = true;
 176   else vnormcontrol = false;
 177   end
 178
 179   %# Implementation of the option NonNegative has been finished. This
 180   %# option can be set by the user to another value than default value.
 181   if (~isempty (vodeoptions.NonNegative))
 182     if (isempty (vodeoptions.Mass)), vhavenonnegative = true;
 183     else
 184       vhavenonnegative = false;
 185       warning ('OdePkg:InvalidOption', ...
 186         'Option "NonNegative" will be ignored if mass matrix is set');
 187     end
 188   else vhavenonnegative = false;
 189   end
 190
 191   %# Implementation of the option OutputFcn has been finished. This
 192   %# option can be set by the user to another value than default value.
 193   if (isempty (vodeoptions.OutputFcn) && nargout == 0)
 194     vodeoptions.OutputFcn = @odeplot;
 195     vhaveoutputfunction = true;
 196   elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
 197   else vhaveoutputfunction = true;
 198   end
 199
 200   %# Implementation of the option OutputSel has been finished. This
 201   %# option can be set by the user to another value than default value.
 202   if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
 203   else vhaveoutputselection = false; end
 204
 205   %# Implementation of the option Refine has been finished. This option
 206   %# can be set by the user to another value than default value.
 207   if (isequal (vodeoptions.Refine, vodetemp.Refine)), vhaverefine = true;
 208   else vhaverefine = false; end
 209
 210   %# Implementation of the option Stats has been finished. This option
 211   %# can be set by the user to another value than default value.
 212
 213   %# Implementation of the option InitialStep has been finished. This
 214   %# option can be set by the user to another value than default value.
 215   if (isempty (vodeoptions.InitialStep) && ~vstepsizegiven)
 216     vodeoptions.InitialStep = abs (vslot(1,1) - vslot(1,2)) / 10;
 217     vodeoptions.InitialStep = vodeoptions.InitialStep / 10^vodeoptions.Refine;
 218     warning ('OdePkg:InvalidOption', ...
 219       'Option "InitialStep" not set, new value %f is used', vodeoptions.InitialStep);
 220   end
 221
 222   %# Implementation of the option MaxStep has been finished. This option
 223   %# can be set by the user to another value than default value.
 224   if (isempty (vodeoptions.MaxStep) && ~vstepsizegiven)
 225     vodeoptions.MaxStep = abs (vslot(1,1) - vslot(1,length (vslot))) / 10;
 226     %# vodeoptions.MaxStep = vodeoptions.MaxStep / 10^vodeoptions.Refine;
 227     warning ('OdePkg:InvalidOption', ...
 228       'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep);
 229   end
 230
 231   %# Implementation of the option Events has been finished. This option
 232   %# can be set by the user to another value than default value.
 233   if (~isempty (vodeoptions.Events)), vhaveeventfunction = true;
 234   else vhaveeventfunction = false; end
 235
 236   %# The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored
 237   %# by this solver because this solver uses an explicit Runge-Kutta
 238   %# method and therefore no Jacobian calculation is necessary
 239   if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian))
 240     warning ('OdePkg:InvalidOption', ...
 241       'Option "Jacobian" will be ignored by this solver');
 242   end
 243   if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
 244     warning ('OdePkg:InvalidOption', ...
 245       'Option "JPattern" will be ignored by this solver');
 246   end
 247   if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
 248     warning ('OdePkg:InvalidOption', ...
 249       'Option "Vectorized" will be ignored by this solver');
 250   end
 251   if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol))
 252     warning ('OdePkg:InvalidArgument', ...
 253       'Option "NewtonTol" will be ignored by this solver');
 254   end
 255   if (~isequal (vodeoptions.MaxNewtonIterations,...
 256                 vodetemp.MaxNewtonIterations))
 257     warning ('OdePkg:InvalidArgument', ...
 258       'Option "MaxNewtonIterations" will be ignored by this solver');
 259   end
 260
 261   %# Implementation of the option Mass has been finished. This option
 262   %# can be set by the user to another value than default value.
 263   if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
 264     vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
 265   elseif (isa (vodeoptions.Mass, 'function_handle'))
 266     vhavemasshandle = true; %# mass defined by a function handle
 267   else %# no mass matrix - creating a diag-matrix of ones for mass
 268     vhavemasshandle = false; %# vmass = diag (ones (length (vinit), 1), 0);
 269   end
 270
 271   %# Implementation of the option MStateDependence has been finished.
 272   %# This option can be set by the user to another value than default
 273   %# value.
 274   if (strcmp (vodeoptions.MStateDependence, 'none'))
 275     vmassdependence = false;
 276   else vmassdependence = true;
 277   end
 278
 279   %# Other options that are not used by this solver. Print a warning
 280   %# message to tell the user that the option(s) is/are ignored.
 281   if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern))
 282     warning ('OdePkg:InvalidOption', ...
 283       'Option "MvPattern" will be ignored by this solver');
 284   end
 285   if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular))
 286     warning ('OdePkg:InvalidOption', ...
 287       'Option "MassSingular" will be ignored by this solver');
 288   end
 289   if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope))
 290     warning ('OdePkg:InvalidOption', ...
 291       'Option "InitialSlope" will be ignored by this solver');
 292   end
 293   if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder))
 294     warning ('OdePkg:InvalidOption', ...
 295       'Option "MaxOrder" will be ignored by this solver');
 296   end
 297   if (~isequal (vodeoptions.BDF, vodetemp.BDF))
 298     warning ('OdePkg:InvalidOption', ...
 299       'Option "BDF" will be ignored by this solver');
 300   end
 301
 302   %# Starting the initialisation of the core solver ode78d
 303   vtimestamp  = vslot(1,1);           %# timestamp = start time
 304   vtimelength = length (vslot);       %# length needed if fixed steps
 305   vtimestop   = vslot(1,vtimelength); %# stop time = last value
 306
 307   if (~vstepsizegiven)
 308     vstepsize = vodeoptions.InitialStep;
 309     vminstepsize = (vtimestop - vtimestamp) / (1/eps);
 310   else %# If step size is given then use the fixed time steps
 311     vstepsize = abs (vslot(1,1) - vslot(1,2));
 312     vminstepsize = eps; %# vslot(1,2) - vslot(1,1) - eps;
 313   end
 314
 315   vretvaltime = vtimestamp; %# first timestamp output
 316   if (vhaveoutputselection) %# first solution output
 317     vretvalresult = vinit(vodeoptions.OutputSel);
 318   else vretvalresult = vinit;
 319   end
 320
 321   %# Initialize the OutputFcn
 322   if (vhaveoutputfunction)
 323     feval (vodeoptions.OutputFcn, vslot', ...
 324       vretvalresult', 'init', vfunarguments{:});
 325   end
 326
 327   %# Initialize the History
 328   if (isnumeric (vhist))
 329     vhmat = vhist;
 330     vhavehistnumeric = true;
 331   else %# it must be a function handle
 332     for vcnt = 1:length (vlags);
 333       vhmat(:,vcnt) = feval (vhist, (vslot(1)-vlags(vcnt)), vfunarguments{:});
 334     end
 335     vhavehistnumeric = false;
 336   end
 337
 338   %# Initialize DDE variables for history calculation
 339   vsaveddetime = [vtimestamp - vlags, vtimestamp]';
 340   vsaveddeinput = [vhmat, vinit']';
 341   vsavedderesult = [vhmat, vinit']';
 342
 343   %# Initialize the EventFcn
 344   if (vhaveeventfunction)
 345     odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
 346       {vretvalresult', vhmat}, 'init', vfunarguments{:});
 347   end
 348
 349   vpow = 1/8;                                  %# MC2001: see p.91 in Ascher & Petzold
 350   va = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;    %# The 7(8) coefficients
 351         1/18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0; %# Coefficients proved, tt 20060827
 352         1/48, 1/16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 353         1/32, 0, 3/32, 0, 0, 0, 0, 0, 0, 0, 0, 0;
 354         5/16, 0, -75/64, 75/64, 0, 0, 0, 0, 0, 0, 0, 0;
 355         3/80, 0, 0, 3/16, 3/20, 0, 0, 0, 0, 0, 0, 0;
 356         29443841/614563906, 0, 0, 77736538/692538347, -28693883/1125000000, ...
 357             23124283/1800000000, 0, 0, 0, 0, 0, 0;
 358         16016141/946692911, 0, 0, 61564180/158732637, 22789713/633445777, ...
 359             545815736/2771057229, -180193667/1043307555, 0, 0, 0, 0, 0;
 360         39632708/573591083, 0, 0, -433636366/683701615, -421739975/2616292301, ...
 361             100302831/723423059, 790204164/839813087, 800635310/3783071287, 0, 0, 0, 0;
 362         246121993/1340847787, 0, 0, -37695042795/15268766246, -309121744/1061227803, ...
 363             -12992083/490766935, 6005943493/2108947869, 393006217/1396673457, ...
 364             123872331/1001029789, 0, 0, 0;
 365         -1028468189/846180014, 0, 0, 8478235783/508512852, 1311729495/1432422823, ...
 366             -10304129995/1701304382, -48777925059/3047939560, 15336726248/1032824649, ...
 367             -45442868181/3398467696, 3065993473/597172653, 0, 0;
 368         185892177/718116043, 0, 0, -3185094517/667107341, -477755414/1098053517, ...
 369             -703635378/230739211, 5731566787/1027545527, 5232866602/850066563, ...
 370             -4093664535/808688257, 3962137247/1805957418, 65686358/487910083, 0;
 371         403863854/491063109, 0, 0, -5068492393/434740067, -411421997/543043805, ...
 372             652783627/914296604, 11173962825/925320556, -13158990841/6184727034, ...
 373             3936647629/1978049680, -160528059/685178525, 248638103/1413531060, 0];
 374   vb7 = [13451932/455176623; 0; 0; 0; 0; -808719846/976000145; ...
 375             1757004468/5645159321; 656045339/265891186; -3867574721/1518517206; ...
 376             465885868/322736535; 53011238/667516719; 2/45; 0];
 377   vb8 = [14005451/335480064; 0; 0; 0; 0; -59238493/1068277825; 181606767/758867731; ...
 378             561292985/797845732; -1041891430/1371343529; 760417239/1151165299; ...
 379             118820643/751138087; -528747749/2220607170; 1/4];
 380   vc = sum (va, 2);
 381
 382   %# The solver main loop - stop if the endpoint has been reached
 383   vcntloop = 2; vcntcycles = 1; vu = vinit; vk = vu' * zeros(1,13);
 384   vcntiter = 0; vunhandledtermination = true;
 385   while ((vtimestamp < vtimestop && vstepsize >= vminstepsize))
 386
 387     %# Hit the endpoint of the time slot exactely
 388     if ((vtimestamp + vstepsize) > vtimestop)
 389       vstepsize = vtimestop - vtimestamp; end
 390
 391     %# Estimate the thirteen results when using this solver
 392     for j = 1:13
 393       vthetime  = vtimestamp + vc(j,1) * vstepsize;
 394       vtheinput = vu' + vstepsize * vk(:,1:j-1) * va(j,1:j-1)';
 395       %# Claculate the history values (or get them from an external
 396       %# function) that are needed for the next step of solving
 397       if (vhavehistnumeric)
 398         for vcnt = 1:length (vlags)
 399           %# Direct implementation of a 'quadrature cubic Hermite interpolation'
 400           %# found at the Faculty for Mathematics of the University of Stuttgart
 401           %# http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage1269
 402           vnumb = find (vthetime - vlags(vcnt) >= vsaveddetime);
 403           velem = min (vnumb(end), length (vsaveddetime) - 1);
 404           vstep = vsaveddetime(velem+1) - vsaveddetime(velem);
 405           vdiff = (vthetime - vlags(vcnt) - vsaveddetime(velem)) / vstep;
 406           vsubs = 1 - vdiff;
 407           %# Calculation of the coefficients for the interpolation algorithm
 408           vua = (1 + 2 * vdiff) * vsubs^2;
 409           vub = (3 - 2 * vdiff) * vdiff^2;
 410           vva = vstep * vdiff * vsubs^2;
 411           vvb = -vstep * vsubs * vdiff^2;
 412           vhmat(:,vcnt) = vua * vsaveddeinput(velem,:)' + ...
 413               vub * vsaveddeinput(velem+1,:)' + ...
 414               vva * vsavedderesult(velem,:)' + ...
 415               vvb * vsavedderesult(velem+1,:)';
 416         end
 417       else %# the history must be a function handle
 418         for vcnt = 1:length (vlags)
 419           vhmat(:,vcnt) = feval ...
 420             (vhist, vthetime - vlags(vcnt), vfunarguments{:});
 421         end
 422       end
 423
 424       if (vhavemasshandle)   %# Handle only the dynamic mass matrix,
 425         if (vmassdependence) %# constant mass matrices have already
 426           vmass = feval ...  %# been set before (if any)
 427             (vodeoptions.Mass, vthetime, vtheinput, vfunarguments{:});
 428         else                 %# if (vmassdependence == false)
 429           vmass = feval ...  %# then we only have the time argument
 430             (vodeoptions.Mass, vthetime, vfunarguments{:});
 431         end
 432         vk(:,j) = vmass \ feval ...
 433           (vfun, vthetime, vtheinput, vhmat, vfunarguments{:});
 434       else
 435         vk(:,j) = feval ...
 436           (vfun, vthetime, vtheinput, vhmat, vfunarguments{:});
 437       end
 438     end
 439
 440     %# Compute the 7th and the 8th order estimation
 441     y7 = vu' + vstepsize * (vk * vb7);
 442     y8 = vu' + vstepsize * (vk * vb8);
 443     if (vhavenonnegative)
 444       vu(vodeoptions.NonNegative) = abs (vu(vodeoptions.NonNegative));
 445       y7(vodeoptions.NonNegative) = abs (y7(vodeoptions.NonNegative));
 446       y8(vodeoptions.NonNegative) = abs (y8(vodeoptions.NonNegative));
 447     end
 448     vSaveVUForRefine = vu;
 449
 450     %# Calculate the absolute local truncation error and the acceptable error
 451     if (~vstepsizegiven)
 452       if (~vnormcontrol)
 453         vdelta = y8 - y7;
 454         vtau = max (vodeoptions.RelTol * vu', vodeoptions.AbsTol);
 455       else
 456         vdelta = norm (y8 - y7, Inf);
 457         vtau = max (vodeoptions.RelTol * max (norm (vu', Inf), 1.0), ...
 458                     vodeoptions.AbsTol);
 459       end
 460     else %# if (vstepsizegiven == true)
 461       vdelta = 1; vtau = 2;
 462     end
 463
 464     %# If the error is acceptable then update the vretval variables
 465     if (all (vdelta <= vtau))
 466       vtimestamp = vtimestamp + vstepsize;
 467       vu = y8'; %# MC2001: the higher order estimation as "local extrapolation"
 468       vretvaltime(vcntloop,:) = vtimestamp;
 469       if (vhaveoutputselection)
 470         vretvalresult(vcntloop,:) = vu(vodeoptions.OutputSel);
 471       else
 472         vretvalresult(vcntloop,:) = vu;
 473       end
 474       vcntloop = vcntloop + 1; vcntiter = 0;
 475
 476       %# Update DDE values for next history calculation
 477       vsaveddetime(end+1) = vtimestamp;
 478       vsaveddeinput(end+1,:) = vtheinput';
 479       vsavedderesult(end+1,:) = vu;
 480
 481       %# Call plot only if a valid result has been found, therefore this
 482       %# code fragment has moved here. Stop integration if plot function
 483       %# returns false
 484       if (vhaveoutputfunction)
 485         if (vhaverefine)                  %# Do interpolation
 486           for vcnt = 0:vodeoptions.Refine %# Approximation between told and t
 487             vapproxtime = (vcnt + 1) * vstepsize / (vodeoptions.Refine + 2);
 488             vapproxvals = vSaveVUForRefine' + vapproxtime * (vk * vb8);
 489             if (vhaveoutputselection)
 490               vapproxvals = vapproxvals(vodeoptions.OutputSel);
 491             end
 492             feval (vodeoptions.OutputFcn, (vtimestamp - vstepsize) + vapproxtime, ...
 493               vapproxvals, [], vfunarguments{:});
 494           end
 495         end
 496         vpltret = feval (vodeoptions.OutputFcn, vtimestamp, ...
 497           vretvalresult(vcntloop-1,:)', [], vfunarguments{:});
 498         if (vpltret), vunhandledtermination = false; break; end
 499       end
 500
 501       %# Call event only if a valid result has been found, therefore this
 502       %# code fragment has moved here. Stop integration if veventbreak is
 503       %# true
 504       if (vhaveeventfunction)
 505         vevent = ...
 506           odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
 507             {vu(:), vhmat}, [], vfunarguments{:});
 508         if (~isempty (vevent{1}) && vevent{1} == 1)
 509           vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
 510           vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
 511           vunhandledtermination = false; break;
 512         end
 513       end
 514     end %# If the error is acceptable ...
 515
 516     %# Update the step size for the next integration step
 517     if (~vstepsizegiven)
 518       %# vdelta may be 0 or even negative - could be an iteration problem
 519       vdelta = max (vdelta, eps);
 520       vstepsize = min (vodeoptions.MaxStep, ...
 521         min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
 522     elseif (vstepsizegiven)
 523       if (vcntloop < vtimelength)
 524         vstepsize = vslot(1,vcntloop-1) - vslot(1,vcntloop-2);
 525       end
 526     end
 527
 528     %# Update counters that count the number of iteration cycles
 529     vcntcycles = vcntcycles + 1; %# Needed for postprocessing
 530     vcntiter = vcntiter + 1;     %# Needed to find iteration problems
 531
 532     %# Stop solving because the last 1000 steps no successful valid
 533     %# value has been found
 534     if (vcntiter >= 5000)
 535       error (['Solving has not been successful. The iterative', ...
 536         ' integration loop exited at time t = %f before endpoint at', ...
 537         ' tend = %f was reached. This happened because the iterative', ...
 538         ' integration loop does not find a valid solution at this time', ...
 539         ' stamp. Try to reduce the value of "InitialStep" and/or', ...
 540         ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
 541     end
 542
 543   end %# The main loop
 544
 545   %# Check if integration of the ode has been successful
 546   if (vtimestamp < vtimestop)
 547     if (vunhandledtermination == true)
 548       error (['Solving has not been successful. The iterative', ...
 549         ' integration loop exited at time t = %f', ...
 550         ' before endpoint at tend = %f was reached. This may', ...
 551         ' happen if the stepsize grows smaller than defined in', ...
 552         ' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
 553         ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
 554     else
 555       warning ('OdePkg:HideWarning', ...
 556         ['Solver has been stopped by a call of "break" in', ...
 557          ' the main iteration loop at time t = %f before endpoint at', ...
 558          ' tend = %f was reached. This may happen because the @odeplot', ...
 559          ' function returned "true" or the @event function returned "true".'], ...
 560          vtimestamp, vtimestop);
 561     end
 562   end
 563
 564   %# Postprocessing, do whatever when terminating integration algorithm
 565   if (vhaveoutputfunction) %# Cleanup plotter
 566     feval (vodeoptions.OutputFcn, vtimestamp, ...
 567       vretvalresult(vcntloop-1,:)', 'done', vfunarguments{:});
 568   end
 569   if (vhaveeventfunction)  %# Cleanup event function handling
 570     odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
 571       {vretvalresult(vcntloop-1,:), vhmat}, 'done', vfunarguments{:});
 572   end
 573
 574   %# Print additional information if option Stats is set
 575   if (strcmp (vodeoptions.Stats, 'on'))
 576     vhavestats = true;
 577     vnsteps    = vcntloop-2;                    %# vcntloop from 2..end
 578     vnfailed   = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
 579     vnfevals   = 13*(vcntcycles-1);             %# number of ode evaluations
 580     vndecomps  = 0;                             %# number of LU decompositions
 581     vnpds      = 0;                             %# number of partial derivatives
 582     vnlinsols  = 0;                             %# no. of solutions of linear systems
 583     %# Print cost statistics if no output argument is given
 584     if (nargout == 0)
 585       vmsg = fprintf (1, 'Number of successful steps: %d', vnsteps);
 586       vmsg = fprintf (1, 'Number of failed attempts:  %d', vnfailed);
 587       vmsg = fprintf (1, 'Number of function calls:   %d', vnfevals);
 588     end
 589   else vhavestats = false;
 590   end
 591
 592   if (nargout == 1)                 %# Sort output variables, depends on nargout
 593     varargout{1}.x = vretvaltime;   %# Time stamps are saved in field x
 594     varargout{1}.y = vretvalresult; %# Results are saved in field y
 595     varargout{1}.solver = 'ode78d';  %# Solver name is saved in field solver
 596     if (vhaveeventfunction)
 597       varargout{1}.ie = vevent{2};  %# Index info which event occured
 598       varargout{1}.xe = vevent{3};  %# Time info when an event occured
 599       varargout{1}.ye = vevent{4};  %# Results when an event occured
 600     end
 601     if (vhavestats)
 602       varargout{1}.stats = struct;
 603       varargout{1}.stats.nsteps   = vnsteps;
 604       varargout{1}.stats.nfailed  = vnfailed;
 605       varargout{1}.stats.nfevals  = vnfevals;
 606       varargout{1}.stats.npds     = vnpds;
 607       varargout{1}.stats.ndecomps = vndecomps;
 608       varargout{1}.stats.nlinsols = vnlinsols;
 609     end
 610   elseif (nargout == 2)
 611     varargout{1} = vretvaltime;     %# Time stamps are first output argument
 612     varargout{2} = vretvalresult;   %# Results are second output argument
 613   elseif (nargout == 5)
 614     varargout{1} = vretvaltime;     %# Same as (nargout == 2)
 615     varargout{2} = vretvalresult;   %# Same as (nargout == 2)
 616     varargout{3} = [];              %# LabMat doesn't accept lines like
 617     varargout{4} = [];              %# varargout{3} = varargout{4} = [];
 618     varargout{5} = [];
 619     if (vhaveeventfunction)
 620       varargout{3} = vevent{3};     %# Time info when an event occured
 621       varargout{4} = vevent{4};     %# Results when an event occured
 622       varargout{5} = vevent{2};     %# Index info which event occured
 623     end
 624   %# else nothing will be returned, varargout{1} undefined
 625   end
 626
 627 %! # We are using a "pseudo-DDE" implementation for all tests that
 628 %! # are done for this function. We also define an Events and a
 629 %! # pseudo-Mass implementation. For further tests we also define a
 630 %! # reference solution (computed at high accuracy) and an OutputFcn.
 631 %!function [vyd] = fexp (vt, vy, vz, varargin)
 632 %!  vyd(1,1) = exp (- vt) - vz(1); %# The DDEs that are
 633 %!  vyd(2,1) = vy(1) - vz(2);      %# used for all examples
 634 %!function [vval, vtrm, vdir] = feve (vt, vy, vz, varargin)
 635 %!  vval = fexp (vt, vy, vz); %# We use the derivatives
 636 %!  vtrm = zeros (2,1);       %# don't stop solving here
 637 %!  vdir = ones (2,1);        %# in positive direction
 638 %!function [vval, vtrm, vdir] = fevn (vt, vy, vz, varargin)
 639 %!  vval = fexp (vt, vy, vz); %# We use the derivatives
 640 %!  vtrm = ones (2,1);        %# stop solving here
 641 %!  vdir = ones (2,1);        %# in positive direction
 642 %!function [vmas] = fmas (vt, vy, vz, varargin)
 643 %!  vmas =  [1, 0; 0, 1];     %# Dummy mass matrix for tests
 644 %!function [vmas] = fmsa (vt, vy, vz, varargin)
 645 %!  vmas = sparse ([1, 0; 0, 1]); %# A dummy sparse matrix
 646 %!function [vref] = fref ()       %# The reference solution
 647 %!  vref = [0.12194462133618, 0.01652432423938];
 648 %!function [vout] = fout (vt, vy, vflag, varargin)
 649 %!  if (regexp (char (vflag), 'init') == 1)
 650 %!    if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
 651 %!  elseif (isempty (vflag))
 652 %!    if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
 653 %!    vout = false;
 654 %!  elseif (regexp (char (vflag), 'done') == 1)
 655 %!    if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
 656 %!  else error ('"fout" invalid vflag');
 657 %!  end
 658 %!
 659 %! %# Turn off output of warning messages for all tests, turn them on
 660 %! %# again if the last test is called
 661 %!error %# input argument number one
 662 %!  warning ('off', 'OdePkg:InvalidOption');
 663 %!  B = ode78d (1, [0 5], [1; 0], 1, [1; 0]);
 664 %!error %# input argument number two
 665 %!  B = ode78d (@fexp, 1, [1; 0], 1, [1; 0]);
 666 %!error %# input argument number three
 667 %!  B = ode78d (@fexp, [0 5], 1, 1, [1; 0]);
 668 %!error %# input argument number four
 669 %!  B = ode78d (@fexp, [0 5], [1; 0], [1; 1], [1; 0]);
 670 %!error %# input argument number five
 671 %!  B = ode78d (@fexp, [0 5], [1; 0], 1, 1);
 672 %!test %# one output argument
 673 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0]);
 674 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 675 %!  assert (isfield (vsol, 'solver'));
 676 %!  assert (vsol.solver, 'ode78d');
 677 %!test %# two output arguments
 678 %!  [vt, vy] = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0]);
 679 %!  assert ([vt(end), vy(end,:)], [5, fref], 0.2);
 680 %!test %# five output arguments and no Events
 681 %!  [vt, vy, vxe, vye, vie] = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0]);
 682 %!  assert ([vt(end), vy(end,:)], [5, fref], 0.2);
 683 %!  assert ([vie, vxe, vye], []);
 684 %!test %# anonymous function instead of real function
 685 %!  faym = @(vt, vy, vz) [exp(-vt) - vz(1); vy(1) - vz(2)];
 686 %!  vsol = ode78d (faym, [0 5], [1; 0], 1, [1; 0]);
 687 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 688 %!test %# extra input arguments passed trhough
 689 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], 'KL');
 690 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 691 %!test %# empty OdePkg structure *but* extra input arguments
 692 %!  vopt = odeset;
 693 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt, 12, 13, 'KL');
 694 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 695 %!error %# strange OdePkg structure
 696 %!  vopt = struct ('foo', 1);
 697 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 698 %!test %# AbsTol option
 699 %!  vopt = odeset ('AbsTol', 1e-5);
 700 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 701 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 702 %!test %# AbsTol and RelTol option
 703 %!  vopt = odeset ('AbsTol', 1e-7, 'RelTol', 1e-7);
 704 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 705 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 706 %!test %# RelTol and NormControl option
 707 %!  vopt = odeset ('AbsTol', 1e-7, 'NormControl', 'on');
 708 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 709 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 710 %!test %# NonNegative for second component
 711 %!  vopt = odeset ('NonNegative', 1);
 712 %!  vsol = ode78d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
 713 %!  assert ([vsol.x(end), vsol.y(end,:)], [2.5, 0.001, 0.237], 0.2);
 714 %!test %# Details of OutputSel and Refine can't be tested
 715 %!  vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5);
 716 %!  vsol = ode78d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
 717 %!test %# Stats must add further elements in vsol
 718 %!  vopt = odeset ('Stats', 'on');
 719 %!  vsol = ode78d (@fexp, [0 2.5], [1; 0], 1, [1; 0], vopt);
 720 %!  assert (isfield (vsol, 'stats'));
 721 %!  assert (isfield (vsol.stats, 'nsteps'));
 722 %!test %# InitialStep option
 723 %!  vopt = odeset ('InitialStep', 1e-8);
 724 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 725 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 726 %!test %# MaxStep option
 727 %!  vopt = odeset ('MaxStep', 1e-2);
 728 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 729 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 730 %!test %# Events option add further elements in vsol
 731 %!  vopt = odeset ('Events', @feve);
 732 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 733 %!  assert (isfield (vsol, 'ie'));
 734 %!  assert (vsol.ie, [1; 1]);
 735 %!  assert (isfield (vsol, 'xe'));
 736 %!  assert (isfield (vsol, 'ye'));
 737 %!test %# Events option, now stop integration
 738 %!  warning ('off', 'OdePkg:HideWarning');
 739 %!  vopt = odeset ('Events', @fevn, 'NormControl', 'on');
 740 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 741 %!  assert ([vsol.ie, vsol.xe, vsol.ye], ...
 742 %!    [1.0000, 2.9219, -0.2127, -0.2671], 0.2);
 743 %!test %# Events option, five output arguments
 744 %!  vopt = odeset ('Events', @fevn, 'NormControl', 'on');
 745 %!  [vt, vy, vxe, vye, vie] = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 746 %!  assert ([vie, vxe, vye], ...
 747 %!    [1.0000, 2.9219, -0.2127, -0.2671], 0.2);
 748 %!
 749 %! %# test for Jacobian option is missing
 750 %! %# test for Jacobian (being a sparse matrix) is missing
 751 %! %# test for JPattern option is missing
 752 %! %# test for Vectorized option is missing
 753 %! %# test for NewtonTol option is missing
 754 %! %# test for MaxNewtonIterations option is missing
 755 %!
 756 %!test %# Mass option as function
 757 %!  vopt = odeset ('Mass', eye (2,2));
 758 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 759 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 760 %!test %# Mass option as matrix
 761 %!  vopt = odeset ('Mass', eye (2,2));
 762 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 763 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 764 %!test %# Mass option as sparse matrix
 765 %!  vopt = odeset ('Mass', sparse (eye (2,2)));
 766 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 767 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 768 %!test %# Mass option as function and sparse matrix
 769 %!  vopt = odeset ('Mass', @fmsa);
 770 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 771 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 772 %!test %# Mass option as function and MStateDependence
 773 %!  vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong');
 774 %!  vsol = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 775 %!  assert ([vsol.x(end), vsol.y(end,:)], [5, fref], 0.2);
 776 %!test %# Set BDF option to something else than default
 777 %!  vopt = odeset ('BDF', 'on');
 778 %!  [vt, vy] = ode78d (@fexp, [0 5], [1; 0], 1, [1; 0], vopt);
 779 %!  assert ([vt(end), vy(end,:)], [5, fref], 0.5);
 780 %!
 781 %! %# test for MvPattern option is missing
 782 %! %# test for InitialSlope option is missing
 783 %! %# test for MaxOrder option is missing
 784 %!
 785 %!  warning ('on', 'OdePkg:InvalidOption');
 786
 787 %# Local Variables: ***
 788 %# mode: octave ***
 789 %# End: ***