octave_packages/odepkg-0.8.2/odebwe.m

   1 %# Copyright (C) 2009-2012, Sebastian Schoeps <schoeps AT math DOT uni-wuppertal DOT de>
   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{}] =} odebwe (@var{@@fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}])
  19 %# @deftypefnx {Command} {[@var{sol}] =} odebwe (@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}]] =} odebwe (@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 stiff ordinary differential equations (stiff ODEs) or stiff differential algebraic equations (stiff DAEs) with the Backward Euler method.
  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 %# vjac = @@(vt,vy) [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
  35 %# vopt = odeset ("RelTol", 1e-3, "AbsTol", 1e-3, \
  36 %#          "NormControl", "on", "OutputFcn", @@odeplot, \
  37 %#          "Jacobian",vjac);
  38 %# odebwe (fvdb, [0 20], [2 0], vopt);
  39 %# @end example
  40 %# @end deftypefn
  41 %#
  42 %# @seealso{odepkg}
  43
  44 function [varargout] = odebwe (vfun, vslot, vinit, varargin)
  45
  46   if (nargin == 0) %# Check number and types of all input arguments
  47     help ('odebwe');
  48     error ('OdePkg:InvalidArgument', ...
  49       'Number of input arguments must be greater than zero');
  50
  51   elseif (nargin < 3)
  52     print_usage;
  53
  54   elseif ~(isa (vfun, 'function_handle') || isa (vfun, 'inline'))
  55     error ('OdePkg:InvalidArgument', ...
  56       'First input argument must be a valid function handle');
  57
  58   elseif (~isvector (vslot) || length (vslot) < 2)
  59     error ('OdePkg:InvalidArgument', ...
  60       'Second input argument must be a valid vector');
  61
  62   elseif (~isvector (vinit) || ~isnumeric (vinit))
  63     error ('OdePkg:InvalidArgument', ...
  64       'Third input argument must be a valid numerical value');
  65
  66   elseif (nargin >= 4)
  67
  68     if (~isstruct (varargin{1}))
  69       %# varargin{1:len} are parameters for vfun
  70       vodeoptions = odeset;
  71       vfunarguments = varargin;
  72
  73     elseif (length (varargin) > 1)
  74       %# varargin{1} is an OdePkg options structure vopt
  75       vodeoptions = odepkg_structure_check (varargin{1}, 'odebwe');
  76       vfunarguments = {varargin{2:length(varargin)}};
  77
  78     else %# if (isstruct (varargin{1}))
  79       vodeoptions = odepkg_structure_check (varargin{1}, 'odebwe');
  80       vfunarguments = {};
  81
  82     end
  83
  84   else %# if (nargin == 3)
  85     vodeoptions = odeset;
  86     vfunarguments = {};
  87   end
  88
  89   %# Start preprocessing, have a look which options are set in
  90   %# vodeoptions, check if an invalid or unused option is set
  91   vslot = vslot(:).';     %# Create a row vector
  92   vinit = vinit(:).';     %# Create a row vector
  93   if (length (vslot) > 2) %# Step size checking
  94     vstepsizefixed = true;
  95   else
  96     vstepsizefixed = false;
  97   end
  98
  99   %# The adaptive method require a second estimate for
 100   %# the comparsion, while the fixed step size algorithm
 101   %# needs only one
 102   if ~vstepsizefixed
 103     vestimators = 2;
 104   else
 105     vestimators = 1;
 106   end
 107
 108   %# Get the default options that can be set with 'odeset' temporarily
 109   vodetemp = odeset;
 110
 111   %# Implementation of the option RelTol has been finished. This option
 112   %# can be set by the user to another value than default value.
 113   if (isempty (vodeoptions.RelTol) && ~vstepsizefixed)
 114     vodeoptions.RelTol = 1e-6;
 115     warning ('OdePkg:InvalidArgument', ...
 116       'Option "RelTol" not set, new value %f is used', vodeoptions.RelTol);
 117   elseif (~isempty (vodeoptions.RelTol) && vstepsizefixed)
 118     warning ('OdePkg:InvalidArgument', ...
 119       'Option "RelTol" will be ignored if fixed time stamps are given');
 120   end
 121
 122   %# Implementation of the option AbsTol has been finished. This option
 123   %# can be set by the user to another value than default value.
 124   if (isempty (vodeoptions.AbsTol) && ~vstepsizefixed)
 125     vodeoptions.AbsTol = 1e-6;
 126     warning ('OdePkg:InvalidArgument', ...
 127       'Option "AbsTol" not set, new value %f is used', vodeoptions.AbsTol);
 128   elseif (~isempty (vodeoptions.AbsTol) && vstepsizefixed)
 129     warning ('OdePkg:InvalidArgument', ...
 130       'Option "AbsTol" will be ignored if fixed time stamps are given');
 131   else
 132     vodeoptions.AbsTol = vodeoptions.AbsTol(:); %# Create column vector
 133   end
 134
 135   %# Implementation of the option NormControl has been finished. This
 136   %# option can be set by the user to another value than default value.
 137   if (strcmp (vodeoptions.NormControl, 'on')) vnormcontrol = true;
 138   else vnormcontrol = false; end
 139
 140   %# Implementation of the option OutputFcn has been finished. This
 141   %# option can be set by the user to another value than default value.
 142   if (isempty (vodeoptions.OutputFcn) && nargout == 0)
 143     vodeoptions.OutputFcn = @odeplot;
 144     vhaveoutputfunction = true;
 145   elseif (isempty (vodeoptions.OutputFcn)), vhaveoutputfunction = false;
 146   else vhaveoutputfunction = true;
 147   end
 148
 149   %# Implementation of the option OutputSel has been finished. This
 150   %# option can be set by the user to another value than default value.
 151   if (~isempty (vodeoptions.OutputSel)), vhaveoutputselection = true;
 152   else vhaveoutputselection = false; end
 153
 154   %# Implementation of the option OutputSave has been finished. This
 155   %# option can be set by the user to another value than default value.
 156   if (isempty (vodeoptions.OutputSave)), vodeoptions.OutputSave = 1;
 157   end
 158
 159   %# Implementation of the option Stats has been finished. This option
 160   %# can be set by the user to another value than default value.
 161
 162   %# Implementation of the option InitialStep has been finished. This
 163   %# option can be set by the user to another value than default value.
 164   if (isempty (vodeoptions.InitialStep) && ~vstepsizefixed)
 165     vodeoptions.InitialStep = (vslot(1,2) - vslot(1,1)) / 10;
 166     warning ('OdePkg:InvalidArgument', ...
 167       'Option "InitialStep" not set, new value %f is used', vodeoptions.InitialStep);
 168   end
 169
 170   %# Implementation of the option MaxNewtonIterations has been finished. This option
 171   %# can be set by the user to another value than default value.
 172   if isempty (vodeoptions.MaxNewtonIterations)
 173     vodeoptions.MaxNewtonIterations = 10;
 174     warning ('OdePkg:InvalidArgument', ...
 175       'Option "MaxNewtonIterations" not set, new value %f is used', vodeoptions.MaxNewtonIterations);
 176   end
 177
 178   %# Implementation of the option NewtonTol has been finished. This option
 179   %# can be set by the user to another value than default value.
 180   if isempty (vodeoptions.NewtonTol)
 181     vodeoptions.NewtonTol = 1e-7;
 182     warning ('OdePkg:InvalidArgument', ...
 183       'Option "NewtonTol" not set, new value %f is used', vodeoptions.NewtonTol);
 184   end
 185
 186   %# Implementation of the option MaxStep has been finished. This option
 187   %# can be set by the user to another value than default value.
 188   if (isempty (vodeoptions.MaxStep) && ~vstepsizefixed)
 189     vodeoptions.MaxStep = (vslot(1,2) - vslot(1,1)) / 10;
 190     warning ('OdePkg:InvalidArgument', ...
 191       'Option "MaxStep" not set, new value %f is used', vodeoptions.MaxStep);
 192   end
 193
 194   %# Implementation of the option Events has been finished. This option
 195   %# can be set by the user to another value than default value.
 196   if (~isempty (vodeoptions.Events)), vhaveeventfunction = true;
 197   else vhaveeventfunction = false; end
 198
 199   %# Implementation of the option Jacobian has been finished. This option
 200   %# can be set by the user to another value than default value.
 201   if (~isempty (vodeoptions.Jacobian) && isnumeric (vodeoptions.Jacobian))
 202     vhavejachandle = false; vjac = vodeoptions.Jacobian; %# constant jac
 203   elseif (isa (vodeoptions.Jacobian, 'function_handle'))
 204     vhavejachandle = true; %# jac defined by a function handle
 205   else %# no Jacobian - we will use numerical differentiation
 206     vhavejachandle = false;
 207   end
 208
 209   %# Implementation of the option Mass has been finished. This option
 210   %# can be set by the user to another value than default value.
 211   if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass))
 212     vhavemasshandle = false; vmass = vodeoptions.Mass; %# constant mass
 213   elseif (isa (vodeoptions.Mass, 'function_handle'))
 214     vhavemasshandle = true; %# mass defined by a function handle
 215   else %# no mass matrix - creating a diag-matrix of ones for mass
 216     vhavemasshandle = false; vmass = sparse (eye (length (vinit)) );
 217   end
 218
 219   %# Implementation of the option MStateDependence has been finished.
 220   %# This option can be set by the user to another value than default
 221   %# value.
 222   if (strcmp (vodeoptions.MStateDependence, 'none'))
 223     vmassdependence = false;
 224   else vmassdependence = true;
 225   end
 226
 227   %# Other options that are not used by this solver. Print a warning
 228   %# message to tell the user that the option(s) is/are ignored.
 229   if (~isequal (vodeoptions.NonNegative, vodetemp.NonNegative))
 230     warning ('OdePkg:InvalidArgument', ...
 231       'Option "NonNegative" will be ignored by this solver');
 232   end
 233   if (~isequal (vodeoptions.Refine, vodetemp.Refine))
 234     warning ('OdePkg:InvalidArgument', ...
 235       'Option "Refine" will be ignored by this solver');
 236   end
 237   if (~isequal (vodeoptions.JPattern, vodetemp.JPattern))
 238     warning ('OdePkg:InvalidArgument', ...
 239       'Option "JPattern" will be ignored by this solver');
 240   end
 241   if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized))
 242     warning ('OdePkg:InvalidArgument', ...
 243       'Option "Vectorized" will be ignored by this solver');
 244   end
 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 odebwe
 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   vdirection  = sign (vtimestop);     %# Flag for direction to solve
 271
 272   if (~vstepsizefixed)
 273     vstepsize = vodeoptions.InitialStep;
 274     vminstepsize = (vtimestop - vtimestamp) / (1/eps);
 275   else %# If step size is given then use the fixed time steps
 276     vstepsize = vslot(1,2) - vslot(1,1);
 277     vminstepsize = sign (vstepsize) * eps;
 278   end
 279
 280   vretvaltime = vtimestamp; %# first timestamp output
 281   vretvalresult = vinit; %# first solution output
 282
 283   %# Initialize the OutputFcn
 284   if (vhaveoutputfunction)
 285     if (vhaveoutputselection)
 286       vretout = vretvalresult(vodeoptions.OutputSel);
 287     else
 288       vretout = vretvalresult;
 289     end
 290     feval (vodeoptions.OutputFcn, vslot.', ...
 291       vretout.', 'init', vfunarguments{:});
 292   end
 293
 294   %# Initialize the EventFcn
 295   if (vhaveeventfunction)
 296     odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
 297       vretvalresult.', 'init', vfunarguments{:});
 298   end
 299
 300   %# Initialize parameters and counters
 301   vcntloop = 2; vcntcycles = 1; vu = vinit; vcntsave = 2;
 302   vunhandledtermination = true; vpow = 1/2; vnpds = 0;
 303   vcntiter = 0; vcntnewt = 0; vndecomps = 0; vnlinsols = 0;
 304
 305   %# the following option enables the simplified Newton method
 306   %# which evaluates the Jacobian only once instead of the
 307   %# standard method that updates the Jacobian in each iteration
 308   vsimplified = false; % or true
 309
 310   %# The solver main loop - stop if the endpoint has been reached
 311   while ((vdirection * (vtimestamp) < vdirection * (vtimestop)) && ...
 312          (vdirection * (vstepsize) >= vdirection * (vminstepsize)))
 313
 314     %# Hit the endpoint of the time slot exactely
 315     if ((vtimestamp + vstepsize) > vdirection * vtimestop)
 316       vstepsize = vtimestop - vdirection * vtimestamp;
 317     end
 318
 319     %# Run the time integration for each estimator
 320     %# from vtimestamp -> vtimestamp+vstepsize
 321     for j = 1:vestimators
 322       %# Initial value (result of the previous timestep)
 323       y0 = vu;
 324       %# Initial guess for Newton-Raphson
 325       y(j,:)  = vu;
 326       %# We do not use a higher order approximation for the
 327       %# comparsion, but two steps by the Backward Euler
 328       %# method
 329       for i = 1:j
 330         % Initialize the time stepping parameters
 331         vthestep = vstepsize / j;
 332         vthetime = vtimestamp + i*vthestep;
 333         vnewtit  = 1;
 334         vresnrm  = inf (1, vodeoptions.MaxNewtonIterations);
 335
 336         %# Start the Newton iteration
 337         while (vnewtit < vodeoptions.MaxNewtonIterations) && ...
 338               (vresnrm (vnewtit) > vodeoptions.NewtonTol)
 339
 340           %# Compute the Jacobian of the non-linear equation,
 341           %# that is the matrix pencil of the mass matrix and
 342           %# the right-hand-side's Jacobian. Perform a (sparse)
 343           %# LU-Decomposition afterwards.
 344           if ( (vnewtit==1) || (~vsimplified) )
 345             %# Get the mass matrix from the left-hand-side
 346             if (vhavemasshandle)   %# Handle only the dynamic mass matrix,
 347               if (vmassdependence) %# constant mass matrices have already
 348                 vmass = feval ...  %# been set before (if any)
 349                   (vodeoptions.Mass, vthetime, y(j,:)', vfunarguments{:});
 350               else                 %# if (vmassdependence == false)
 351                 vmass = feval ...  %# then we only have the time argument
 352                   (vodeoptions.Mass, y(j,:)', vfunarguments{:});
 353               end
 354             end
 355             %# Get the Jacobian of the right-hand-side's function
 356             if (vhavejachandle)    %# Handle only the dynamic jacobian
 357               vjac = feval(vodeoptions.Jacobian, vthetime,...
 358                 y(j,:)', vfunarguments{:});
 359             elseif isempty(vodeoptions.Jacobian) %# If no Jacobian is given
 360               vjac = feval(@jacobian, vfun, vthetime,y(j,:)',...
 361                 vfunarguments);    %# then we differentiate
 362             end
 363             vnpds = vnpds + 1;
 364             vfulljac  = vmass/vthestep - vjac;
 365             %# one could do a matrix decomposition of vfulljac here,
 366             %# but the choice of decomposition depends on the problem
 367             %# and therefore we use the backslash-operator in row 374
 368           end
 369
 370           %# Compute the residual
 371           vres = vmass/vthestep*(y(j,:)-y0)' - feval(vfun,vthetime,y(j,:)',vfunarguments{:});
 372           vresnrm(vnewtit+1) = norm(vres,inf);
 373           %# Solve the linear system
 374           y(j,:) = vfulljac\(-vres+vfulljac*y(j,:)');
 375           %# the backslash operator decomposes the matrix
 376           %# and solves the system in a single step.
 377           vndecomps = vndecomps + 1;
 378           vnlinsols = vnlinsols + 1;
 379           %# Prepare next iteration
 380           vnewtit = vnewtit + 1;
 381         end %# while Newton
 382
 383         %# Leave inner loop if Newton diverged
 384         if vresnrm(vnewtit)>vodeoptions.NewtonTol
 385           break;
 386         end
 387         %# Save intermediate solution as initial value
 388         %# for the next intermediate step
 389         y0 = y(j,:);
 390         %# Count all Newton iterations
 391         vcntnewt = vcntnewt + (vnewtit-1);
 392       end %# for steps
 393
 394       %# Leave outer loop if Newton diverged
 395       if vresnrm(vnewtit)>vodeoptions.NewtonTol
 396         break;
 397       end
 398     end %# for estimators
 399
 400     % if all Newton iterations converged
 401     if vresnrm(vnewtit)<vodeoptions.NewtonTol
 402       %# First order approximation using step size h
 403       y1 = y(1,:);
 404       %# If adaptive: first order approximation using step
 405       %# size h/2, if fixed: y1=y2=y3
 406       y2 = y(vestimators,:);
 407       %# Second order approximation by ("Richardson")
 408       %# extrapolation using h and h/2
 409       y3 = y2 + (y2-y1);
 410     end
 411
 412     %# If Newton did not converge, repeat step with reduced
 413     %# step size, otherwise calculate the absolute local
 414     %# truncation error and the acceptable error
 415     if vresnrm(vnewtit)>vodeoptions.NewtonTol
 416       vdelta = 2; vtau = 1;
 417     elseif (~vstepsizefixed)
 418       if (~vnormcontrol)
 419         vdelta = abs (y3 - y1)';
 420         vtau = max (vodeoptions.RelTol * abs (vu.'), vodeoptions.AbsTol);
 421       else
 422         vdelta = norm ((y3 - y1)', Inf);
 423         vtau = max (vodeoptions.RelTol * max (norm (vu.', Inf), 1.0), ...
 424                     vodeoptions.AbsTol);
 425       end
 426     else %# if (vstepsizefixed == true)
 427       vdelta = 1; vtau = 2;
 428     end
 429
 430     %# If the error is acceptable then update the vretval variables
 431     if (all (vdelta <= vtau))
 432       vtimestamp = vtimestamp + vstepsize;
 433       vu = y2; % or y3 if we want the extrapolation....
 434
 435       %# Save the solution every vodeoptions.OutputSave steps
 436       if (mod (vcntloop-1,vodeoptions.OutputSave) == 0)
 437         vretvaltime(vcntsave,:) = vtimestamp;
 438         vretvalresult(vcntsave,:) = vu;
 439         vcntsave = vcntsave + 1;
 440       end
 441       vcntloop = vcntloop + 1; vcntiter = 0;
 442
 443       %# Call plot only if a valid result has been found, therefore this
 444       %# code fragment has moved here. Stop integration if plot function
 445       %# returns false
 446       if (vhaveoutputfunction)
 447         if (vhaveoutputselection)
 448           vpltout = vu(vodeoptions.OutputSel);
 449         else
 450           vpltout = vu;
 451         end
 452         vpltret = feval (vodeoptions.OutputFcn, vtimestamp, ...
 453           vpltout.', [], vfunarguments{:});
 454         if vpltret %# Leave loop
 455           vunhandledtermination = false; break;
 456         end
 457       end
 458
 459       %# Call event only if a valid result has been found, therefore this
 460       %# code fragment has moved here. Stop integration if veventbreak is
 461       %# true
 462       if (vhaveeventfunction)
 463         vevent = ...
 464           odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
 465             vu(:), [], vfunarguments{:});
 466         if (~isempty (vevent{1}) && vevent{1} == 1)
 467           vretvaltime(vcntloop-1,:) = vevent{3}(end,:);
 468           vretvalresult(vcntloop-1,:) = vevent{4}(end,:);
 469           vunhandledtermination = false; break;
 470         end
 471       end
 472     end %# If the error is acceptable ...
 473
 474     %# Update the step size for the next integration step
 475     if (~vstepsizefixed)
 476       %# 20080425, reported by Marco Caliari
 477       %# vdelta cannot be negative (because of the absolute value that
 478       %# has been introduced) but it could be 0, then replace the zeros
 479       %# with the maximum value of vdelta
 480       vdelta(find (vdelta == 0)) = max (vdelta);
 481       %# It could happen that max (vdelta) == 0 (ie. that the original
 482       %# vdelta was 0), in that case we double the previous vstepsize
 483       vdelta(find (vdelta == 0)) = max (vtau) .* (0.4 ^ (1 / vpow));
 484
 485       if (vdirection == 1)
 486         vstepsize = min (vodeoptions.MaxStep, ...
 487            min (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
 488       else
 489         vstepsize = max (vodeoptions.MaxStep, ...
 490           max (0.8 * vstepsize * (vtau ./ vdelta) .^ vpow));
 491       end
 492
 493     else %# if (vstepsizefixed)
 494       if (vresnrm(vnewtit)>vodeoptions.NewtonTol)
 495         vunhandledtermination = true;
 496         break;
 497       elseif (vcntloop <= vtimelength)
 498         vstepsize = vslot(vcntloop) - vslot(vcntloop-1);
 499       else %# Get out of the main integration loop
 500         break;
 501       end
 502     end
 503
 504     %# Update counters that count the number of iteration cycles
 505     vcntcycles = vcntcycles + 1; %# Needed for cost statistics
 506     vcntiter = vcntiter + 1;     %# Needed to find iteration problems
 507
 508     %# Stop solving because the last 1000 steps no successful valid
 509     %# value has been found
 510     if (vcntiter >= 5000)
 511       error (['Solving has not been successful. The iterative', ...
 512         ' integration loop exited at time t = %f before endpoint at', ...
 513         ' tend = %f was reached. This happened because the iterative', ...
 514         ' integration loop does not find a valid solution at this time', ...
 515         ' stamp. Try to reduce the value of "InitialStep" and/or', ...
 516         ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
 517     end
 518
 519   end %# The main loop
 520
 521   %# Check if integration of the ode has been successful
 522   if (vdirection * vtimestamp < vdirection * vtimestop)
 523     if (vunhandledtermination == true)
 524       error ('OdePkg:InvalidArgument', ...
 525         ['Solving has not been successful. The iterative', ...
 526          ' integration loop exited at time t = %f', ...
 527          ' before endpoint at tend = %f was reached. This may', ...
 528          ' happen if the stepsize grows smaller than defined in', ...
 529          ' vminstepsize. Try to reduce the value of "InitialStep" and/or', ...
 530          ' "MaxStep" with the command "odeset".\n'], vtimestamp, vtimestop);
 531     else
 532       warning ('OdePkg:InvalidArgument', ...
 533         ['Solver has been stopped by a call of "break" in', ...
 534          ' the main iteration loop at time t = %f before endpoint at', ...
 535          ' tend = %f was reached. This may happen because the @odeplot', ...
 536          ' function returned "true" or the @event function returned "true".'], ...
 537          vtimestamp, vtimestop);
 538     end
 539   end
 540
 541   %# Postprocessing, do whatever when terminating integration algorithm
 542   if (vhaveoutputfunction) %# Cleanup plotter
 543     feval (vodeoptions.OutputFcn, vtimestamp, ...
 544       vu.', 'done', vfunarguments{:});
 545   end
 546   if (vhaveeventfunction)  %# Cleanup event function handling
 547     odepkg_event_handle (vodeoptions.Events, vtimestamp, ...
 548       vu.', 'done', vfunarguments{:});
 549   end
 550   %# Save the last step, if not already saved
 551   if (mod (vcntloop-2,vodeoptions.OutputSave) ~= 0)
 552     vretvaltime(vcntsave,:) = vtimestamp;
 553     vretvalresult(vcntsave,:) = vu;
 554   end
 555
 556   %# Print additional information if option Stats is set
 557   if (strcmp (vodeoptions.Stats, 'on'))
 558     vhavestats = true;
 559     vnsteps    = vcntloop-2;                    %# vcntloop from 2..end
 560     vnfailed   = (vcntcycles-1)-(vcntloop-2)+1; %# vcntcycl from 1..end
 561     vnfevals   = vcntnewt;                      %# number of rhs evaluations
 562     if isempty(vodeoptions.Jacobian)            %# additional evaluations due
 563       vnfevals = vnfevals + vcntnewt*(1+length(vinit)); %# to differentiation
 564     end
 565     %# Print cost statistics if no output argument is given
 566     if (nargout == 0)
 567       vmsg = fprintf (1, 'Number of successful steps: %d\n', vnsteps);
 568       vmsg = fprintf (1, 'Number of failed attempts:  %d\n', vnfailed);
 569       vmsg = fprintf (1, 'Number of function calls:   %d\n', vnfevals);
 570     end
 571   else
 572     vhavestats = false;
 573   end
 574
 575   if (nargout == 1)                 %# Sort output variables, depends on nargout
 576     varargout{1}.x = vretvaltime;   %# Time stamps are saved in field x
 577     varargout{1}.y = vretvalresult; %# Results are saved in field y
 578     varargout{1}.solver = 'odebwe'; %# Solver name is saved in field solver
 579     if (vhaveeventfunction)
 580       varargout{1}.ie = vevent{2};  %# Index info which event occured
 581       varargout{1}.xe = vevent{3};  %# Time info when an event occured
 582       varargout{1}.ye = vevent{4};  %# Results when an event occured
 583     end
 584     if (vhavestats)
 585       varargout{1}.stats = struct;
 586       varargout{1}.stats.nsteps   = vnsteps;
 587       varargout{1}.stats.nfailed  = vnfailed;
 588       varargout{1}.stats.nfevals  = vnfevals;
 589       varargout{1}.stats.npds     = vnpds;
 590       varargout{1}.stats.ndecomps = vndecomps;
 591       varargout{1}.stats.nlinsols = vnlinsols;
 592     end
 593   elseif (nargout == 2)
 594     varargout{1} = vretvaltime;     %# Time stamps are first output argument
 595     varargout{2} = vretvalresult;   %# Results are second output argument
 596   elseif (nargout == 5)
 597     varargout{1} = vretvaltime;     %# Same as (nargout == 2)
 598     varargout{2} = vretvalresult;   %# Same as (nargout == 2)
 599     varargout{3} = [];              %# LabMat doesn't accept lines like
 600     varargout{4} = [];              %# varargout{3} = varargout{4} = [];
 601     varargout{5} = [];
 602     if (vhaveeventfunction)
 603       varargout{3} = vevent{3};     %# Time info when an event occured
 604       varargout{4} = vevent{4};     %# Results when an event occured
 605       varargout{5} = vevent{2};     %# Index info which event occured
 606     end
 607   end
 608 end
 609
 610 function df = jacobian(vfun,vthetime,vtheinput,vfunarguments);
 611 %# Internal function for the numerical approximation of a jacobian
 612   vlen   = length(vtheinput);
 613   vsigma = sqrt(eps);
 614   vfun0 = feval(vfun,vthetime,vtheinput,vfunarguments{:});
 615   df = zeros(vlen,vlen);
 616   for j = 1:vlen
 617      vbuffer = vtheinput(j);
 618      if (vbuffer==0)
 619         vh = vsigma;
 620      elseif (abs(vbuffer)>1)
 621         vh = vsigma*vbuffer;
 622      else
 623         vh = sign(vbuffer)*vsigma;
 624      end
 625      vtheinput(j) = vbuffer + vh;
 626      df(:,j) = (feval(vfun,vthetime,vtheinput,...
 627                 vfunarguments{:}) - vfun0) / vh;
 628      vtheinput(j) = vbuffer;
 629   end
 630 end
 631 %! # We are using the "Van der Pol" implementation for all tests that
 632 %! # are done for this function. We also define a Jacobian, Events,
 633 %! # pseudo-Mass implementation. For further tests we also define a
 634 %! # reference solution (computed at high accuracy) and an OutputFcn
 635 %!function [ydot] = fpol (vt, vy, varargin) %# The Van der Pol
 636 %!  ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
 637 %!function [vjac] = fjac (vt, vy, varargin) %# its Jacobian
 638 %!  vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2];
 639 %!function [vjac] = fjcc (vt, vy, varargin) %# sparse type
 640 %!  vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]);
 641 %!function [vval, vtrm, vdir] = feve (vt, vy, varargin)
 642 %!  vval = fpol (vt, vy, varargin); %# We use the derivatives
 643 %!  vtrm = zeros (2,1);             %# that's why component 2
 644 %!  vdir = ones (2,1);              %# seems to not be exact
 645 %!function [vval, vtrm, vdir] = fevn (vt, vy, varargin)
 646 %!  vval = fpol (vt, vy, varargin); %# We use the derivatives
 647 %!  vtrm = ones (2,1);              %# that's why component 2
 648 %!  vdir = ones (2,1);              %# seems to not be exact
 649 %!function [vmas] = fmas (vt, vy)
 650 %!  vmas = [1, 0; 0, 1];            %# Dummy mass matrix for tests
 651 %!function [vmas] = fmsa (vt, vy)
 652 %!  vmas = sparse ([1, 0; 0, 1]);   %# A sparse dummy matrix
 653 %!function [vref] = fref ()         %# The computed reference sol
 654 %!  vref = [0.32331666704577, -1.83297456798624];
 655 %!function [vout] = fout (vt, vy, vflag, varargin)
 656 %!  if (regexp (char (vflag), 'init') == 1)
 657 %!    if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end
 658 %!  elseif (isempty (vflag))
 659 %!    if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end
 660 %!    vout = false;
 661 %!  elseif (regexp (char (vflag), 'done') == 1)
 662 %!    if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end
 663 %!  else error ('"fout" invalid vflag');
 664 %!  end
 665 %!
 666 %! %# Turn off output of warning messages for all tests, turn them on
 667 %! %# again if the last test is called
 668 %!error %# input argument number one
 669 %!  warning ('off', 'OdePkg:InvalidArgument');
 670 %!  B = odebwe (1, [0 25], [3 15 1]);
 671 %!error %# input argument number two
 672 %!  B = odebwe (@fpol, 1, [3 15 1]);
 673 %!error %# input argument number three
 674 %!  B = odebwe (@flor, [0 25], 1);
 675 %!test %# one output argument
 676 %!  vsol = odebwe (@fpol, [0 2], [2 0]);
 677 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 678 %!  assert (isfield (vsol, 'solver'));
 679 %!  assert (vsol.solver, 'odebwe');
 680 %!test %# two output arguments
 681 %!  [vt, vy] = odebwe (@fpol, [0 2], [2 0]);
 682 %!  assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
 683 %!test %# five output arguments and no Events
 684 %!  [vt, vy, vxe, vye, vie] = odebwe (@fpol, [0 2], [2 0]);
 685 %!  assert ([vt(end), vy(end,:)], [2, fref], 1e-3);
 686 %!  assert ([vie, vxe, vye], []);
 687 %!test %# anonymous function instead of real function
 688 %!  fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
 689 %!  vsol = odebwe (fvdb, [0 2], [2 0]);
 690 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 691 %!test %# extra input arguments passed trhough
 692 %!  vsol = odebwe (@fpol, [0 2], [2 0], 12, 13, 'KL');
 693 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 694 %!test %# empty OdePkg structure *but* extra input arguments
 695 %!  vopt = odeset;
 696 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL');
 697 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 698 %!error %# strange OdePkg structure
 699 %!  vopt = struct ('foo', 1);
 700 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 701 %!test %# Solve vdp in fixed step sizes
 702 %!  vsol = odebwe (@fpol, [0:0.001:2], [2 0]);
 703 %!  assert (vsol.x(:), [0:0.001:2]');
 704 %!  assert (vsol.y(end,:), fref, 1e-2);
 705 %!test %# Solve in backward direction starting at t=0
 706 %! %# vref = [-1.2054034414, 0.9514292694];
 707 %!  vsol = odebwe (@fpol, [0 -2], [2 0]);
 708 %! %# assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3);
 709 %!test %# Solve in backward direction starting at t=2
 710 %! %#  vref = [-1.2154183302, 0.9433018000];
 711 %!  vsol = odebwe (@fpol, [2 -2], [0.3233166627 -1.8329746843]);
 712 %! %#  assert ([vsol.x(end), vsol.y(end,:)], [-2, fref], 1e-3);
 713 %!test %# AbsTol option
 714 %!  vopt = odeset ('AbsTol', 1e-5);
 715 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 716 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-2);
 717 %!test %# AbsTol and RelTol option
 718 %!  vopt = odeset ('AbsTol', 1e-6, 'RelTol', 1e-6);
 719 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 720 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-2);
 721 %!test %# RelTol and NormControl option -- higher accuracy
 722 %!  vopt = odeset ('RelTol', 1e-6, 'NormControl', 'on');
 723 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 724 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 725 %!
 726 %! %# test for NonNegative option is missing
 727 %! %# test for OutputSel and Refine option is missing
 728 %!
 729 %!test %# Details of OutputSave can't be tested
 730 %!  vopt = odeset ('OutputSave', 1, 'OutputSel', 1);
 731 %!  vsla = odebwe (@fpol, [0 2], [2 0], vopt);
 732 %!  vopt = odeset ('OutputSave', 2);
 733 %!  vslb = odebwe (@fpol, [0 2], [2 0], vopt);
 734 %!  assert (length (vsla.x) > length (vslb.x))
 735 %!test %# Stats must add further elements in vsol
 736 %!  vopt = odeset ('Stats', 'on');
 737 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 738 %!  assert (isfield (vsol, 'stats'));
 739 %!  assert (isfield (vsol.stats, 'nsteps'));
 740 %!test %# InitialStep option
 741 %!  vopt = odeset ('InitialStep', 1e-8);
 742 %!  vsol = odebwe (@fpol, [0 0.2], [2 0], vopt);
 743 %!  assert ([vsol.x(2)-vsol.x(1)], [1e-8], 1e-9);
 744 %!test %# MaxStep option
 745 %!  vopt = odeset ('MaxStep', 1e-2);
 746 %!  vsol = odebwe (@fpol, [0 0.2], [2 0], vopt);
 747 %!  assert ([vsol.x(5)-vsol.x(4)], [1e-2], 1e-2);
 748 %!test %# Events option add further elements in vsol
 749 %!  vopt = odeset ('Events', @feve);
 750 %!  vsol = odebwe (@fpol, [0 10], [2 0], vopt);
 751 %!  assert (isfield (vsol, 'ie'));
 752 %!  assert (vsol.ie, [2; 1; 2; 1]);
 753 %!  assert (isfield (vsol, 'xe'));
 754 %!  assert (isfield (vsol, 'ye'));
 755 %!test %# Events option, now stop integration
 756 %!  vopt = odeset ('Events', @fevn, 'NormControl', 'on');
 757 %!  vsol = odebwe (@fpol, [0 10], [2 0], vopt);
 758 %!  assert ([vsol.ie, vsol.xe, vsol.ye], ...
 759 %!    [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
 760 %!test %# Events option, five output arguments
 761 %!  vopt = odeset ('Events', @fevn, 'NormControl', 'on');
 762 %!  [vt, vy, vxe, vye, vie] = odebwe (@fpol, [0 10], [2 0], vopt);
 763 %!  assert ([vie, vxe, vye], ...
 764 %!    [2.0, 2.496110, -0.830550, -2.677589], 1e-3);
 765 %!test %# Jacobian option
 766 %!  vopt = odeset ('Jacobian', @fjac);
 767 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 768 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 769 %!test %# Jacobian option and sparse return value
 770 %!  vopt = odeset ('Jacobian', @fjcc);
 771 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 772 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 773 %!
 774 %! %# test for JPattern option is missing
 775 %! %# test for Vectorized option is missing
 776 %!
 777 %!test %# Mass option as function
 778 %!  vopt = odeset ('Mass', @fmas);
 779 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 780 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 781 %!test %# Mass option as matrix
 782 %!  vopt = odeset ('Mass', eye (2,2));
 783 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 784 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 785 %!test %# Mass option as sparse matrix
 786 %!  vopt = odeset ('Mass', sparse (eye (2,2)));
 787 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 788 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 789 %!test %# Mass option as function and sparse matrix
 790 %!  vopt = odeset ('Mass', @fmsa);
 791 %!  vsol = odebwe (@fpol, [0 2], [2 0], vopt);
 792 %!  assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3);
 793 %!
 794 %! %# test for MStateDependence option is missing
 795 %! %# test for MvPattern option is missing
 796 %! %# test for InitialSlope option is missing
 797 %! %# test for MaxOrder option is missing
 798 %! %# test for BDF option is missing
 799 %!
 800 %!  warning ('on', 'OdePkg:InvalidArgument');
 801
 802 %# Local Variables: ***
 803 %# mode: octave ***
 804 %# End: ***