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