Mercurial > octave-nkf
changeset 20568:fcb792acab9b
Moving ode45, odeset, odeget, and levenshtein from odepkg to core.
* libinterp/corefcn/levenshtein.cc: move function from odepkg into core
* libinterp/corefcn/module.mk: include levenshtein.cc
* scripts/ode: move ode45, odeset, odeget, and all dependencies
from odepkg into core
* scripts/module.mk: include them
* doc/interpreter/diffeq.txi: add documentation for ode45,
odeset, odeget
* NEWS: announce functions included with this changeset
* scripts/help/__unimplemented__.m: removed new functions
| author | jcorno <jacopo.corno@gmail.com> |
|---|---|
| date | Thu, 24 Sep 2015 12:58:46 +0200 |
| parents | 2480bbcd1333 |
| children | b70cc4bd8109 |
| files | NEWS doc/interpreter/diffeq.txi libinterp/corefcn/levenshtein.cc libinterp/corefcn/module.mk scripts/help/__unimplemented__.m scripts/module.mk scripts/ode/module.mk scripts/ode/ode45.m scripts/ode/odeget.m scripts/ode/odeset.m scripts/ode/private/AbsRel_Norm.m scripts/ode/private/fuzzy_compare.m scripts/ode/private/hermite_quartic_interpolation.m scripts/ode/private/integrate_adaptive.m scripts/ode/private/integrate_const.m scripts/ode/private/integrate_n_steps.m scripts/ode/private/kahan.m scripts/ode/private/ode_struct_value_check.m scripts/ode/private/odepkg_event_handle.m scripts/ode/private/odepkg_structure_check.m scripts/ode/private/runge_kutta_45_dorpri.m scripts/ode/private/starting_stepsize.m |
| diffstat | 22 files changed, 3898 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Thu Oct 01 13:37:03 2015 -0400 +++ b/NEWS Thu Sep 24 12:58:46 2015 +0200 @@ -43,6 +43,9 @@ ** Other new functions added in 4.2: psi + odeset + odeget + ode45 ** Deprecated functions.
--- a/doc/interpreter/diffeq.txi Thu Oct 01 13:37:03 2015 -0400 +++ b/doc/interpreter/diffeq.txi Thu Sep 24 12:58:46 2015 +0200 @@ -115,6 +115,41 @@ Octave distribution in the examples directory under the name @file{oregonator.m}. +@menu +* Matlab-compatible solvers:: +@end menu + +@node Matlab-compatible solvers +@subsection Matlab-compatible solvers + +Octave also provides a set of solvers for initial value problems for Ordinary +Differential Equations that have a Matlab-compatible interface. +The options for this class of methods are set using the functions. +@itemize + @item @code{odeset} + @item @code{odeget} +@end itemize + +Currently implemented solvers are: +@itemize + @item Runge-Kutta methods + @itemize + @item @code{ode45} Integrates a system of non--stiff ordinary differential equations + (non--stiff ODEs and DAEs) using second order Dormand-Prince method. + This is a fourth--order accurate integrator therefore the local error + normally expected is @math{O(h^5)}. This solver requires six function + evaluations per integration step. + @end itemize +@end itemize + + +@DOCSTRING(ode45) + +@DOCSTRING(odeset) + +@DOCSTRING(odeget) + + @node Differential-Algebraic Equations @section Differential-Algebraic Equations
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/libinterp/corefcn/levenshtein.cc Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,109 @@ +/* + +Copyright (C) 2014 Jacopo Corno +Copyright (C) 2014 Carlo de Falco +Copyright (C) 2013 Roberto Porcu +Copyright (C) 2013 Mattia Penati + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 3 of the License, or (at +your option) any later version. + +Octave is distributed in the hope that it will be useful, but +WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +General Public License for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun.h" + +#define MIN2(A,B) ((A)<(B)?(A):(B)) +#define MIN3(A,B,C) (MIN2(MIN2((A),(B)),(C))) + +DEFUN (levenshtein, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{d} =} levenshtein (@var{str1}, @var{str2}, [@var{u_bound}])\n \ +@deftypefnx {Loadable Function} {[@var{d}, @var{mat}] =} levenshtein (@var{str1}, @var{str2}, [@var{u_bound}])\n \ +This function file computes the Levenshtein distance between two strings. More details at @url{http://en.wikipedia.org/wiki/Levenshtein_distance}.\n \ +This function can be called with one or two output arguments: @var{d} is the distance between the two strings and @var{mat} is the matrix computed by Levenshtein algorithm.\n \ +The first and the second input arguments, @var{str1} and @var{str2}, are the two strings to be compared. This comparison is case-sensitive.\n \ +The third argument @var{u_bound} is optional and fixes an upper bound for the distance. If the distance is greater than this limit then the function ends and returns a value equal to Inf.\n \ +@seealso{odeset}\n\ +@end deftypefn") +{ + octave_value_list retval; + int nargin = args.length (); + octave_idx_type ub; + + if (nargin < 2 + || nargin > 3 + || nargout > 2) + { + print_usage (); + return retval; + } + + if (nargin == 3) + ub = args(2).idx_type_value (); + else + ub = std::numeric_limits<int32_t>::max (); + + Array<char> s1 (args(0).char_array_value ()); + char *s1_p = s1.fortran_vec (); + + Array<char> s2 (args(1).char_array_value ()); + char *s2_p = s2.fortran_vec (); + + octave_idx_type + s1len = s1.length (), + s2len = s2.length (), + ii, jj; + + if (! error_state) + { + Array<octave_idx_type> + dist (dim_vector (s1len + 1, s2len + 1), 0); + + for (ii = 1; ii <= s1len; ii++) + dist.xelem (ii, 0) = ii; + + for (jj = 1; jj <= s2len; jj++) + dist.xelem (0, jj) = jj; + + for (jj = 1; jj <= s2len; jj++) + { + for (ii = 1; ii <= s1len; ii++) + if (s1_p[ii-1] == s2_p[jj-1]) + dist.xelem (ii, jj) = dist.xelem (ii-1, jj-1); + else + dist.xelem (ii, jj) = + MIN3(dist.xelem (ii-1, jj) + 1, + dist.xelem (ii, jj-1) + 1, + dist.xelem (ii-1, jj-1) + 1); + + if (dist(MIN2(jj, s1len), jj) > ub) + { + retval(0) = std::numeric_limits<int32_t>::max (); + if (nargout == 2) + retval(1) = Matrix (); + return retval; + } + } + retval(0) = dist.xelem (s1len, s2len); + if (nargout == 2) + retval(1) = dist; + } + return retval; +}
--- a/libinterp/corefcn/module.mk Thu Oct 01 13:37:03 2015 -0400 +++ b/libinterp/corefcn/module.mk Thu Sep 24 12:58:46 2015 +0200 @@ -174,6 +174,7 @@ libinterp/corefcn/input.cc \ libinterp/corefcn/inv.cc \ libinterp/corefcn/kron.cc \ + libinterp/corefcn/levenshtein.cc \ libinterp/corefcn/load-path.cc \ libinterp/corefcn/load-save.cc \ libinterp/corefcn/lookup.cc \
--- a/scripts/help/__unimplemented__.m Thu Oct 01 13:37:03 2015 -0400 +++ b/scripts/help/__unimplemented__.m Thu Sep 24 12:58:46 2015 +0200 @@ -82,10 +82,9 @@ txt = ["matlabrc is not implemented. ", ... 'Octave uses the file ".octaverc" instead.']; - case {"ode113", "ode15i", "ode15s", "ode23", "ode23s", "ode23t", ... - "ode23tb", "ode45", "odeget", "odeset"} - txt = ["Octave provides lsode for solving differential equations. ", ... - "For more information try @code{help lsode}. ", ... + case {"ode113", "ode15i", "ode15s", "ode23", "ode23s", "ode23t", "ode23tb"} + txt = ["Octave provides lsode and ode45 for solving differential equations. ", ... + "For more information try @code{help lsode}, @code{help ode45}. ", ... "Matlab-compatible ODE functions are provided by the odepkg ", ... "package. See @url{http://octave.sourceforge.net/odepkg/}."];
--- a/scripts/module.mk Thu Oct 01 13:37:03 2015 -0400 +++ b/scripts/module.mk Thu Sep 24 12:58:46 2015 +0200 @@ -47,6 +47,7 @@ include scripts/java/module.mk include scripts/linear-algebra/module.mk include scripts/miscellaneous/module.mk +include scripts/ode/module.mk include scripts/optimization/module.mk include scripts/path/module.mk include scripts/pkg/module.mk
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/module.mk Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,39 @@ +FCN_FILE_DIRS += \ + scripts/ode \ + scripts/ode/private + +scripts_ode_PRIVATE_FCN_FILES = \ + scripts/ode/private/AbsRel_Norm.m \ + scripts/ode/private/fuzzy_compare.m \ + scripts/ode/private/hermite_quartic_interpolation.m \ + scripts/ode/private/integrate_adaptive.m \ + scripts/ode/private/integrate_const.m \ + scripts/ode/private/integrate_n_steps.m \ + scripts/ode/private/kahan.m \ + scripts/ode/private/odepkg_event_handle.m \ + scripts/ode/private/odepkg_structure_check.m \ + scripts/ode/private/ode_struct_value_check.m \ + scripts/ode/private/runge_kutta_45_dorpri.m \ + scripts/ode/private/starting_stepsize.m + + +scripts_ode_FCN_FILES = \ + scripts/ode/ode45.m \ + scripts/ode/odeset.m \ + scripts/ode/odeget.m + +scripts_odedir = $(fcnfiledir)/ode + +scripts_ode_DATA = $(scripts_ode_FCN_FILES) + +scripts_ode_privatedir = $(fcnfiledir)/ode/private + +scripts_ode_private_DATA = $(scripts_ode_PRIVATE_FCN_FILES) + +FCN_FILES += \ + $(scripts_ode_FCN_FILES) \ + $(scripts_ode_PRIVATE_FCN_FILES) + +PKG_ADD_FILES += scripts/ode/PKG_ADD + +DIRSTAMP_FILES += scripts/ode/$(octave_dirstamp)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/ode45.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,657 @@ +## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2014, Jacopo Corno <jacopo.corno@gmail.com> +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{sol}] =} ode45 (@var{fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}]) +## @deftypefnx {Function File} {[@var{t}, @var{y}, [@var{xe}, @var{ye}, @var{ie}]] =} ode45 (@var{fun}, @var{slot}, @var{init}, [@var{opt}], [@var{par1}, @var{par2}, @dots{}]) +## +## This function file can be used to solve a set of non--stiff ordinary +## differential equations (non--stiff ODEs) with the well known explicit +## Dormand-Prince method of order 4. +## +## This function can be called with two output arguments: @var{t} and @var{y}. +## Variable @var{t} is a column vector and contains the time stamps, instead +## @var{y} is a matrix in which each column refers to a different unknown of +## the problem and the rows number is the same of @var{t} rows number so that each +## row of @var{y} contains the values of all unknowns at the time value contained +## in the corresponding row in @var{t}. +## +## The first input argument must be a function_handle or an inline function that +## defines the set of ODE: @code{y' = f(t,y)}. As described above, this function +## must take two input arguments, where the first is the time and the second +## the unknowns, and must have just one output argument. +## +## The second input argument must contain time informations. Usually it should +## be a vector with at least two elements which define the initial and the final +## time instants; if the elements are more than two, then the solution will be +## evaluated also at these intermediate time instants unless the integrate function +## called is the @command{integrate_n_steps}. If there is only one time value, +## then it will give an error unless the options structure has no empty fields +## named @var{"TimeStepNumber"} and @var{"TimeStepSize"}. If the option +## @var{"TimeStepSize"} is not empty, then the stepper called will be +## @command{integrate_const}, if also @var{"TimeStepNumber"} is not empty it will +## be called the integrate function @command{integrate_n_steps}, otherwise it will +## be called @command{integrate_adaptive}. For this last possibility the user can +## set the tolerance for the timestep computation by setting a value to the option +## @var{"Tau"}, that as default value has @math{1.e-6}. +## +## The third input argument must contain the initial value for the unknown. +## If this is a vector then the solution @var{y} will be a matrix in which each +## column is the solution for the corresponding initial value in @var{init}. +## +## The fourth input argument is not mandatory and it should contain a structure +## with valid ODE fields. +## +## For example, solve an anonymous implementation of the Van der Pol equation +## @example +## fvdp = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; +## [T,Y] = ode45 (fvdp, [0 20], [2 0]); +## @end example +## @end deftypefn +## + +function [varargout] = ode45 (vfun, vslot, vinit, varargin) + + vorder = 5; % runge_kutta_45_dorpri uses local extrapolation + vsolver = "ode45"; + + if (nargin == 0) ## Check number and types of all input arguments + help (vsolver); + error ("OdePkg:InvalidArgument", ... + "number of input arguments must be greater than zero"); + endif + + if (nargin < 3) + print_usage; + endif + + if (nargin >= 4) + if (~isstruct (varargin{1})) + ## varargin{1:len} are parameters for vfun + vodeoptions = odeset; + vodeoptions.vfunarguments = varargin; + elseif (length (varargin) > 1) + ## varargin{1} is an OdePkg options structure vopt + vodeoptions = odepkg_structure_check (varargin{1}, "ode45"); + vodeoptions.vfunarguments = {varargin{2:length(varargin)}}; + else ## if (isstruct (varargin{1})) + vodeoptions = odepkg_structure_check (varargin{1}, "ode45"); + vodeoptions.vfunarguments = {}; + endif + else ## if (nargin == 3) + vodeoptions = odeset; + vodeoptions.vfunarguments = {}; + endif + + if (~isvector (vslot) || ~isnumeric (vslot)) + error ("OdePkg:InvalidArgument", ... + "second input argument must be a valid vector"); + endif + + if (length (vslot) < 2 && ... + (isempty (vodeoptions.TimeStepSize) ... + || isempty (vodeoptions.TimeStepNumber))) + error ("OdePkg:InvalidArgument", ... + "second input argument must be a valid vector"); + elseif (vslot(2) == vslot(1)) + error ("OdePkg:InvalidArgument", ... + "second input argument must be a valid vector"); + else + vodeoptions.vdirection = sign (vslot(2) - vslot(1)); + endif + vslot = vslot(:); + + if (~isvector (vinit) || ~isnumeric (vinit)) + error ("OdePkg:InvalidArgument", ... + "third input argument must be a valid numerical value"); + endif + vinit = vinit(:); + + if ~(isa (vfun, "function_handle") || isa (vfun, "inline")) + error ("OdePkg:InvalidArgument", ... + "first input argument must be a valid function handle"); + endif + + ## Start preprocessing, have a look which options are set in + ## vodeoptions, check if an invalid or unused option is set + if (isempty (vodeoptions.TimeStepNumber) ... + && isempty (vodeoptions.TimeStepSize)) + integrate_func = "adaptive"; + vodeoptions.vstepsizefixed = false; + elseif (~isempty (vodeoptions.TimeStepNumber) ... + && ~isempty (vodeoptions.TimeStepSize)) + integrate_func = "n_steps"; + vodeoptions.vstepsizefixed = true; + if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection) + warning ("OdePkg:InvalidArgument", ... + "option ''TimeStepSize'' has a wrong sign", ... + "it will be corrected automatically"); + vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize; + endif + elseif (isempty (vodeoptions.TimeStepNumber) && ~isempty (vodeoptions.TimeStepSize)) + integrate_func = "const"; + vodeoptions.vstepsizefixed = true; + if (sign (vodeoptions.TimeStepSize) != vodeoptions.vdirection) + warning ("OdePkg:InvalidArgument", ... + "option ''TimeStepSize'' has a wrong sign", ... + "it will be corrected automatically"); + vodeoptions.TimeStepSize = (-1)*vodeoptions.TimeStepSize; + endif + else + warning ("OdePkg:InvalidArgument", ... + "assuming an adaptive integrate function"); + integrate_func = "adaptive"; + endif + + ## Get the default options that can be set with "odeset" temporarily + vodetemp = odeset; + + ## Implementation of the option RelTol has been finished. This option + ## can be set by the user to another value than default value. + if (isempty (vodeoptions.RelTol) && ~vodeoptions.vstepsizefixed) + vodeoptions.RelTol = 1e-3; + warning ("OdePkg:InvalidArgument", ... + "Option ''RelTol'' not set, new value %f is used", ... + vodeoptions.RelTol); + elseif (~isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed) + warning ("OdePkg:InvalidArgument", ... + "Option ''RelTol'' will be ignored if fixed time stamps are given"); + endif + + ## Implementation of the option AbsTol has been finished. This option + ## can be set by the user to another value than default value. + if (isempty (vodeoptions.AbsTol) && ~vodeoptions.vstepsizefixed) + vodeoptions.AbsTol = 1e-6; + warning ("OdePkg:InvalidArgument", ... + "Option ''AbsTol'' not set, new value %f is used", ... + vodeoptions.AbsTol); + elseif (~isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed) + warning ("OdePkg:InvalidArgument", ... + "Option ''AbsTol'' will be ignored if fixed time stamps are given"); + else + vodeoptions.AbsTol = vodeoptions.AbsTol(:); ## Create column vector + endif + + ## Implementation of the option NormControl has been finished. This + ## option can be set by the user to another value than default value. + if (strcmp (vodeoptions.NormControl, "on")) + vodeoptions.vnormcontrol = true; + else + vodeoptions.vnormcontrol = false; + endif + + ## Implementation of the option NonNegative has been finished. This + ## option can be set by the user to another value than default value. + if (~isempty (vodeoptions.NonNegative)) + if (isempty (vodeoptions.Mass)) + vodeoptions.vhavenonnegative = true; + else + vodeoptions.vhavenonnegative = false; + warning ("OdePkg:InvalidArgument", ... + "Option 'NonNegative' will be ignored if mass matrix is set"); + endif + else + vodeoptions.vhavenonnegative = false; + endif + + ## Implementation of the option OutputFcn has been finished. This + ## option can be set by the user to another value than default value. + if (isempty (vodeoptions.OutputFcn) && nargout == 0) + vodeoptions.OutputFcn = @odeplot; + vodeoptions.vhaveoutputfunction = true; + elseif (isempty (vodeoptions.OutputFcn)) + vodeoptions.vhaveoutputfunction = false; + else + vodeoptions.vhaveoutputfunction = true; + endif + + ## Implementation of the option OutputSel has been finished. This + ## option can be set by the user to another value than default value. + if (~isempty (vodeoptions.OutputSel)) + vodeoptions.vhaveoutputselection = true; + else + vodeoptions.vhaveoutputselection = false; + endif + + ## Implementation of the option OutputSave has been finished. This + ## option can be set by the user to another value than default value. + if (isempty (vodeoptions.OutputSave)) + vodeoptions.OutputSave = 1; + endif + + ## Implementation of the option Refine has been finished. This option + ## can be set by the user to another value than default value. + if (vodeoptions.Refine > 0) + vodeoptions.vhaverefine = true; + else + vodeoptions.vhaverefine = false; + endif + + ## Implementation of the option Stats has been finished. This option + ## can be set by the user to another value than default value. + + ## Implementation of the option InitialStep has been finished. This + ## option can be set by the user to another value than default value. + if (isempty (vodeoptions.InitialStep) && strcmp (integrate_func, "adaptive")) + vodeoptions.InitialStep = vodeoptions.vdirection* ... + starting_stepsize (vorder, vfun, vslot(1), vinit, vodeoptions.AbsTol, ... + vodeoptions.RelTol, vodeoptions.vnormcontrol); + warning ("OdePkg:InvalidArgument", ... + "option ''InitialStep'' not set, estimated value %f is used", ... + vodeoptions.InitialStep); + elseif(isempty (vodeoptions.InitialStep)) + vodeoptions.InitialStep = odeget (vodeoptions, "TimeStepSize"); + endif + + ## Implementation of the option MaxStep has been finished. This option + ## can be set by the user to another value than default value. + if (isempty (vodeoptions.MaxStep) && ~vodeoptions.vstepsizefixed) + vodeoptions.MaxStep = abs (vslot(end) - vslot(1)) / 10; + warning ("OdePkg:InvalidArgument", ... + "Option ''MaxStep'' not set, new value %f is used", ... + vodeoptions.MaxStep); + endif + + ## Implementation of the option Events has been finished. This option + ## can be set by the user to another value than default value. + if (~isempty (vodeoptions.Events)) + vodeoptions.vhaveeventfunction = true; + else + vodeoptions.vhaveeventfunction = false; + endif + + ## The options 'Jacobian', 'JPattern' and 'Vectorized' will be ignored + ## by this solver because this solver uses an explicit Runge-Kutta + ## method and therefore no Jacobian calculation is necessary + if (~isequal (vodeoptions.Jacobian, vodetemp.Jacobian)) + warning ("OdePkg:InvalidArgument", ... + "option ''Jacobian'' will be ignored by this solver"); + endif + + if (~isequal (vodeoptions.JPattern, vodetemp.JPattern)) + warning ("OdePkg:InvalidArgument", ... + "option ''JPattern'' will be ignored by this solver"); + endif + + if (~isequal (vodeoptions.Vectorized, vodetemp.Vectorized)) + warning ("OdePkg:InvalidArgument", ... + "option ''Vectorized'' will be ignored by this solver"); + endif + + if (~isequal (vodeoptions.NewtonTol, vodetemp.NewtonTol)) + warning ("OdePkg:InvalidArgument", ... + "option ''NewtonTol'' will be ignored by this solver"); + endif + + if (~isequal (vodeoptions.MaxNewtonIterations,vodetemp.MaxNewtonIterations)) + warning ("OdePkg:InvalidArgument", ... + "option ''MaxNewtonIterations'' will be ignored by this solver"); + endif + + ## Implementation of the option Mass has been finished. This option + ## can be set by the user to another value than default value. + if (~isempty (vodeoptions.Mass) && isnumeric (vodeoptions.Mass)) + vhavemasshandle = false; + vmass = vodeoptions.Mass; ## constant mass + elseif (isa (vodeoptions.Mass, "function_handle")) + vhavemasshandle = true; ## mass defined by a function handle + else ## no mass matrix - creating a diag-matrix of ones for mass + vhavemasshandle = false; ## vmass = diag (ones (length (vinit), 1), 0); + endif + + ## Implementation of the option MStateDependence has been finished. + ## This option can be set by the user to another value than default + ## value. + if (strcmp (vodeoptions.MStateDependence, "none")) + vmassdependence = false; + else + vmassdependence = true; + endif + + ## Other options that are not used by this solver. Print a warning + ## message to tell the user that the option(s) is/are ignored. + if (~isequal (vodeoptions.MvPattern, vodetemp.MvPattern)) + warning ("OdePkg:InvalidArgument", ... + "option ''MvPattern'' will be ignored by this solver"); + endif + + if (~isequal (vodeoptions.MassSingular, vodetemp.MassSingular)) + warning ("OdePkg:InvalidArgument", ... + "option ''MassSingular'' will be ignored by this solver"); + endif + + if (~isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope)) + warning ("OdePkg:InvalidArgument", ... + "option ''InitialSlope'' will be ignored by this solver"); + endif + + if (~isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder)) + warning ("OdePkg:InvalidArgument", ... + "option ''MaxOrder'' will be ignored by this solver"); + endif + + if (~isequal (vodeoptions.BDF, vodetemp.BDF)) + warning ("OdePkg:InvalidArgument", ... + "option ''BDF'' will be ignored by this solver"); + endif + + ## Starting the initialisation of the core solver ode45 + SubOpts = vodeoptions; + + if (vhavemasshandle) ## Handle only the dynamic mass matrix, + if (vmassdependence) ## constant mass matrices have already + vmass = @(t,x) vodeoptions.Mass (t, x, vodeoptions.vfunarguments{:}); + vfun = @(t,x) vmass (t, x, vodeoptions.vfunarguments{:}) ... + \ vfun (t, x, vodeoptions.vfunarguments{:}); + else ## if (vmassdependence == false) + vmass = @(t) vodeoptions.Mass (t, vodeoptions.vfunarguments{:}); + vfun = @(t,x) vmass (t, vodeoptions.vfunarguments{:}) ... + \ vfun (t, x, vodeoptions.vfunarguments{:}); + endif + endif + + switch integrate_func + case "adaptive" + solution = integrate_adaptive (@runge_kutta_45_dorpri, ... + vorder, vfun, vslot, vinit, SubOpts); + case "n_steps" + solution = integrate_n_steps (@runge_kutta_45_dorpri, ... + vfun, vslot(1), vinit, ... + vodeoptions.TimeStepSize, ... + vodeoptions.TimeStepNumber, SubOpts); + case "const" + solution = integrate_const (@runge_kutta_45_dorpri, ... + vfun, vslot, vinit, ... + vodeoptions.TimeStepSize, SubOpts); + endswitch + + ## Postprocessing, do whatever when terminating integration algorithm + if (vodeoptions.vhaveoutputfunction) ## Cleanup plotter + feval (vodeoptions.OutputFcn, solution.t(end), ... + solution.x(end,:)', "done", vodeoptions.vfunarguments{:}); + endif + if (vodeoptions.vhaveeventfunction) ## Cleanup event function handling + odepkg_event_handle (vodeoptions.Events, solution.t(end), ... + solution.x(end,:)', "done", vodeoptions.vfunarguments{:}); + endif + + ## Print additional information if option Stats is set + if (strcmp (vodeoptions.Stats, "on")) + vhavestats = true; + vnsteps = solution.vcntloop-2; ## vcntloop from 2..end + vnfailed = (solution.vcntcycles-1)-(solution.vcntloop-2)+1; ## vcntcycl from 1..end + vnfevals = 7*(solution.vcntcycles-1); ## number of ode evaluations + vndecomps = 0; ## number of LU decompositions + vnpds = 0; ## number of partial derivatives + vnlinsols = 0; ## no. of solutions of linear systems + ## Print cost statistics if no output argument is given + if (nargout == 0) + vmsg = fprintf (1, "Number of successful steps: %d\n", vnsteps); + vmsg = fprintf (1, "Number of failed attempts: %d\n", vnfailed); + vmsg = fprintf (1, "Number of function calls: %d\n", vnfevals); + endif + else + vhavestats = false; + endif + + if (nargout == 1) ## Sort output variables, depends on nargout + varargout{1}.x = solution.t; ## Time stamps are saved in field x + varargout{1}.y = solution.x; ## Results are saved in field y + varargout{1}.solver = vsolver; ## Solver name is saved in field solver + if (vodeoptions.vhaveeventfunction) + varargout{1}.ie = solution.vevent{2}; ## Index info which event occured + varargout{1}.xe = solution.vevent{3}; ## Time info when an event occured + varargout{1}.ye = solution.vevent{4}; ## Results when an event occured + endif + if (vhavestats) + varargout{1}.stats = struct; + varargout{1}.stats.nsteps = vnsteps; + varargout{1}.stats.nfailed = vnfailed; + varargout{1}.stats.nfevals = vnfevals; + varargout{1}.stats.npds = vnpds; + varargout{1}.stats.ndecomps = vndecomps; + varargout{1}.stats.nlinsols = vnlinsols; + endif + elseif (nargout == 2) + varargout{1} = solution.t; ## Time stamps are first output argument + varargout{2} = solution.x; ## Results are second output argument + elseif (nargout == 5) + varargout{1} = solution.t; ## Same as (nargout == 2) + varargout{2} = solution.x; ## Same as (nargout == 2) + varargout{3} = []; ## LabMat doesn't accept lines like + varargout{4} = []; ## varargout{3} = varargout{4} = []; + varargout{5} = []; + if (vodeoptions.vhaveeventfunction) + varargout{3} = solution.vevent{3}; ## Time info when an event occured + varargout{4} = solution.vevent{4}; ## Results when an event occured + varargout{5} = solution.vevent{2}; ## Index info which event occured + endif + endif + +endfunction + +%! # We are using the "Van der Pol" implementation for all tests that +%! # are done for this function. +%! # For further tests we also define a reference solution (computed at high accuracy) +%!function [ydot] = fpol (vt, vy) ## The Van der Pol +%! ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; +%!endfunction +%!function [vref] = fref () ## The computed reference sol +%! vref = [0.32331666704577, -1.83297456798624]; +%!endfunction +%!function [vjac] = fjac (vt, vy, varargin) ## its Jacobian +%! vjac = [0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]; +%!function [vjac] = fjcc (vt, vy, varargin) ## sparse type +%! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]) +%!function [vval, vtrm, vdir] = feve (vt, vy, varargin) +%! vval = fpol (vt, vy, varargin); ## We use the derivatives +%! vtrm = zeros (2,1); ## that's why component 2 +%! vdir = ones (2,1); ## seems to not be exact +%!function [vval, vtrm, vdir] = fevn (vt, vy, varargin) +%! vval = fpol (vt, vy, varargin); ## We use the derivatives +%! vtrm = ones (2,1); ## that's why component 2 +%! vdir = ones (2,1); ## seems to not be exact +%!function [vmas] = fmas (vt, vy, varargin) +%! vmas = [1, 0; 0, 1]; ## Dummy mass matrix for tests +%!function [vmas] = fmsa (vt, vy, varargin) +%! vmas = sparse ([1, 0; 0, 1]); ## A sparse dummy matrix +%!function [vout] = fout (vt, vy, vflag, varargin) +%! if (regexp (char (vflag), 'init') == 1) +%! if (any (size (vt) ~= [2, 1])) error ('"fout" step "init"'); end +%! elseif (isempty (vflag)) +%! if (any (size (vt) ~= [1, 1])) error ('"fout" step "calc"'); end +%! vout = false; +%! elseif (regexp (char (vflag), 'done') == 1) +%! if (any (size (vt) ~= [1, 1])) error ('"fout" step "done"'); end +%! else error ('"fout" invalid vflag'); +%! end +%! +%! ## Turn off output of warning messages for all tests, turn them on +%! ## again if the last test is called +%!error ## ouput argument +%! warning ('off', 'OdePkg:InvalidArgument'); +%! B = ode45 (1, [0 25], [3 15 1]); +%!error ## input argument number one +%! [vt, vy] = ode45 (1, [0 25], [3 15 1]); +%!error ## input argument number two +%! [vt, vy] = ode45 (@fpol, 1, [3 15 1]); +%!test ## two output arguments +%! [vt, vy] = ode45 (@fpol, [0 2], [2 0]); +%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); +%!test ## not too many steps +%! [vt, vy] = ode45 (@fpol, [0 2], [2 0]); +%! assert (size (vt) < 20); +%!test ## anonymous function instead of real function +%! fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)]; +%! [vt, vy] = ode45 (fvdb, [0 2], [2 0]); +%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); +%!test ## extra input arguments passed through +%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], 12, 13, 'KL'); +%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); +%!test ## empty OdePkg structure *but* extra input arguments +%! vopt = odeset; +%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt, 12, 13, 'KL'); +%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); +%!error ## strange OdePkg structure +%! vopt = struct ('foo', 1); +%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt); +%!test ## Solve vdp in fixed step sizes +%! vopt = odeset('TimeStepSize', 0.1); +%! [vt, vy] = ode45 (@fpol, [0,2], [2 0], vopt); +%! assert (vt(:), [0:0.1:2]', 1e-2); +%!test ## Solve another anonymous function below zero +%! vref = [0, 14.77810590694212]; +%! [vt, vy] = ode45 (@(t,y) y, [-2 0], 2); +%! assert ([vt(end), vy(end,:)], vref, 1e-1); +%!test ## InitialStep option +%! vopt = odeset ('InitialStep', 1e-8); +%! [vt, vy] = ode45 (@fpol, [0 0.2], [2 0], vopt); +%! assert ([vt(2)-vt(1)], [1e-8], 1e-9); +%!test ## MaxStep option +%! vopt = odeset ('MaxStep', 1e-3); +%! vsol = ode45 (@fpol, [0 0.2], [2 0], vopt); +%! assert ([vsol.x(5)-vsol.x(4)], [1e-3], 1e-3); +%!test ## Solve with intermidiate step +%! vsol = ode45 (@fpol, [0 1 2], [2 0]); +%! assert (any((vsol.x-1) == 0)); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## Solve in backward direction starting at t=0 +%! vref = [-1.205364552835178, 0.951542399860817]; +%! vsol = ode45 (@fpol, [0 -2], [2 0]); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2); +%!test ## Solve in backward direction starting at t=2 +%! vref = [-1.205364552835178, 0.951542399860817]; +%! vsol = ode45 (@fpol, [2 -2], fref); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2); +%!test ## Solve in backward direction starting at t=2, with intermidiate step +%! vref = [-1.205364552835178, 0.951542399860817]; +%! vsol = ode45 (@fpol, [2 0 -2], fref); +%! idx = find(vsol.x < 0, 1, 'first') - 1; +%! assert ([vsol.x(idx), vsol.y(idx,:)], [0 2 0], 1e-2); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-2); +%!test ## Solve another anonymous function in backward direction +%! vref = [-1, 0.367879437558975]; +%! vsol = ode45 (@(t,y) y, [0 -1], 1); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); +%!test ## Solve another anonymous function below zero +%! vref = [0, 14.77810590694212]; +%! vsol = ode45 (@(t,y) y, [-2 0], 2); +%! assert ([vsol.x(end), vsol.y(end,:)], vref, 1e-3); +%!test ## Solve in backward direction starting at t=0 with MaxStep option +%! vref = [-1.205364552835178, 0.951542399860817]; +%! vopt = odeset ('MaxStep', 1e-3); +%! vsol = ode45 (@fpol, [0 -2], [2 0], vopt); +%! assert ([abs(vsol.x(8)-vsol.x(7))], [1e-3], 1e-3); +%! assert ([vsol.x(end), vsol.y(end,:)], [-2, vref], 1e-3); +%!test ## AbsTol option +%! vopt = odeset ('AbsTol', 1e-5); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## AbsTol and RelTol option +%! vopt = odeset ('AbsTol', 1e-8, 'RelTol', 1e-8); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## RelTol and NormControl option -- higher accuracy +%! vopt = odeset ('RelTol', 1e-8, 'NormControl', 'on'); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-5); +%!test ## Keeps initial values while integrating +%! vopt = odeset ('NonNegative', 2); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, 2, 0], 0.5); +%!test ## Details of OutputSel and Refine can't be tested +%! vopt = odeset ('OutputFcn', @fout, 'OutputSel', 1, 'Refine', 5); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%!test ## Details of OutputSave can't be tested +%! vopt = odeset ('OutputSave', 1, 'OutputSel', 1); +%! vsla = ode45 (@fpol, [0 2], [2 0], vopt); +%! vopt = odeset ('OutputSave', 2); +%! vslb = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert (length (vsla.x) + 1 >= 2 * length (vslb.x)) +%!test ## Stats must add further elements in vsol +%! vopt = odeset ('Stats', 'on'); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert (isfield (vsol, 'stats')); +%! assert (isfield (vsol.stats, 'nsteps')); +%!test ## Events option add further elements in vsol +%! vopt = odeset ('Events', @feve); +%! vsol = ode45 (@fpol, [0 10], [2 0], vopt); +%! assert (isfield (vsol, 'ie')); +%! assert (vsol.ie(1), 2); +%! assert (isfield (vsol, 'xe')); +%! assert (isfield (vsol, 'ye')); +%!test ## Events option, now stop integration +%! vopt = odeset ('Events', @fevn, 'NormControl', 'on'); +%! vsol = ode45 (@fpol, [0 10], [2 0], vopt); +%! assert ([vsol.ie, vsol.xe, vsol.ye], ... +%! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); +%!test ## Events option, five output arguments +%! vopt = odeset ('Events', @fevn, 'NormControl', 'on'); +%! [vt, vy, vxe, vye, vie] = ode45 (@fpol, [0 10], [2 0], vopt); +%! assert ([vie, vxe, vye], ... +%! [2.0, 2.496110, -0.830550, -2.677589], 6e-1); +%!test ## Jacobian option +%! vopt = odeset ('Jacobian', @fjac); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## Jacobian option and sparse return value +%! vopt = odeset ('Jacobian', @fjcc); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%! +%! ## test for JPattern option is missing +%! ## test for Vectorized option is missing +%! ## test for NewtonTol option is missing +%! ## test for MaxNewtonIterations option is missing +%! +%!test ## Mass option as function +%! vopt = odeset ('Mass', @fmas); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## Mass option as matrix +%! vopt = odeset ('Mass', eye (2,2)); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## Mass option as sparse matrix +%! vopt = odeset ('Mass', sparse (eye (2,2))); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## Mass option as function and sparse matrix +%! vopt = odeset ('Mass', @fmsa); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## Mass option as function and MStateDependence +%! vopt = odeset ('Mass', @fmas, 'MStateDependence', 'strong'); +%! vsol = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vsol.x(end), vsol.y(end,:)], [2, fref], 1e-3); +%!test ## Set BDF option to something else than default +%! vopt = odeset ('BDF', 'on'); +%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt); +%! assert ([vt(end), vy(end,:)], [2, fref], 1e-3); +%! +%! ## test for MvPattern option is missing +%! ## test for InitialSlope option is missing +%! ## test for MaxOrder option is missing +%! +%! warning ('on', 'OdePkg:InvalidArgument'); + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/odeget.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,174 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{}] =} odeget (@var{ode_opt}, @var{field}, [@var{default}], [@var{opt}]) +## @deftypefnx {Function File} {[@var{res}] =} odeget (@var{ode_opt}, @var{field}, [@var{default}], [@var{opt}]) +## +## If this function is called with two input arguments and the first input +## argument @var{ode_opt} is of type structure array and the second input +## argument @var{field} is of type string then return the option +## value @var{res} that is specified by the option name @var{field} in +## the ODE option structure @var{ode_opt}. Optionally if this function +## is called with a third input argument then return the default value +## @var{default} if @var{field} is not set in the structure @var{ode_opt}. +## @end deftypefn +## + +## Note: 20061022, Thomas Treichl +## We cannot create a function of the form odeget (@var{odestruct}, +## @var{name1}, @var{name2}) because we would get a mismatch with +## the function form 1 like described above. + +function res = odeget (ode_opt, field, default, opt) + + ## Check number and types of input arguments + if (nargin == 1) + error ("OdePkg:InvalidArgument", + "input arguments number must be at least 2") + endif + + if (isempty (ode_opt)) + if (nargin == 2) + res = []; + return + endif + res = default; + return + endif + + if (! isstruct (ode_opt)) + error ("OdePkg:InvalidArgument", + "first input argument must be a valid ODE_STRUCT."); + endif + + if (! ischar (field)) + error ("OdePkg:InvalidArgument", + "second input argument must be a string."); + endif + + if (nargin == 2) + default = []; + endif + + if ((nargin == 4) + && strcmp (lower (deblank (opt)), "fast")) + if (isfield (ode_opt, field)) + res = ode_opt.(field); + else + res = default; + endif + return + endif + + if ((nargin == 4) + && strcmp (lower (deblank (opt)), "fast_not_empty")) + if (isfield (ode_opt, field) + && ! isempty (ode_opt.(field)) ) + res = ode_opt.(field); + else + res = default; + endif + return + endif + + ## check if the given struct is a valid OdePkg struct + ode_struct_value_check (ode_opt); + + ## define all the possible OdePkg fields + options = ["AbsTol"; "Algorithm"; "BDF"; "Choice"; "Eta"; "Events"; + "Explicit"; "InexactSolver"; "InitialSlope"; "InitialStep"; + "Jacobian";"JConstant";"JPattern";"Mass"; "MassConstant"; + "MassSingular"; "MaxNewtonIterations"; "MaxOrder"; "MaxStep"; + "MStateDependence"; "MvPattern"; "NewtonTol"; "NonNegative"; + "NormControl"; "OutputFcn"; "OutputSave"; "OutputSel"; + "PolynomialDegree"; "QuadratureOrder"; "Refine"; "RelTol"; + "Restart"; "Stats"; "TimeStepNumber"; "TimeStepSize"; + "UseJacobian"; "Vectorized"]; + + while (1) + pos = fuzzy_compare (field, options); + + if (size (pos, 1) == 0) # no matching for the given option + if (nargin == 2) + error ("OdePkg:InvalidArgument", + "invalid property. No property found with name ''%s''.", field); + endif + warning ("odeget:NoExactMatching", + "no property found with name ''%s''. ", + "Assuming default value.", field); + res = default; + return + endif + + if (size (pos, 1) == 1) # one matching + if (! strcmp (lower (deblank (field)), + lower (deblank (options(pos,:)))) ) + warning ("OdePkg:InvalidArgument", + "no exact matching for ''%s''. ", + "Assuming you was intending ''%s''.", + field, deblank (options(pos,:))); + endif + res = ode_opt.(deblank (options(pos,:))); + if (isempty (res)) + res = default; + return + endif + return + endif + + ## if there are more matching, ask the user to be more precise + warning ("OdePkg:InvalidArgument", ... + "no exact matching for ''%s''. %d possible fields were found.", + field, size (pos, 1)); + for j = 1:(size (pos, 1)) + fprintf ("%s\n", deblank (options(pos(j),:))); + endfor + do + fprintf ("Please insert field name again.\n"); + field = input ("New field name: "); + until (ischar (field)) + endwhile + +endfunction + +%! ## Turn off output of warning messages for all tests, turn them on +%! ## again if the last test is called +%! warning ('off', 'OdePkg:InvalidArgument'); +%!test assert (odeget (odeset (), 'RelTol'), []); +%!test assert (odeget (odeset (), 'RelTol', 10), 10); +%!test assert (odeget (odeset (), 'Stats'), []); +%!test assert (odeget (odeset (), 'Stats', 'on'), 'on'); +%!test assert (odeget (odeset (), 'AbsTol', 1.e-6, 'fast'), []); +%!test assert (odeget (odeset (), 'AbsTol', 1.e-6, 'fast_not_empty'), 1.e-6); +%!test assert (odeget (odeset (), 'AbsTol', 1e-9), 1e-9); +%! +%! warning ('on', 'OdePkg:InvalidArgument'); + +%!demo +%! # Return the manually changed value RelTol of the OdePkg options +%! # strutcure A. If RelTol wouldn't have been changed then an +%! # empty matrix value would have been returned. +%! +%! A = odeset ('RelTol', 1e-1, 'AbsTol', 1e-2); +%! odeget (A, 'RelTol', []) + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/odeset.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,366 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{}] =} odeset () +## @deftypefnx {Function File} {[@var{odestruct}] =} odeset (@var{"field1"}, @var{value1}, @var{"field2"}, @var{value2}, @dots{}) +## @deftypefnx {Function File} {[@var{odestruct}] =} odeset (@var{oldstruct}, @var{"field1"}, @var{value1}, @var{"field2"}, @var{value2}, @dots{}) +## @deftypefnx {Function File} {[@var{odestruct}] =} odeset (@var{oldstruct}, +## @var{newstruct}) +## +## If this function is called without an input argument then return a new +## ODE options structure array that contains all the necessary fields and +## sets the values of all fields to default values. +## +## If this function is called with string input arguments @var{"field1"}, +## @var{"field2"}, @dots{} identifying valid ODE options then return a +## new ODE options structure with all necessary fields and set the values +## of the fields @var{"field1"}, @var{"field2"}, @dots{} to the values +## @var{value1}, @var{value2}, @dots{} +## +## If this function is called with a first input argument @var{oldstruct} +## of type structure array then overwrite all values of the options +## @var{"field1"}, @var{"field2"}, @dots{} of the structure @var{oldstruct} +## with new values @var{value1}, @var{value2}, @dots{} and return the +## modified structure array. +## +## If this function is called with two input arguments @var{oldstruct} +## and @var{newstruct} of type structure array then overwrite all values +## in the fields from the structure @var{oldstruct} with new values of the +## fields from the structure @var{newstruct}. Empty values of @var{newstruct} +## will not overwrite values in @var{oldstruct}. +## @end deftypefn +## + +function opt = odeset (varargin) + + ## Check number and types of all input arguments + if ((nargin == 0) + && (nargout == 0)) + print_options; + return + endif + + ## Creating a vector of OdePkg possible fields + fields = ["AbsTol"; "Algorithm"; "BDF"; "Choice"; "Eta"; "Events"; + "Explicit"; "InexactSolver"; "InitialSlope"; "InitialStep"; + "Jacobian";"JConstant";"JPattern";"Mass"; "MassConstant"; + "MassSingular"; "MaxNewtonIterations"; "MaxOrder"; "MaxStep"; + "MStateDependence"; "MvPattern"; "NewtonTol"; "NonNegative"; + "NormControl"; "OutputFcn"; "OutputSave"; "OutputSel"; + "PolynomialDegree"; "QuadratureOrder"; "Refine"; "RelTol"; + "Restart"; "Stats"; "TimeStepNumber"; "TimeStepSize"; + "UseJacobian"; "Vectorized"]; + + fields_nb = size (fields, 1); + + ## initialize output + opt = []; + for i = 1:1:fields_nb + opt.(deblank (fields(i,:))) = []; + endfor + + opt.Refine = 0; + opt.OutputSave = 1; + + if ((nargin == 0) + && (nargout == 1)) + return + endif + + ode_fields = fieldnames (opt); + + if (isstruct (varargin{1})) + ode_struct_value_check (varargin{1}); + + optA_fields = fieldnames (varargin{1}); + optA_f_nb = length (optA_fields); + + ## loop on first struct options for updating + for i = 1:optA_f_nb + name = lower (deblank (optA_fields{i})); + + while (1) + pos = fuzzy_compare (name, fields); + if (size (pos, 1) == 0) + warning ("OdePkg:InvalidArgument", + "no property found with name ''%s''", name); + endif + + if (size (pos, 1) == 1) + if (! strcmp (lower (deblank (name)), + lower (deblank (fields(pos,:))))) + warning ("OdePkg:InvalidArgument", "no exact matching for ", + "''%s''. Assuming you were intending ''%s''", + name, deblank (fields(pos,:))); + endif + + opt.(deblank (fields(pos,:))) = varargin{1}.(optA_fields{i}); + break + endif + + ## if there are more matching, ask the user to be more precise + warning ("OdePkg:InvalidArgument", + "no exact matching for ''%s''. %d possible fields were found", + name, size(pos, 1)); + for j = 1:(size (pos, 1)) + fprintf ("%s\n", deblank (fields(pos(j),:))); + endfor + do + fprintf ("Please insert field name again\n"); + name = input ("New field name: "); + until (ischar (name)) + endwhile + endfor + + if ((nargin == 2) + && (isstruct (varargin{2}))) + ode_struct_value_check (varargin{2}); + + optB_fields = fieldnames (varargin{2}); + optB_f_nb = length (optB_fields); + + ## update the first struct with the values in the second one + for i = 1:optB_f_nb + name = lower (deblank (optB_fields{i})); + while (1) + pos = fuzzy_compare (name, fields); + + if (size (pos, 1) == 0) + warning ("OdePkg:InvalidArgument", ... + "no property found with name ''%s''", name); + endif + + if (size(pos, 1) == 1) + if (! strcmp (lower (deblank (name)), lower (deblank (fields(pos,:))))) + warning ("OdePkg:InvalidArgument", "no exact matching for ", + "''%s''. Assuming you were intending ''%s''", + name, deblank (fields(pos,:))); + endif + opt.(deblank (fields(pos,:))) = varargin{2}.(optB_fields{i}); + break + endif + + ## if there are more matching, ask the user to be more precise + warning ("OdePkg:InvalidArgument", "no exact matching for ''%s''. ", + "%d possible fields were found", + name, size (pos, 1)); + for j = 1:(size (pos, 1)) + fprintf ("%s\n", deblank (fields(pos(j),:))); + endfor + do + fprintf ("Please insert field name again\n"); + name = input ("New field name: "); + until (ischar (name)) + endwhile + endfor + return + endif + + ## if the second argument is not a struct, + ## pass new values of the OdePkg options to the first struct + if (mod (nargin, 2) != 1) + error ("OdePkg:InvalidArgument", + "odeset expects an odd number of input arguments", + " when the first is a ODE_STRUCT"); + endif + + ## loop on the input arguments + for i = 2:2:(nargin - 1) + + if (! ischar(varargin{i})) + error ("OdePkg:InvalidArgument", + "not all odd input arguments are strings"); + endif + + name = varargin{i}; + + while (1) + pos = fuzzy_compare (name, fields); + + if (size (pos, 1) == 0) + error ("OdePkg:InvalidArgument", + "no property found with name ''%s''", name); + endif + + if (size (pos, 1) == 1) + if (! strcmp (lower (deblank (name)), lower (deblank (fields(pos,:))))) + warning ("OdePkg:InvalidArgument", "no exact matching for ''%s''. ", + "%d possible fields were found", + name, size (pos, 1)); + endif + opt.(deblank (fields(pos,:))) = varargin{i+1}; + break + endif + + ## if there are more matching, ask the user to be more precise + warning ("OdePkg:InvalidArgument", "no exact matching for ''%s''. ", + "%d possible fields were found", + name, size (pos, 1)); + for j = 1:(size (pos, 1)) + fprintf ("%s\n", deblank (fields(pos(j),:))); + endfor + do + fprintf ("Please insert field name again\n"); + name = input ("New field name: "); + until (ischar (name)) + endwhile + endfor + + ## check if all has been done gives a valid OdePkg struct + ode_struct_value_check (opt); + return + endif + + ## first input argument was not a struct + if (mod (nargin, 2) != 0) + error ("OdePkg:InvalidArgument", "odeset expects an even number ", + "of input arguments when the first is a string"); + endif + + for i = 1:2:(nargin-1) + if (! ischar (varargin{i})) + error ("OdePkg:InvalidArgument", + "not all even input arguments are strings"); + endif + + name = varargin{i}; + + while (1) + pos = fuzzy_compare (name, fields); + + if (size (pos, 1) == 0) + error ("OdePkg:InvalidArgument", + "invalid property. No property found with name ''%s''", name); + endif + + if (size (pos, 1) == 1) + if (! strcmp (lower (deblank (name)), + lower (deblank (fields(pos,:))))) + warning ("OdePkg:InvalidArgument", "no exact matching for ", + "''%s''. Assuming you were intending ''%s''", + name, deblank (fields(pos,:))); + endif + opt.(deblank (fields(pos,:))) = varargin{i+1}; + break + endif + + ## if there are more matching, ask the user to be more precise + warning ("OdePkg:InvalidArgument", "no exact matching for ''%s''. ", + "%d possible fields were found", + name, size (pos, 1)); + for j = 1:(size (pos, 1)) + fprintf ("%s\n", deblank (fields(pos(j),:))); + endfor + do + fprintf ("Please insert field name again\n"); + name = input ("New field name: "); + until (ischar (name)) + endwhile + endfor + + ## check if all has been done gives a valid OdePkg struct + ode_struct_value_check (opt); +endfunction + +## function useful to print all the possible options +function print_options () + + fprintf ("These following are all possible options.\n", + "Default values are put in square brackets.\n\n"); + + disp (" AbsTol: scalar or vector, >0, [1.e-6]"); + disp (" Algorithm: string, {[''gmres''], ''pcg'', ''bicgstab''}"); + disp (" BDF: binary, {''on'', [''off'']}"); + disp (" Choice: switch, {[1], 2}"); + disp (" Eta: scalar, >=0, <1, [0.5]"); + disp (" Events: function_handle, []"); + disp (" Explicit: binary, {''yes'', [''no'']}"); + disp (" InexactSolver: string, {''inexact_newton'', ''fsolve'', []}"); + disp (" InitialSlope: vector, []"); + disp (" InitialStep: scalar, >0, []"); + disp (" Jacobian: matrix or function_handle, []"); + disp (" JConstant: binary, {''on'', [''off'']}"); + disp (" JPattern: sparse matrix, []"); + disp (" Mass: matrix or function_handle, []"); + disp (" MassConstant: binary, {''on'', [''off'']}"); + disp (" MassSingular: switch, {''yes'', [''maybe''], ''no''}"); + disp ("MaxNewtonIterations: scalar, integer, >0, [1.e3]"); + disp (" MaxOrder: switch, {1, 2, 3, 4, [5]}"); + disp (" MaxStep: scalar, >0, []"); + disp (" MStateDependence: switch, {''none'', [''weak''], ''strong''}"); + disp (" MvPattern: sparse matrix, []"); + disp (" NewtonTol: scalar or vector, >0, []"); + disp (" NonNegative: vector of integers, []"); + disp (" NormControl: binary, {''on'', [''off'']}"); + disp (" OutputFcn: function_handle, []"); + disp (" OutputSave: scalar, integer, >0, []"); + disp (" OutputSel: scalar or vector, []"); + disp (" PolynomialDegree: scalar, integer, >0, []"); + disp (" QuadratureOrder: scalar, integer, >0, []"); + disp (" Refine: scalar, integer, >0, []"); + disp (" RelTol: scalar, >0, [1.e-3]"); + disp (" Restart: scalar, integer, >0, [20]"); + disp (" Stats: binary, {''on'', [''off'']}"); + disp (" TimeStepNumber: scalar, integer, >0, []"); + disp (" TimeStepSize: scalar, >0, []"); + disp (" UseJacobian: binary, {''yes'', [''no'']}"); + disp (" Vectorized: binary, {''on'', [''off'']}"); + +endfunction + +## All tests that are needed to check if a correct resp. valid option +## has been set are implemented in ode_struct_value_check.m. +%! ## Turn off output of warning messages for all tests, turn them on +%! ## again if the last test is called +%! warning ('off', 'OdePkg:InvalidArgument'); +%!test odeoptA = odeset (); +%!test odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1); +%! if (odeoptB.AbsTol != 1e-2), error; endif +%! if (odeoptB.RelTol != 1e-1), error; endif +%!test odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1); +%! odeoptC = odeset (odeoptB, 'NormControl', 'on'); +%!test odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1); +%! odeoptC = odeset (odeoptB, 'NormControl', 'on'); +%! odeoptD = odeset (odeoptC, odeoptB); +%! +%! warning ('on', 'OdePkg:InvalidArgument'); + +%!demo +%! # A new OdePkg options structure with default values is created. +%! +%! odeoptA = odeset (); +%! +%!demo +%! # A new OdePkg options structure with manually set options +%! # "AbsTol" and "RelTol" is created. +%! +%! odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1); +%! +%!demo +%! # A new OdePkg options structure from odeoptB is created with +%! # a modified value for option "NormControl". +%! +%! odeoptB = odeset ('AbsTol', 1e-2, 'RelTol', 1e-1); +%! odeoptC = odeset (odeoptB, 'NormControl', 'on'); + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/AbsRel_Norm.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,31 @@ +function res = AbsRel_Norm (x, x_old, AbsTol, RelTol, normcontrol, y) + + n = length (x); + + if (nargin == 5) + y = zeros (size (x)); + elseif (nargin != 6) + error ("OdePkg:InvalidArgument", + "invalid number of input arguments"); + endif + + if (length (x_old) != n + || length (y) != n) + error ("OdePkg:InvalidArgument", + "invalid dimensions of input arguments"); + endif + + if ((length (AbsTol) != 1 && length (AbsTol) != n) + || (length (RelTol) != 1 && length (RelTol) != n)) + error ("OdePkg:InvalidArgument", + "invalid dimensions of input arguments"); + endif + + sc = AbsTol + max (abs (x), abs (x_old)) .* RelTol; + if (normcontrol) + res = max (abs (x - y) ./ sc); + else + res = sqrt ((1 / n) * sum (((x - y) ./ sc).^2)); + endif + +endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/fuzzy_compare.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,172 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Command} {[@var{res}] =} fuzzy_compare @ +## (@var{"string1"}, @var{string_set}, [@var{correctness}]) +## +## Compare a string with a set of strings and returns the positions +## in the set of strings at which there are the fields that best fit +## the one we are comparing. +## +## The distance used to compare the words is the Levenshtein distance +## and for more details see +## @url{http://en.wikipedia.org/wiki/Levenshtein_distance}. +## +## This function must be called with one output argument @var{res} +## which contains the positions of the elements in @var{string_set} +## which best fit the given word. The tolerance that is used to +## determine if a field of the list fits or not the given word is a +## function of the length of the word and of the minimum distance of +## the word from all the elements of the list. The more the length, +## the more the tolerance. The less the minimum, the less the +## tolerance but if the minimum is close to the length of the word, +## the tolerance must be small because it means that no field in the +## list is really fitting the given word. So that the function is: +## +## @ifhtml +## @example +## @math{tolerance = 2 * (length-minimum) * minimum / length} +## @end example +## @end ifhtml +## @ifnothtml +## @math{tolerance = 2 * (length-minimum) * minimum / length}. +## @end ifnothtml +## +## The first input argument must be a string containing the word to +## compare. +## +## The second input argument must be a vector of strings or a +## cell_array of strings and should contain the fields to use for the +## comparison. +## +## The third input argument is optional and represents a fixed +## tolerance that will replace the implemented one. +## @end deftypefn +## +## @seealso{odeset,odeget,levenshtein} + +function res = fuzzy_compare (string1, string_set, correctness) + + ## check on output arguments + if (nargout > 1) + error ("OdePkg:InvalidArgument", "too many output arguments"); + endif + + ## check on input arguments + if (nargin < 2 || nargin > 3) + error ("OdePkg:InvalidArgument", "wrong input arguments number"); + endif + + if (! ischar (string1) + || (! iscellstr (string_set) + && ! ischar (string_set))) + error ("OdePkg:InvalidArgument", + "first argument must be a string, second argument ", + "must be an array of strings or a cell array of strings"); + endif + + if (nargin == 3) + if ((! isnumeric (correctness) || ! isscalar (correctness)) + && (! ischar (correctness))) + error ("OdePkg:InvalidArgument", + "third input argument must be a positive ", + "integer or a string"); + endif + + if (isnumeric (correctness) + && ( correctness < 0 || mod (correctness, 1) != 0)) + error ("OdePkg:InvalidArgument", + "third input argument must be a positive integer"); + endif + endif + + res = []; + + m = length (string1); + fields_nb = size (string_set, 1); + + values = inf .* ones (fields_nb, 1); + + string1 = deblank (string1); + string2 = []; + + minimus = inf; + ## loop on every field of the list + for i = 1:fields_nb + if (iscellstr (string_set)) + string2 = deblank (string_set{i}); + else + string2 = deblank (string_set(i,:)); + endif + ## compute Levenshtein distance (not case sensitive) + values(i) = levenshtein (lower (string1), + lower (string2), + minimus); + ## update the upper_bound to speedup the computation + minimus = min (minimus, values(i)); + endfor + + positions = find (values == minimus); + + if (minimus == 0) # exact match + if (size (positions, 1) != 1) + error ("OdePkg:InvalidArgument", + "there are %d strings perfectly matching ''%s''", + size (positions, 1), string1); + endif + res = positions; + return + endif + + ## determine the tolerance with the formula described in the + ## textinfo section it is a downwards parable with zeros in 0 and m + ## and with a maximum in m/2 of value m/2 + tolerance = m * (-(minimus - m) * minimus * (2 / (m*m))); + + ## if the degree of correctness is fixed by the user, it will + ## replace the tolerance + if (nargin == 3) + if ((isnumeric (correctness) + && isscalar (correctness) + && correctness == 0) + || (ischar (correctness) + && strcmp (lower (deblank (correctness)), 'exact'))) + error ("OdePkg:InvalidArgument", + "no exact matching for string ''%s''", string1); + endif + if (isnumeric (correctness) + && isscalar (correctness)) + tolerance = correctness; + endif + endif + + ## returning the positions of the fields whose distance is lower + ## than the tolerance + for i = 1:1:fields_nb + if (values(i) <= tolerance) + res = [res; i]; + endif + endfor + +endfunction + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/hermite_quartic_interpolation.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,68 @@ +## Copyright (C) 2014, Jacopo Corno <jacopo.corno@gmail.com> +## OdePkg - A package for solving ordinary differential equations and more +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Command} {[@var{x_out}] =} +## hermite_quartic_interpolation (@var{t}, @var{x}, @var{der}, @var{t_out}) +## +## This function file can be called by a ODE solver function in order to +## interpolate the solution at the time @var{t_out} using 4th order +## hermite interpolation. +## +## This function must be called with one output arguments: @var{x_out} +## which contains the evaluation at @var{t_out} of the hermite interpolant. +## +## The first input argument is the vector with two given times. +## +## The second input argument is the vector with the values of the function +## to interpolate at the times specified in @var{t} and at the middle point. +## +## The third input argument is the value of the derivatives of the function +## evaluated at the two extreme points. +## +## @end deftypefn +## +## @seealso{linear_interpolation, quadratic_interpolation, +## hermite_cubic_interpolation, hermite_quintic_interpolation, +## dorpri_interpolation} + +function x_out = hermite_quartic_interpolation (t, x, der, t_out) + + # Rescale time on [0,1] + s = (t_out - t(1)) / (t(2) - t(1)); + + # Hermite basis functions + # H0 = 1 - 11*s.^2 + 18*s.^3 - 8*s.^4; + # H1 = s - 4*s.^2 + 5*s.^3 - 2*s.^4; + # H2 = 16*s.^2 - 32*s.^3 + 16*s.^4; + # H3 = - 5*s.^2 + 14*s.^3 - 8*s.^4; + # H4 = s.^2 - 3*s.^3 + 2*s.^4; + + x_out = zeros (size (x, 1), length (t_out)); + for ii = 1:size (x, 1) + x_out(ii,:) = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4)*x(ii,1) ... + + ( s - 4*s.^2 + 5*s.^3 - 2*s.^4)*(t(2)-t(1))*der(ii,1) ... + + ( 16*s.^2 - 32*s.^3 + 16*s.^4)*x(ii,2) ... + + ( - 5*s.^2 + 14*s.^3 - 8*s.^4)*x(ii,3) ... + + ( s.^2 - 3*s.^3 + 2*s.^4)*(t(2)-t(1))*der(ii,2); + endfor + +endfunction + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/integrate_adaptive.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,387 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## OdePkg - A package for solving ordinary differential equations and more +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Command} {[@var{t}, @var{y}] =} +## integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@fun}, @var{tspan}, +## @var{x0}, @var{options}) +## +## This function file can be called by a ODE solver function in order to +## integrate the set of ODEs on the interval @var{[t0,t1]} with an +## adaptive timestep. +## +## This function must be called with two output arguments: @var{t} and @var{y}. +## Variable @var{t} is a column vector and contains the time stamps, instead +## @var{y} is a matrix in which each column refers to a different unknown +## of the problem and the rows number is the same of @var{t} rows number so +## that each row of @var{y} contains the values of all unknowns at the time +## value contained in the corresponding row in @var{t}. +## +## The first input argument must be a function_handle or an inline function +## representing the stepper, that is the function responsible for step-by-step +## integration. This function discriminates one method from the others. +## +## The second input argument is the order of the stepper. It is needed +## to compute the adaptive timesteps. +## +## The third input argument is a function_handle or an inline function that +## defines the set of ODE: +## @ifhtml +## @example +## @math{y' = f(t,y)} +## @end example +## @end ifhtml +## @ifnothtml +## @math{y' = f(t,y)}. +## @end ifnothtml +## +## The fourth input argument is the time vector which defines integration +## interval, that is @var{[tspan(1),tspan(end)]} and all the intermediate +## elements are taken as times at which the solution is required. +## +## The fifth argument represents the initial conditions for the ODEs and the +## last input argument contains some options that may be needed for the stepper. +## +## @end deftypefn +## +## @seealso{integrate_const, integrate_n_steps} + +function solution = integrate_adaptive (stepper, order, func, tspan, x0, options) + + solution = struct; + + ## first values for time and solution + t = tspan(1); + x = x0(:); + + ## get first initial timestep + dt = odeget (options, "InitialStep", + starting_stepsize (order, func, t, x, options.AbsTol, + options.RelTol, options.vnormcontrol), + "fast_not_empty"); + vdirection = odeget (options, "vdirection", [], "fast"); + if (sign (dt) != vdirection) + dt = -dt; + endif + dt = vdirection * min (abs (dt), options.MaxStep); + + ## set parameters + k = length (tspan); + counter = 2; + comp = 0.0; + tk = tspan(1); + options.comp = comp; + + ## factor multiplying the stepsize guess + facmin = 0.8; + fac = 0.38^(1/(order+1)); ## formula taken from Hairer + t_caught = false; + + + ## Initialize the OutputFcn + if (options.vhaveoutputfunction) + if (options.vhaveoutputselection) + solution.vretout = x(options.OutputSel,end); + else + solution.vretout = x; + endif + feval (options.OutputFcn, tspan, solution.vretout, + "init", options.vfunarguments{:}); + endif + + ## Initialize the EventFcn + if (options.vhaveeventfunction) + odepkg_event_handle (options.Events, t(end), x, + "init", options.vfunarguments{:}); + endif + + solution.vcntloop = 2; + solution.vcntcycles = 1; + vcntiter = 0; + solution.vunhandledtermination = true; + solution.vcntsave = 2; + + z = t; + u = x; + + k_vals = feval (func, t , x, options.vfunarguments{:}); + + while (counter <= k) + facmax = 1.5; + + ## compute integration step from t to t+dt + [s, y, y_est, k_vals] = stepper (func, z(end), u(:,end), + dt, options, k_vals); + + if (options.vhavenonnegative) + x(options.NonNegative,end) = abs (x(options.NonNegative,end)); + y(options.NonNegative,end) = abs (y(options.NonNegative,end)); + y_est(options.NonNegative,end) = abs (y_est(options.NonNegative,end)); + endif + + if (options.vhaveoutputfunction && options.vhaverefine) + vSaveVUForRefine = u(:,end); + endif + + err = AbsRel_Norm (y(:,end), u(:,end), options.AbsTol, options.RelTol, + options.vnormcontrol, y_est(:,end)); + + ## solution accepted only if the error is less or equal to 1.0 + if (err <= 1) + + [tk, comp] = kahan (tk, comp, dt); + options.comp = comp; + s(end) = tk; + + ## values on this interval for time and solution + z = [z(end);s]; + u = [u(:,end),y]; + + ## if next tspan value is caught, update counter + if ((z(end) == tspan(counter)) + || (abs (z(end) - tspan(counter)) / + (max (abs (z(end)), abs (tspan(counter)))) < 8*eps) ) + counter++; + + ## if there is an element in time vector at which the solution is required + ## the program must compute this solution before going on with next steps + elseif (vdirection * z(end) > vdirection * tspan(counter)) + ## initialize counter for the following cycle + i = 2; + while (i <= length (z)) + + ## if next tspan value is caught, update counter + if ((counter <= k) + && ((z(i) == tspan(counter)) + || (abs (z(i) - tspan(counter)) / + (max (abs (z(i)), abs (tspan(counter)))) < 8*eps)) ) + counter++; + endif + ## else, loop until there are requested values inside this subinterval + while ((counter <= k) + && (vdirection * z(i) > vdirection * tspan(counter))) + ## choose interpolation scheme according to order of the solver + switch order + case 1 + u_interp = linear_interpolation ([z(i-1) z(i)], + [u(:,i-1) u(:,i)], + tspan(counter)); + case 2 + if (! isempty (k_vals)) + der = k_vals(:,1); + else + der = feval (func, z(i-1) , u(:,i-1), + options.vfunarguments{:}); + endif + u_interp = quadratic_interpolation ([z(i-1) z(i)], + [u(:,i-1) u(:,i)], + der, tspan(counter)); + case 3 + u_interp = ... + hermite_cubic_interpolation ([z(i-1) z(i)], + [u(:,i-1) u(:,i)], + [k_vals(:,1) k_vals(:,end)], + tspan(counter)); + case 4 + ## if ode45 is used without local extrapolation this function + ## doesn't require a new function evaluation. + u_interp = dorpri_interpolation ([z(i-1) z(i)], + [u(:,i-1) u(:,i)], + k_vals, tspan(counter)); + case 5 + ## ode45 with Dormand-Prince scheme: + ## 4th order approximation of y in t+dt/2 as proposed by + ## Shampine in Lawrence, Shampine, "Some Practical + ## Runge-Kutta Formulas", 1986. + u_half = u(:,i-1) ... + + 1/2*dt*((6025192743/30085553152) * k_vals(:,1) + + (51252292925/65400821598) * k_vals(:,3) + - (2691868925/45128329728) * k_vals(:,4) + + (187940372067/1594534317056) * k_vals(:,5) + - (1776094331/19743644256) * k_vals(:,6) + + (11237099/235043384) * k_vals(:,7)); + u_interp = ... + hermite_quartic_interpolation ([z(i-1) z(i)], + [u(:,i-1) u_half u(:,i)], + [k_vals(:,1) k_vals(:,end)], + tspan(counter)); + + ## it is also possible to do a new function evaluation and use + ## the quintic hermite interpolator + ## f_half = feval (func, t+1/2*dt, u_half, + ## options.vfunarguments{:}); + ## u_interp = + ## hermite_quintic_interpolation ([z(i-1) z(i)], + ## [u(:,i-1) u_half u(:,i)], + ## [k_vals(:,1) f_half k_vals(:,end)], + ## tspan(counter)); + otherwise + warning ("High order interpolation not yet implemented: ", + "using cubic iterpolation instead"); + der(:,1) = feval (func, z(i-1) , u(:,i-1), + options.vfunarguments{:}); + der(:,2) = feval (func, z(i) , u(:,i), + options.vfunarguments{:}); + u_interp = ... + hermite_cubic_interpolation ([z(i-1) z(i)], + [u(:,i-1) u(:,i)], + der, tspan(counter)); + endswitch + + ## add the interpolated value of the solution + u = [u(:,1:i-1), u_interp, u(:,i:end)]; + + ## add the time requested + z = [z(1:i-1);tspan(counter);z(i:end)]; + + ## update counters + counter++; + i++; + endwhile + + ## if new time requested is not out of this interval + if ((counter <= k) + && (vdirection * z(end) > vdirection * tspan(counter))) + ## update the counter + i++; + else + ## stop the cycle and go on with the next iteration + i = length (z) + 1; + endif + + endwhile + endif + + if (mod (solution.vcntloop-1, options.OutputSave) == 0) + x = [x,u(:,2:end)]; + t = [t;z(2:end)]; + solution.vcntsave = solution.vcntsave + 1; + endif + solution.vcntloop = solution.vcntloop + 1; + vcntiter = 0; + + ## Call plot only if a valid result has been found, therefore this + ## code fragment has moved here. Stop integration if plot function + ## returns false + if (options.vhaveoutputfunction) + for vcnt = 0:options.Refine # Approximation between told and t + if (options.vhaverefine) # Do interpolation + vapproxtime = (vcnt + 1) / (options.Refine + 2); + vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ... + + (vapproxtime) * y(:,end); + vapproxtime = s(end) + vapproxtime * dt; + else + vapproxvals = x(:,end); + vapproxtime = t(end); + endif + if (options.vhaveoutputselection) + vapproxvals = vapproxvals(options.OutputSel); + endif + vpltret = feval (options.OutputFcn, vapproxtime, + vapproxvals, [], options.vfunarguments{:}); + if (vpltret) # Leave refinement loop + break + endif + endfor + if (vpltret) # Leave main loop + solution.vunhandledtermination = false; + break + endif + endif + + ## Call event only if a valid result has been found, therefore this + ## code fragment has moved here. Stop integration if veventbreak is + ## true + if (options.vhaveeventfunction) + solution.vevent = odepkg_event_handle (options.Events, t(end), + x(:,end), [], options.vfunarguments{:}); + if (! isempty (solution.vevent{1}) + && solution.vevent{1} == 1) + t(solution.vcntloop-1,:) = solution.vevent{3}(end,:); + x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)'; + solution.vunhandledtermination = false; + break + endif + endif + + else + facmax = 1.0; + endif + + ## Compute next timestep, formula taken from Hairer + err += eps; # adding an eps to avoid divisions by zero + dt = dt * min (facmax, max (facmin, + fac * (1 / err)^(1 / (order + 1)))); + dt = vdirection * min (abs (dt), options.MaxStep); + + ## Update counters that count the number of iteration cycles + solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics + vcntiter = vcntiter + 1; # Needed to find iteration problems + + ## Stop solving because in the last 1000 steps no successful valid + ## value has been found + if (vcntiter >= 5000) + error (["Solving has not been successful. The iterative", + " integration loop exited at time t = %f before endpoint at", + " tend = %f was reached. This happened because the iterative", + " integration loop does not find a valid solution at this time", + " stamp. Try to reduce the value of ''InitialStep'' and/or", + " ''MaxStep'' with the command ''odeset''.\n"], + s(end), tspan(end)); + endif + + ## if this is the last iteration, save the length of last interval + if (counter > k) + j = length (z); + endif + endwhile + + ## Check if integration of the ode has been successful + if (vdirection * z(end) < vdirection * tspan(end)) + if (solution.vunhandledtermination == true) + error ("OdePkg:InvalidArgument", + ["Solving has not been successful. The iterative", + " integration loop exited at time t = %f", + " before endpoint at tend = %f was reached. This may", + " happen if the stepsize grows smaller than defined in", + " vminstepsize. Try to reduce the value of ''InitialStep''", + " and/or ''MaxStep'' with the command ''odeset''.\n"], + z(end), tspan(end)); + else + warning ("OdePkg:InvalidArgument", + ["Solver has been stopped by a call of ''break'' in the main", + " iteration loop at time t = %f before endpoint at tend = %f ", + " was reached. This may happen because the @odeplot function", + " returned ''true'' or the @event function returned", + " ''true''.\n"], + z(end), tspan(end)); + endif + endif + + ## Compute how many values are out of time inerval + d = vdirection * t((end-(j-1)):end) > vdirection * tspan(end)*ones (j, 1); + f = sum (d); + + ## Remove not-requested values of time and solution + solution.t = t(1:end-f); + solution.x = x(:,1:end-f)'; + +endfunction + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/integrate_const.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,289 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## OdePkg - A package for solving ordinary differential equations and more +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Command} {[@var{t}, @var{y}] =} integrate_const (@var{@@stepper}, +## @var{@@fun}, @var{tspan}, @var{x0}, @var{dt}, @var{options}) +## +## This function file can be called by an ODE solver function in order to +## integrate the set of ODEs on the interval @var{[t0,t1]} with a +## constant timestep @var{dt}. +## +## This function must be called with two output arguments: @var{t} and @var{y}. +## Variable @var{t} is a column vector and contains the time stamps, +## instead @var{y} is a matrix in which each column refers to a different +## unknown of the problem and the rows number is the same of @var{t} rows +## number so that each row of @var{y} contains the values of all unknowns at +## the time value contained in the corresponding row in @var{t}. +## +## The first input argument must be a function_handle or an inline function +## representing the stepper, that is the function responsible for step-by-step +## integration. This function discriminates one method from the others. +## +## The second input argument is the order of the stepper. It is needed to +## compute the adaptive timesteps. +## +## The third input argument is a function_handle or an inline function +## that defines the set of ODE: +## +## @ifhtml +## @example +## @math{y' = f(t,y)} +## @end example +## @end ifhtml +## @ifnothtml +## @math{y' = f(t,y)}. +## @end ifnothtml +## +## The third input argument is the time vector which defines integration +## interval, that is @var{[tspan(1),tspan(end)]} and all the intermediate +## elements are taken as times at which the solution is required. +## +## The fourth argument contains the initial conditions for the ODEs. +## +## The fifth input argument represents the fixed timestep and the last input +## argument contains some options that may be needed for the stepper. +## @end deftypefn +## +## @seealso{integrate_adaptive, integrate_n_steps} + +function solution = integrate_const (stepper, func, tspan, x0, dt, options) + + solution = struct; + + ## first values for time and solution + t = tspan(1); + x = x0(:); + + vdirection = odeget (options, "vdirection", [], "fast"); + if (sign (dt) != vdirection) + error ("OdePkg:InvalidArgument", + "option ''InitialStep'' has a wrong sign"); + endif + + ## setting parameters + k = length (tspan); + counter = 2; + comp = 0.0; + tk = tspan(1); + options.comp = comp; + + ## Initialize the OutputFcn + if (options.vhaveoutputfunction) + if (options.vhaveoutputselection) + solution.vretout = x(options.OutputSel,end); + else + solution.vretout = x; + endif + feval (options.OutputFcn, tspan, solution.vretout, "init", + options.vfunarguments{:}); + endif + + ## Initialize the EventFcn + if (options.vhaveeventfunction) + odepkg_event_handle (options.Events, t(end), x, "init", + options.vfunarguments{:}); + endif + + solution.vcntloop = 2; + solution.vcntcycles = 1; + #vu = vinit; + #vk = vu.' * zeros(1,6); + vcntiter = 0; + solution.vunhandledtermination = true; + solution.vcntsave = 2; + + z = t; + u = x; + k_vals = feval (func, t , x, options.vfunarguments{:}); + + while (counter <= k) + ## computing the integration step from t to t+dt + [s, y, ~, k_vals] = stepper (func, z(end), u(:,end), dt, options, k_vals); + + [tk, comp] = kahan (tk,comp, dt); + options.comp = comp; + s(end) = tk; + + if (options.vhavenonnegative) + x(options.NonNegative,end) = abs (x(options.NonNegative,end)); + y(options.NonNegative,end) = abs (y(options.NonNegative,end)); + y_est(options.NonNegative,end) = abs (y_est(options.NonNegative,end)); + endif + + if (options.vhaveoutputfunction && options.vhaverefine) + vSaveVUForRefine = u(:,end); + endif + + ## values on this interval for time and solution + z = [t(end);s]; + u = [x(:,end),y]; + + ## if next tspan value is caught, update counter + if ((z(end) == tspan(counter)) + || (abs (z(end) - tspan(counter)) / + (max(abs (z(end)), abs(tspan(counter)))) < 8*eps) ) + counter++; + + ## if there is an element in time vector at which the solution is required + ## the program must compute this solution before going on with next steps + elseif (vdirection * z(end) > vdirection * tspan(counter) ) + ## initializing counter for the following cycle + i = 2; + while (i <= length (z)) + + ## if next tspan value is caught, update counter + if ((counter <= k) + && (((z(i) == tspan(counter)) + || (abs (z(i) - tspan(counter)) / + (max(abs (z(i)), abs (tspan(counter)))) < 8*eps))) ) + counter++; + endif + ## else, loop until there are requested values inside this subinterval + while ((counter <= k) + && vdirection * z(i) > vdirection * tspan(counter) ) + ## add the interpolated value of the solution + u = [u(:,1:i-1),u(:,i-1) + (tspan(counter)-z(i-1))/(z(i)-z(i-1))* ... + (u(:,i)-u(:,i-1)),u(:,i:end)]; + ## add the time requested + z = [z(1:i-1);tspan(counter);z(i:end)]; + + ## update counters + counter++; + i++; + endwhile + + ## if new time requested is not out of this interval + if ((counter <= k) + && vdirection * z(end) > vdirection * tspan(counter)) + ## update the counter + i++; + else + ## else, stop the cycle and go on with the next iteration + i = length (z)+1; + endif + + endwhile + endif + + if (mod (solution.vcntloop-1, options.OutputSave) == 0) + x = [x,u(:,2:end)]; + t = [t;z(2:end)]; + solution.vcntsave = solution.vcntsave + 1; + endif + solution.vcntloop = solution.vcntloop + 1; + vcntiter = 0; + + ## Call plot only if a valid result has been found, therefore this + ## code fragment has moved here. Stop integration if plot function + ## returns false + if (options.vhaveoutputfunction) + for vcnt = 0:options.Refine # Approximation between told and t + if (options.vhaverefine) # Do interpolation + vapproxtime = (vcnt + 1) / (options.Refine + 2); + vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ... + + (vapproxtime) * y(:,end); + vapproxtime = s(end) + vapproxtime*dt; + else + vapproxvals = x(:,end); + vapproxtime = t(end); + endif + if (options.vhaveoutputselection) + vapproxvals = vapproxvals(options.OutputSel); + endif + vpltret = feval (options.OutputFcn, vapproxtime, vapproxvals, [], + options.vfunarguments{:}); + if (vpltret) # Leave refinement loop + break + endif + endfor + if (vpltret) # Leave main loop + solution.vunhandledtermination = false; + break + endif + endif + + ## Call event only if a valid result has been found, therefore this + ## code fragment has moved here. Stop integration if veventbreak is true + if (options.vhaveeventfunction) + solution.vevent = odepkg_event_handle (options.Events, t(end), x(:,end), + [], options.vfunarguments{:}); + if (! isempty (solution.vevent{1}) + && solution.vevent{1} == 1) + t(solution.vcntloop-1,:) = solution.vevent{3}(end,:); + x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)'; + solution.vunhandledtermination = false; + break + endif + endif + + ## Update counters that count the number of iteration cycles + solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics + vcntiter = vcntiter + 1; # Needed to find iteration problems + + ## Stop solving because the last 1000 steps no successful valid + ## value has been found + if (vcntiter >= 5000) + error (["Solving has not been successful. The iterative", + " integration loop exited at time t = %f before endpoint at", + " tend = %f was reached. This happened because the iterative", + " integration loop does not find a valid solution at this time", + " stamp. Try to reduce the value of ''InitialStep'' and/or", + " ''MaxStep'' with the command ''odeset''.\n"], + s(end), tspan(end)); + endif + + ## if this is the last iteration, save the length of last interval + if (counter > k) + j = length (z); + endif + endwhile + + ## Check if integration of the ode has been successful + if (vdirection * z(end) < vdirection * tspan(end)) + if (solution.vunhandledtermination == true) + error ("OdePkg:InvalidArgument", + ["Solving has not been successful. The iterative integration loop", + " exited at time t = %f before endpoint at tend = %f was", + " reached. This may happen if the stepsize grows smaller than", + " defined in vminstepsize. Try to reduce the value of", + " ''InitialStep'' and/or ''MaxStep'' with the command", + " ''odeset''.\n"], z(end), tspan(end)); + else + warning ("OdePkg:InvalidArgument", + ["Solver has been stopped by a call of ''break'' in the main", + " iteration loop at time t = %f before endpoint at tend = %f", + " was reached. This may happen because the @odeplot function", + " returned ''true'' or the @event function returned", + " ''true''.\n"], + z(end), tspan(end)); + endif + endif + + ## compute how many values are out of time inerval + d = vdirection * t((end-(j-1)):end) > vdirection * tspan(end) * ones (j, 1); + f = sum (d); + + ## remove not-requested values of time and solution + solution.t = t(1:end-f); + solution.x = x(:,1:end-f)'; + +endfunction + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/integrate_n_steps.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,231 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## OdePkg - A package for solving ordinary differential equations and more +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Command} {[@var{t}, @var{y}] =} integrate_n_steps +## (@var{@@stepper}, @var{@@fun}, @var{t0}, @var{x0}, @var{dt}, @var{n}, +## @var{options}) +## +## This function file can be called by an ODE solver function in order to +## integrate the set of ODEs on the interval @var{[t0,t0 + n*dt]} with a +## constant timestep dt and on a fixed number of steps. +## +## This function must be called with two output arguments: @var{t} and @var{y}. +## Variable @var{t} is a column vector and contains the time stamps, instead +## @var{y} is a matrix in which each column refers to a different unknown of +## the problem and the rows number is the same of @var{t} rows number so that +## each row of @var{y} contains the values of all unknowns at the time +## value contained in the corresponding row in @var{t}. +## +## The first input argument must be a function_handle or an inline function +## representing the stepper, that is the function responsible for step-by-step +## integration. This function discriminates one method from the others. +## +## The second input argument is the order of the stepper. It is needed +## to compute the adaptive timesteps. +## +## The third input argument is a function_handle or an inline function +## that defines the set of ODE: +## +## @ifhtml +## @example +## @math{y' = f(t,y)} +## @end example +## @end ifhtml +## @ifnothtml +## @math{y' = f(t,y)}. +## @end ifnothtml +## +## The third input argument is the starting point for the integration. +## +## The fourth argument contains the initial conditions for the ODEs. +## +## The fifth input argument represents the fixed timestep while the sixth +## contains the number of integration steps. +## +## The last argument is a struct with the options that may be +## needed by the stepper. +## @end deftypefn +## +## @seealso{integrate_adaptive, integrate_const} + +function solution = integrate_n_steps (stepper, func, t0, x0, dt, n, options) + + solution = struct; + + ## first values for time and solution + x = x0(:); + t = t0; + + vdirection = odeget (options, "vdirection", [], "fast"); + if (sign (dt) != vdirection) + error("OdePkg:InvalidArgument", + "option ''InitialStep'' has a wrong sign"); + endif + + comp = 0.0; + tk = t0; + options.comp = comp; + + ## Initialize the OutputFcn + if (options.vhaveoutputfunction) + if (options.vhaveoutputselection) + solution.vretout = x(options.OutputSel,end); + else + solution.vretout = x; + endif + feval (options.OutputFcn, tspan, solution.vretout, "init", + options.vfunarguments{:}); + endif + + ## Initialize the EventFcn + if (options.vhaveeventfunction) + odepkg_event_handle (options.Events, t(end), x, "init", + options.vfunarguments{:}); + endif + + solution.vcntloop = 2; + solution.vcntcycles = 1; + #vu = vinit; + #vk = vu.' * zeros(1,6); + vcntiter = 0; + solution.vunhandledtermination = true; + solution.vcntsave = 2; + + z = t; + u = x; + k_vals = feval (func, t , x, options.vfunarguments{:}); + + for i = 1:n + ## Compute the integration step from t to t+dt + [s, y, ~, k_vals] = stepper (func, z(end), u(:,end), dt, options, k_vals); + + [tk, comp] = kahan (tk, comp, dt); + options.comp = comp; + s(end) = tk; + + if (options.vhavenonnegative) + x(options.NonNegative,end) = abs (x(options.NonNegative,end)); + y(options.NonNegative,end) = abs (y(options.NonNegative,end)); + endif + + if (options.vhaveoutputfunction && options.vhaverefine) + vSaveVUForRefine = u(:,end); + endif + + ## values on this interval for time and solution + z = [t(end);s]; + u = [x(:,end),y]; + + if (mod (solution.vcntloop-1, options.OutputSave) == 0) + x = [x,u(:,2:end)]; + t = [t;z(2:end)]; + solution.vcntsave = solution.vcntsave + 1; + endif + solution.vcntloop = solution.vcntloop + 1; + vcntiter = 0; + + ## Call plot only if a valid result has been found, therefore this code + ## fragment has moved here. Stop integration if plot function returns false + if (options.vhaveoutputfunction) + for vcnt = 0:options.Refine # Approximation between told and t + if (options.vhaverefine) # Do interpolation + vapproxtime = (vcnt + 1) / (options.Refine + 2); + vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ... + + (vapproxtime) * y(:,end); + vapproxtime = s(end) + vapproxtime*dt; + else + vapproxvals = x(:,end); + vapproxtime = t(end); + endif + if (options.vhaveoutputselection) + vapproxvals = vapproxvals(options.OutputSel); + endif + vpltret = feval (options.OutputFcn, vapproxtime, vapproxvals, [], + options.vfunarguments{:}); + if (vpltret) # Leave refinement loop + break + endif + endfor + if (vpltret) # Leave main loop + solution.vunhandledtermination = false; + break + endif + endif + + ## Call event only if a valid result has been found, therefore this + ## code fragment has moved here. Stop integration if veventbreak is + ## true + if (options.vhaveeventfunction) + solution.vevent = odepkg_event_handle (options.Events, t(end), x(:,end), + [], options.vfunarguments{:}); + if (! isempty (solution.vevent{1}) + && solution.vevent{1} == 1) + t(solution.vcntloop-1,:) = solution.vevent{3}(end,:); + x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)'; + solution.vunhandledtermination = false; + break + endif + endif + + ## Update counters that count the number of iteration cycles + solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics + vcntiter = vcntiter + 1; # Needed to find iteration problems + + ## Stop solving because the last 1000 steps no successful valid + ## value has been found + if (vcntiter >= 5000) + error (["Solving has not been successful. The iterative", + " integration loop exited at time t = %f before endpoint at", + " tend = %f was reached. This happened because the iterative", + " integration loop does not find a valid solution at this time", + " stamp. Try to reduce the value of ''InitialStep'' and/or", + " ''MaxStep'' with the command ''odeset''.\n"], + s(end), tspan(end)); + endif + endfor + + ## Check if integration of the ode has been successful + #if (vdirection * z(end) < vdirection * tspan(end)) + # if (solution.vunhandledtermination == true) + # error ("OdePkg:InvalidArgument", + # ["Solving has not been successful. The iterative", + # " integration loop exited at time t = %f", + # " before endpoint at tend = %f was reached. This may", + # " happen if the stepsize grows smaller than defined in", + # " vminstepsize. Try to reduce the value of ''InitialStep''", + # " and/or ''MaxStep'' with the command ''odeset''.\n"], + # z(end), tspan(end)); + # else + # warning ("OdePkg:InvalidArgument", + # ["Solver has been stopped by a call of ''break'' in the main", + # " iteration loop at time t = %f before endpoint at tend = %f", + # " was reached. This may happen because the @odeplot function", + # " returned ''true'' or the @event function returned", + # " ''true''.\n"], + # z(end), tspan(end)); + # endif + #endif + + solution.t = t; + solution.x = x'; + +endfunction + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/kahan.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,55 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## OdePkg - A package for solving ordinary differential equations and more +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Command} {[@var{sum}] =} kahan (@var{sum}, +## @var{comp}, @var{temp}) +## @deftypefnx {Command} {[@var{sum}, @var{comp}] =} kahan (@var{sum}, +## @var{comp}, @var{temp}) +## +## This function is the implementation of the Kahan summation algorithm, +## also known as compensated summation. It significantly reduces the numerical +## error in the total obtained by adding a sequence of finite precision +## floating point numbers, compared to the obvious approach. For more details +## see @url{http://en.wikipedia.org/wiki/Kahan_summation_algorithm}. +## This function is called in @command{integrate_adaptive} and in +## @command{integrate_const} to better catch equality comparisons. +## +## The first input argument is the variable that will contain the summation, +## so that is also returned as first output argument in order to reuse it in +## next calls to kahan function. +## +## The second input argument contains the compensation term and it is returned +## as second output argument so that it can be reused in the next calls +## of the same computation. +## +## The third input argument is the variable that contains the term to +## be added to @var{Sum}. +## @end deftypefn + +function [Sum, comp] = kahan (Sum, comp, term) + + y = term - comp; + t = Sum + y; + comp = (t - Sum) - y; + Sum = t; + +endfunction + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/ode_struct_value_check.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,495 @@ +## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## OdePkg - A package for solving ordinary differential equations and more +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{}] =} +## ode_struct_value_check (@var{arg}, [@var{"solver"}]) +## +## If this function is called with one input argument of type structure array +## then check the field names and the field values of the OdePkg structure +## @var{arg}. Optionally if this function is called with a second input +## argument @var{"solver"} of type string that specifies the name of a valid +## OdePkg solver then a higher level error detection is performed. The function +## does not modify any of the field names or field values but terminates with +## an error if an invalid option or value is found. +## +## This function is an OdePkg internal helper function therefore it should +## never be necessary that this function is called directly by a user. +## @end deftypefn +## +## @seealso{odeset, odeget} + +function ode_struct_value_check (arg, solver) + + ## Check the number of input arguments + if (nargin == 0 || nargin > 2) + error ("OdePkg:InvalidArgument", + "wrong input arguments number"); + endif + + if (! isstruct (arg)) + error ("OdePkg:InvalidArgument", + "first input argument is not a struct"); + endif + + if (nargin == 1) + solver = []; + elseif (! ischar (solver) ) + error ("OdePkg:InvalidArgument", + "second input argument is not a string"); + endif + + fields = fieldnames (arg); + fields_nb = length (fields); + + for i = 1:1:fields_nb # Run through the number of given structure field names + switch (fields{i}) + + case "AbsTol" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i})) + || any (arg.(fields{i}) <= 0) + || ~isreal (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (! isvector (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif ( any (arg.(fields{i}) <= 0)) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Algorithm" + if (! isempty (arg.(fields{i}))) + if (! ischar (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "BDF" + if (! isempty (arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "on") + && ! strcmp (arg.(fields{i}), "off")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Choice" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (arg.(fields{i}) != 1 && arg.(fields{i}) != 2) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Eta" + if ( ~isempty (arg.(fields{i})) ) + if ( ~isreal (arg.(fields{i})) ) + error ("OdePkg:InvalidArgument", ... + "value assigned to field %s is not a valid one", fields{i}); + elseif ( arg.(fields{i})<0 || arg.(fields{i})>=1 ) + error ("OdePkg:InvalidArgument", ... + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Events" + if (! isempty (arg.(fields{i}))) + if (! isa (arg.(fields{i}), "function_handle")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Explicit" + if (! isempty (arg.(fields{i}))) + if (! ischar (arg.(fields{i})) + || (! strcmp (arg.(fields{i}), "yes") + && ! strcmp (arg.(fields{i}), "no"))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "InexactSolver" + if (! isempty (arg.(fields{i}))) + if (! ischar (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "InitialSlope" + if (! isempty (arg.(fields{i}))) + if (! ischar (arg.(fields{i})) + && (! isnumeric (arg.(fields{i})) + || (! isvector (arg.(fields{i})) + && ! isreal (arg.(fields{i}))))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "InitialStep" + if (! isempty (arg.(fields{i})) ) + if (! isnumeric (arg.(fields{i})) + || ! isscalar (arg.(fields{i})) + || ! isreal (arg.(fields{i})) + || arg.(fields{i}) <=0) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Jacobian" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i}))) + if (! isa (arg.(fields{i}), "function_handle") + && ! iscell (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", + fields{i}); + endif + endif + endif + + case "JConstant" + if (! isempty(arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "on") + && ! strcmp (arg.(fields{i}), "off")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "JPattern" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i})) + && ! isscalar (arg.(fields{i})) + && ! isvector (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Mass" + if (! isempty (arg.(fields{i}))) + if ((! isnumeric (arg.(fields{i})) + || ~ismatrix (arg.(fields{i}))) + && ! isa (arg.(fields{i}), "function_handle")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "MassConstant" + if (! isempty (arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "on") + && ! strcmp (arg.(fields{i}), "off")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "MassSingular" + if (! isempty (arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "yes") + && ! strcmp (arg.(fields{i}), "no") + && ! strcmp (arg.(fields{i}), "maybe")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "MaxNewtonIterations" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (mod (arg.(fields{i}), 1) != 0 + || arg.(fields{i}) <= 0) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "MaxOrder" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (mod (arg.(fields{i}), 1) != 0 + || arg.(fields{i}) <= 0 + || arg.(fields{i}) >= 8 ) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "MaxStep" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i})) + || ! isscalar (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (arg.(fields{i}) <= 0) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "MStateDependence" + if (! isempty (arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "none") + && ! strcmp (arg.(fields{i}), "weak") + && ! strcmp (arg.(fields{i}), "strong")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "MvPattern" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i})) + && ! isvector (arg.(fields{i})) ) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "NewtonTol" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i})) + || ! isreal (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (any (arg.(fields{i}) <= 0)) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "NonNegative" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i})) + || ! isvector (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (any (arg.(fields{i}) <= 0) + || any (mod (arg.(fields{i}), 1) != 0)) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "NormControl" + if (! isempty (arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "on") + && ! strcmp (arg.(fields{i}), "off")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "OutputFcn" + if (! isempty (arg.(fields{i}))) + if (! isa (arg.(fields{i}), "function_handle")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "OutputSave" + if (! isempty (arg.(fields{i}))) + if (! isscalar (arg.(fields{i})) + && arg.(fields{i}) != Inf) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif ((mod (arg.(fields{i}), 1) != 0 || arg.(fields{i}) <= 0) + && arg.(fields{i}) != Inf) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "OutputSel" + if (! isempty (arg.(fields{i}))) + if (! isscalar (arg.(fields{i})) ) + if (! isnumeric (arg.(fields{i})) + || ! isvector (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", + fields{i}); + endif + endif + endif + + case "PolynomialDegree" + if (! isempty (arg.(fields{i})) ) + if (! isnumeric (arg.(fields{i})) + || ! isvector (arg.(fields{i})) ) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (any (arg.(fields{i}) <= 0)) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "QuadratureOrder" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i})) + || ! isvector (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (any(arg.(fields{i}) <= 0)) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Refine" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i})) + || ! isscalar (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (mod (arg.(fields{i}), 1) != 0 + || arg.(fields{i}) < 0 + || arg.(fields{i}) > 5 ) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "RelTol" + if (! isempty (arg.(fields{i})) ) + if (! isnumeric (arg.(fields{i})) ) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (! isreal (arg.(fields{i})) + || any (arg.(fields{i}) <= 0)) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + switch (solver) + case {"ode23", "ode23d", "ode45", "ode45d", + "ode54", "ode54d", "ode78", "ode78d"} + if (! isempty (arg.(fields{i}))) + if (! isscalar (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "for this type of solver, value assigned to field %s ", + "is not a valid one", fields{i}); + endif + endif + otherwise + endswitch + + case "Restart" + if (! isempty (arg.(fields{i}))) + if (! isnumeric (arg.(fields{i}))) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + elseif (mod (arg.(fields{i}), 1) != 0 + || arg.(fields{i}) <=0 ) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Stats" + if (! isempty (arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "on") + && ! strcmp (arg.(fields{i}), "off")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "TimeStepNumber" + if (! isempty (arg.(fields{i}))) + if (mod (arg.(fields{i}), 1) != 0 + || arg.(fields{i}) <= 0) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "TimeStepSize" + if (! isempty (arg.(fields{i}))) + if (! isreal (arg.(fields{i})) + || (arg.(fields{i}) == 0)) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "UseJacobian" + if (! isempty (arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "yes") + && ! strcmp (arg.(fields{i}), "no")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + case "Vectorized" + if (! isempty (arg.(fields{i}))) + if (! strcmp (arg.(fields{i}), "on") + && ! strcmp (arg.(fields{i}), "off")) + error ("OdePkg:InvalidArgument", + "value assigned to field %s is not a valid one", fields{i}); + endif + endif + + otherwise + warning ("OdePkg:InvalidArgument", + "no fields with name %s in ODE options.", fields{i}); + endswitch + endfor + +endfunction + +%!demo +%! # Return the checked OdePkg options structure that is created by +%! # the command odeset. +%! +%! ode_struct_value_check (odeset); +%! +%!demo +%! # Create the OdePkg options structure A with odeset and check it +%! # with odepkg_structure_check. This actually is unnecessary +%! # because odeset automtically calls odepkg_structure_check before +%! # returning. +%! +%! A = odeset (); ode_struct_value_check (A); + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/odepkg_event_handle.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,136 @@ +%# Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +%# OdePkg - A package for solving ordinary differential equations and more +%# +%# This program is free software; you can redistribute it and/or modify +%# it under the terms of the GNU General Public License as published by +%# the Free Software Foundation; either version 2 of the License, or +%# (at your option) any later version. +%# +%# This program is distributed in the hope that it will be useful, +%# but WITHOUT ANY WARRANTY; without even the implied warranty of +%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +%# GNU General Public License for more details. +%# +%# You should have received a copy of the GNU General Public License +%# along with this program; If not, see <http://www.gnu.org/licenses/>. + +%# -*- texinfo -*- +%# @deftypefn {Function File} {[@var{sol}] =} odepkg_event_handle (@var{@@fun}, @var{time}, @var{y}, @var{flag}, [@var{par1}, @var{par2}, @dots{}]) +%# +%# Return the solution of the event function that is specified as the first input argument @var{@@fun} in form of a function handle. The second input argument @var{time} is of type double scalar and specifies the time of the event evaluation, the third input argument @var{y} either is of type double column vector (for ODEs and DAEs) and specifies the solutions or is of type cell array (for IDEs and DDEs) and specifies the derivatives or the history values, the third input argument @var{flag} is of type string and can be of the form +%# @table @option +%# @item @code{"init"} +%# then initialize internal persistent variables of the function @command{odepkg_event_handle} and return an empty cell array of size 4, +%# @item @code{"calc"} +%# then do the evaluation of the event function and return the solution @var{sol} as type cell array of size 4, +%# @item @code{"done"} +%# then cleanup internal variables of the function @command{odepkg_event_handle} and return an empty cell array of size 4. +%# @end table +%# Optionally if further input arguments @var{par1}, @var{par2}, @dots{} of any type are given then pass these parameters through @command{odepkg_event_handle} to the event function. +%# +%# This function is an OdePkg internal helper function therefore it should never be necessary that this function is called directly by a user. There is only little error detection implemented in this function file to achieve the highest performance. +%# @end deftypefn +%# +%# @seealso{odepkg} + +function [vretval] = odepkg_event_handle (vevefun, vt, vy, vflag, varargin) + + %# No error handling has been implemented in this function to achieve + %# the highest performance available. + + %# vretval{1} is true or false; either to terminate or to continue + %# vretval{2} is the index information for which event occured + %# vretval{3} is the time information column vector + %# vretval{4} is the line by line result information matrix + + %# These persistent variables are needed to store the results and the + %# time value from the processing in the time stamp before, veveold + %# are the results from the event function, vtold the time stamp, + %# vretcell the return values cell array, vyold the result of the ode + %# and vevecnt the counter for how often this event handling + %# has been called + persistent veveold; persistent vtold; + persistent vretcell; persistent vyold; + persistent vevecnt; + + %# Call the event function if an event function has been defined to + %# initialize the internal variables of the event function an to get + %# a value for veveold + if (strcmp (vflag, 'init')) + + if (~iscell (vy)) + vinpargs = {vevefun, vt, vy}; + else + vinpargs = {vevefun, vt, vy{1}, vy{2}}; + vy = vy{1}; %# Delete cell element 2 + end + if (nargin > 4) + vinpargs = {vinpargs{:}, varargin{:}}; + end + [veveold, vterm, vdir] = feval (vinpargs{:}); + + %# We assume that all return values must be column vectors + veveold = veveold(:)'; vterm = vterm(:)'; vdir = vdir(:)'; + vtold = vt; vyold = vy; vevecnt = 1; vretcell = cell (1,4); + + %# Process the event, find the zero crossings either for a rising + %# or for a falling edge + elseif (isempty (vflag)) + + if (~iscell (vy)) + vinpargs = {vevefun, vt, vy}; + else + vinpargs = {vevefun, vt, vy{1}, vy{2}}; + vy = vy{1}; %# Delete cell element 2 + end + if (nargin > 4) + vinpargs = {vinpargs{:}, varargin{:}}; + end + [veve, vterm, vdir] = feval (vinpargs{:}); + + %# We assume that all return values must be column vectors + veve = veve(:)'; vterm = vterm(:)'; vdir = vdir(:)'; + + %# Check if one or more signs of the event has changed + vsignum = (sign (veveold) ~= sign (veve)); + if (any (vsignum)) %# One or more values have changed + vindex = find (vsignum); %# Get the index of the changed values + + if (any (vdir(vindex) == 0)) + %# Rising or falling (both are possible) + %# Don't change anything, keep the index + elseif (any (vdir(vindex) == sign (veve(vindex)))) + %# Detected rising or falling, need a new index + vindex = find (vdir == sign (veve)); + else + %# Found a zero crossing but must not be notified + vindex = []; + end + + %# Create new output values if a valid index has been found + if (~isempty (vindex)) + %# Change the persistent result cell array + vretcell{1} = any (vterm(vindex)); %# Stop integration or not + vretcell{2}(vevecnt,1) = vindex(1,1); %# Take first event found + %# Calculate the time stamp when the event function returned 0 and + %# calculate new values for the integration results, we do both by + %# a linear interpolation + vtnew = vt - veve(1,vindex) * (vt - vtold) / (veve(1,vindex) - veveold(1,vindex)); + vynew = (vy - (vt - vtnew) * (vy - vyold) / (vt - vtold))'; + vretcell{3}(vevecnt,1) = vtnew; + vretcell{4}(vevecnt,:) = vynew; + vevecnt = vevecnt + 1; + end %# if (~isempty (vindex)) + + end %# Check for one or more signs ... + veveold = veve; vtold = vt; vretval = vretcell; vyold = vy; + + elseif (strcmp (vflag, 'done')) %# Clear this event handling function + clear ('veveold', 'vtold', 'vretcell', 'vyold', 'vevecnt'); + vretcell = cell (1,4); + + end + +%# Local Variables: *** +%# mode: octave *** +%# End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/odepkg_structure_check.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,455 @@ +%# Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +%# Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +%# OdePkg - A package for solving ordinary differential equations and more +%# +%# This program is free software; you can redistribute it and/or modify +%# it under the terms of the GNU General Public License as published by +%# the Free Software Foundation; either version 2 of the License, or +%# (at your option) any later version. +%# +%# This program is distributed in the hope that it will be useful, +%# but WITHOUT ANY WARRANTY; without even the implied warranty of +%# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +%# GNU General Public License for more details. +%# +%# You should have received a copy of the GNU General Public License +%# along with this program; If not, see <http://www.gnu.org/licenses/>. + +%# -*- texinfo -*- +%# @deftypefn {Function File} {[@var{newstruct}] =} odepkg_structure_check (@var{oldstruct}, [@var{"solver"}]) +%# +%# If this function is called with one input argument of type structure array then check the field names and the field values of the OdePkg structure @var{oldstruct} and return the structure as @var{newstruct} if no error is found. Optionally if this function is called with a second input argument @var{"solver"} of type string taht specifies the name of a valid OdePkg solver then a higher level error detection is performed. The function does not modify any of the field names or field values but terminates with an error if an invalid option or value is found. +%# +%# This function is an OdePkg internal helper function therefore it should never be necessary that this function is called directly by a user. There is only little error detection implemented in this function file to achieve the highest performance. +%# +%# Run examples with the command +%# @example +%# demo odepkg_structure_check +%# @end example +%# @end deftypefn +%# +%# @seealso{odepkg} + +function [vret] = odepkg_structure_check (varargin) + + %# Check the number of input arguments + if (nargin == 0) + help ('odepkg_structure_check'); + error ('OdePkg:InvalidArgument', ... + 'Number of input arguments must be greater than zero'); + elseif (nargin > 2) + print_usage; + elseif (nargin == 1 && isstruct (varargin{1})) + vret = varargin{1}; + vsol = ''; + vfld = fieldnames (vret); + vlen = length (vfld); + elseif (nargin == 2 && isstruct (varargin{1}) && ischar (varargin{2})) + vret = varargin{1}; + vsol = varargin{2}; + vfld = fieldnames (vret); + vlen = length (vfld); + end + + for vcntarg = 1:vlen %# Run through the number of given structure field names + + switch (vfld{vcntarg}) + + case 'RelTol' + if (isempty (vret.(vfld{vcntarg})) || ... + (isnumeric (vret.(vfld{vcntarg})) && ... + isreal (vret.(vfld{vcntarg})) && ... + all (vret.(vfld{vcntarg}) > 0))) %# 'all' is a MatLab need + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + switch (vsol) + case {'ode23', 'ode45', 'ode54', 'ode78', ... + 'ode23d', 'ode45d', 'ode54d', 'ode78d',} + if (~isscalar (vret.(vfld{vcntarg})) && ... + ~isempty (vret.(vfld{vcntarg}))) + error ('OdePkg:InvalidParameter', ... + 'Value of option "RelTol" must be a scalar'); + end + otherwise + + end + + case 'AbsTol' + if (isempty (vret.(vfld{vcntarg})) || ... + (isnumeric (vret.(vfld{vcntarg})) && ... + isreal (vret.(vfld{vcntarg})) && ... + all (vret.(vfld{vcntarg}) > 0))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'NormControl' + if (isempty (vret.(vfld{vcntarg})) || ... + (strcmp (vret.(vfld{vcntarg}), 'on') || ... + strcmp (vret.(vfld{vcntarg}), 'off'))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'NonNegative' + if (isempty (vret.(vfld{vcntarg})) || ... + (isnumeric (vret.(vfld{vcntarg})) && isvector (vret.(vfld{vcntarg})))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'OutputFcn' + if (isempty (vret.(vfld{vcntarg})) || ... + isa (vret.(vfld{vcntarg}), 'function_handle')) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'OutputSel' + if (isempty (vret.(vfld{vcntarg})) || ... + (isnumeric (vret.(vfld{vcntarg})) && isvector (vret.(vfld{vcntarg}))) || ... + isscalar (vret.(vfld{vcntarg}))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'OutputSave' + if (isempty (vret.(vfld{vcntarg})) || ... + (isscalar (vret.(vfld{vcntarg})) && ... + mod (vret.(vfld{vcntarg}), 1) == 0 && ... + vret.(vfld{vcntarg}) > 0) || ... + vret.(vfld{vcntarg}) == Inf) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Refine' + if (isempty (vret.(vfld{vcntarg})) || ... + (isscalar (vret.(vfld{vcntarg})) && ... + mod (vret.(vfld{vcntarg}), 1) == 0 && ... + vret.(vfld{vcntarg}) >= 0 && ... + vret.(vfld{vcntarg}) <= 5)) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Stats' + if (isempty (vret.(vfld{vcntarg})) || ... + (strcmp (vret.(vfld{vcntarg}), 'on') || ... + strcmp (vret.(vfld{vcntarg}), 'off'))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'InitialStep' + if (isempty (vret.(vfld{vcntarg})) || ... + (isscalar (vret.(vfld{vcntarg})) && ... + isreal (vret.(vfld{vcntarg})))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'MaxStep' + if (isempty (vret.(vfld{vcntarg})) || ... + (isscalar (vret.(vfld{vcntarg})) && ... + vret.(vfld{vcntarg}) > 0) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Events' + if (isempty (vret.(vfld{vcntarg})) || ... + isa (vret.(vfld{vcntarg}), 'function_handle')) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Jacobian' + if (isempty (vret.(vfld{vcntarg})) || ... + isnumeric (vret.(vfld{vcntarg})) || ... + isa (vret.(vfld{vcntarg}), 'function_handle') || ... + iscell (vret.(vfld{vcntarg}))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'JPattern' + if (isempty (vret.(vfld{vcntarg})) || ... + isvector (vret.(vfld{vcntarg})) || ... + isnumeric (vret.(vfld{vcntarg}))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Vectorized' + if (isempty (vret.(vfld{vcntarg})) || ... + (strcmp (vret.(vfld{vcntarg}), 'on') || ... + strcmp (vret.(vfld{vcntarg}), 'off'))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Mass' + if (isempty (vret.(vfld{vcntarg})) || ... + (isnumeric (vret.(vfld{vcntarg})) || ... + isa (vret.(vfld{vcntarg}), 'function_handle'))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'MStateDependence' + if (isempty (vret.(vfld{vcntarg})) || ... + (strcmp (vret.(vfld{vcntarg}), 'none') || ... + strcmp (vret.(vfld{vcntarg}), 'weak') || ... + strcmp (vret.(vfld{vcntarg}), 'strong'))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'MvPattern' + if (isempty (vret.(vfld{vcntarg})) || ... + (isvector (vret.(vfld{vcntarg})) || ... + isnumeric (vret.(vfld{vcntarg})))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'MassSingular' + if (isempty (vret.(vfld{vcntarg})) || ... + (strcmp (vret.(vfld{vcntarg}), 'yes') || ... + strcmp (vret.(vfld{vcntarg}), 'no') || ... + strcmp (vret.(vfld{vcntarg}), 'maybe'))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'InitialSlope' + if (isempty (vret.(vfld{vcntarg})) || ... + isvector (vret.(vfld{vcntarg}))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'MaxOrder' + if (isempty (vret.(vfld{vcntarg})) || ... + (mod (vret.(vfld{vcntarg}), 1) == 0 && ... + vret.(vfld{vcntarg}) > 0 && ... + vret.(vfld{vcntarg}) < 8)) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'BDF' + if (isempty (vret.(vfld{vcntarg})) || ... + (strcmp (vret.(vfld{vcntarg}), 'on') || ... + strcmp (vret.(vfld{vcntarg}), 'off'))) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'NewtonTol' + if (isempty (vret.(vfld{vcntarg})) || ... + (isnumeric (vret.(vfld{vcntarg})) && ... + isreal (vret.(vfld{vcntarg})) && ... + all (vret.(vfld{vcntarg}) > 0))) %# 'all' is a MatLab need + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'MaxNewtonIterations' + if (isempty (vret.(vfld{vcntarg})) || ... + (mod (vret.(vfld{vcntarg}), 1) == 0 && ... + vret.(vfld{vcntarg}) > 0)) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + %# new fields added + case 'Algorithm' + if( isempty(vret.(vfld{vcntarg})) || ischar(vret.(vfld{vcntarg})) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Choice' + if( isempty(vret.(vfld{vcntarg})) || (isnumeric(vret.(vfld{vcntarg})) && (vret.(vfld{vcntarg})==1) || vret.(vfld{vcntarg})==2 ) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Eta' + if( isempty(vret.(vfld{vcntarg})) || ( isreal(vret.(vfld{vcntarg})) && vret.(vfld{vcntarg})>=0 && vret.(vfld{vcntarg})<1) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Explicit' + if( isempty(vret.(vfld{vcntarg})) || (ischar(vret.(vfld{vcntarg})) && (strcmp(vret.(vfld{vcntarg}),'yes') || strcmp(vret.(vfld{vcntarg}),'no'))) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'InexactSolver' + if( isempty(vret.(vfld{vcntarg})) || ischar(vret.(vfld{vcntarg})) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'InitialSlope' + if( isempty(vret.(vfld{vcntarg})) || ( ischar(vret.(vfld{vcntarg})) || (isnumeric(vret.(vfld{vcntarg})) && (isvector(vret.(vfld{vcntarg})) || isreal(vret.(vfld{vcntarg}))))) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'JConstant' + if( isempty(vret.(vfld{vcntarg})) || (ischar(vret.(vfld{vcntarg})) && (strcmp(vret.(vfld{vcntarg}),'yes') || strcmp(vret.(vfld{vcntarg}),'no'))) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'MassConstant' + if( isempty(vret.(vfld{vcntarg})) || (strcmp(vret.(vfld{vcntarg}),'on') || strcmp(vret.(vfld{vcntarg}),'off')) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'PolynomialDegree' + if( isempty(vret.(vfld{vcntarg})) || (isnumeric(vret.(vfld{vcntarg})) && mod(vret.(vfld{vcntarg}),1)==0 && vret.(vfld{vcntarg})>0) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'QuadratureOrder' + if( isempty(vret.(vfld{vcntarg})) || (isnumeric(vret.(vfld{vcntarg})) && mod(vret.(vfld{vcntarg}),1)==0 && vret.(vfld{vcntarg})>0) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'Restart' + if( isempty(vret.(vfld{vcntarg})) || (isnumeric(vret.(vfld{vcntarg})) && mod(vret.(vfld{vcntarg}),1)==0 && vret.(vfld{vcntarg})>0) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'TimeStepNumber' + if( isempty(vret.(vfld{vcntarg})) || (isnumeric(vret.(vfld{vcntarg})) && mod(vret.(vfld{vcntarg}),1)==0 && vret.(vfld{vcntarg})>0) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'TimeStepSize' + if( isempty(vret.(vfld{vcntarg})) || ( isreal(vret.(vfld{vcntarg})) && vret.(vfld{vcntarg})!=0) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + case 'UseJacobian' + if( isempty(vret.(vfld{vcntarg})) || (ischar(vret.(vfld{vcntarg})) && (strcmp(vret.(vfld{vcntarg}),'yes') || strcmp(vret.(vfld{vcntarg}),'no'))) ) + else + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s" or no valid parameter value', ... + vfld{vcntarg}); + end + + otherwise + error ('OdePkg:InvalidParameter', ... + 'Unknown parameter name "%s"', ... + vfld{vcntarg}); + + end %# switch + + end %# for + +%!demo +%! # Return the checked OdePkg options structure that is created by +%! # the command odeset. +%! +%! odepkg_structure_check (odeset); +%! +%!demo +%! # Create the OdePkg options structure A with odeset and check it +%! # with odepkg_structure_check. This actually is unnecessary +%! # because odeset automtically calls odepkg_structure_check before +%! # returning. +%! +%! A = odeset (); odepkg_structure_check (A); + +%# Local Variables: *** +%# mode: octave *** +%# End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/runge_kutta_45_dorpri.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,117 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## OdePkg - A package for solving ordinary differential equations and more +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {Command} {[@var{t_next}, @var{x_next}] =} +## runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, +## [@var{options}, @var{k_vals_in}]) +## @deftypefnx {Command} {[@var{t_next}, @var{x_next}, @var{x_est}] =} +## runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, +## [@var{options}, @var{k_vals_in}]) +## @deftypefnx {Command} {[@var{t_next}, @var{x_next}, @var{x_est}, +## @var{k_vals_out}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, +## @var{dt}, [@var{options}, @var{k_vals_in}]) +## +## This function can be used to integrate a system of ODEs with a given initial +## condition @var{x} from @var{t} to @var{t+dt}, with the Dormand-Prince method. +## For the definition of this method see +## @url{http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method}. +## +## First output argument is the final integration time value. +## +## Second output parameter is the higher order computed solution at time +## @var{t_next} (local extrapolation). +## +## Third output parameter is a lower order solution for the estimation +## of the error. +## +## Fourth output parameter is matrix containing the Runge-Kutta evaluations +## to use in a FSAL scheme or for dense output. +## +## First input argument is the function describing the system of ODEs +## to be integrated. +## +## Second input parameter is the first extreme of integration interval. +## +## Third input argument is the initial condition of the system. +## +## Fourth input argument is the timestep, that is the length of the +## integration interval. +## +## Fifth input parameter is optional and describes a set of options useful +## to adapt the computation to what is needed. +## +## Sixth input parameter is optional and describes the Runge-Kutta evaluations +## of the previous step to use in a FSAL scheme. +## @end deftypefn +## +## @seealso{odepkg} + +function varargout = runge_kutta_45_dorpri (f, t, x, dt, varargin) + + options = varargin{1}; + k = zeros (size (x, 1), 4); + + if (nargin == 5) # only the options are passed + k(:,1) = feval (f, t , x, options.vfunarguments{:}); + elseif (nargin == 6) # both the options and the k values are passed + k(:,1) = varargin{2}(:,end); # FSAL property + endif + k(:,1) = feval (f, t, x, options.vfunarguments{:}); + k(:,2) = feval (f, t + (1/5)*dt, ... + x + dt * (1/5)*k(:,1), ... + options.vfunarguments{:}); + k(:,3) = feval (f, t + (3/10)*dt, ... + x + dt * ((3/40)*k(:,1) + (9/40)*k(:,2)), ... + options.vfunarguments{:}); + k(:,4) = feval (f, t + (4/5)*dt, ... + x + dt * ((44/45)*k(:,1) - (56/15)*k(:,2) + (32/9)*k(:,3)), ... + options.vfunarguments{:}); + k(:,5) = feval (f, t + (8/9)*dt, ... + x + dt * ((19372/6561)*k(:,1) - (25360/2187)*k(:,2) ... + + (64448/6561)*k(:,3) - (212/729)*k(:,4)), ... + options.vfunarguments{:}); + k(:,6) = feval (f, t + dt, ... + x + dt * ((9017/3168)*k(:,1) - (355/33)*k(:,2) ... + + (46732/5247)*k(:,3) + (49/176)*k(:,4) ... + - (5103/18656)*k(:,5)), ... + options.vfunarguments{:}); + + ## compute new time and new values for the unkwnowns + varargout{1} = t + dt; + varargout{2} = x + dt * ((35/384)*k(:,1) + (500/1113)*k(:,3) + + (125/192)*k(:,4) - (2187/6784)*k(:,5) + + (11/84)*k(:,6)); # 5th order approximation + + ## if the estimation of the error is required + if (nargout >= 3) + ## new solution to be compared with the previous one + k(:,7) = feval (f, t + dt, varargout{2}, options.vfunarguments{:}); + ## x_est according to Shampine 1986: + ## varargout{3} = x + dt * ((1951/21600)*k(:,1) + (22642/50085)*k(:,3) + ## + (451/720)*k(:,4) - (12231/42400)*k(:,5) + ## + (649/6300)*k(:,6) + (1/60)*k(:,7)); + varargout{3} = x + dt * ((5179/57600)*k(:,1) + (7571/16695)*k(:,3) + + (393/640)*k(:,4) - (92097/339200)*k(:,5) + + (187/2100)*k(:,6) + (1/40)*k(:,7)); # x_est + varargout{4} = k; + endif + +endfunction + +## Local Variables: *** +## mode: octave *** +## End: ***
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/starting_stepsize.m Thu Sep 24 12:58:46 2015 +0200 @@ -0,0 +1,74 @@ +## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## OdePkg - A package for solving ordinary differential equations and more +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; If not, see <http://www.gnu.org/licenses/>. + + +## -*- texinfo -*- +## @deftypefn {Command} {[@var{h}] =} starting_stepsize (@var{order}, +## @var{@@fun}, @var{t0}, @var{x0}) +## +## This function file can be used to determine a good initial step for an ODE +## solver of order @var{order}. The algorithm is that one described in [1]. +## +## Second input argument, which is @var{@@fun}, is the function describing +## the differential equations, @var{t0} is the initial time and @var{x0} +## is the initial condition. +## +## This function returns a good guess for the initial timestep @var{h}. +## +## References: +## [1] E. Hairer, S.P. Norsett and G. Wanner, +## "Solving Ordinary Differential Equations I: Nonstiff Problems", Springer. +## @end deftypefn +## +## @seealso{odepkg} + +function h = starting_stepsize (order, func, t0, x0, + AbsTol, RelTol, normcontrol) + + ## compute norm of initial conditions + d0 = AbsRel_Norm (x0, x0, AbsTol, RelTol, normcontrol); + + ## compute norm of the function evaluated at initial conditions + y = func (t0, x0); + d1 = AbsRel_Norm (y, y, AbsTol, RelTol, normcontrol); + + if (d0 < 1.e-5 || d1 < 1.e-5) + h0 = 1.e-6; + else + h0 = .01 * (d0 / d1); + endif + + ## compute one step of Explicit-Euler + x1 = x0 + h0 * y; + + ## approximate the derivative norm + d2 = (1 / h0) * ... + AbsRel_Norm (func (t0+h0, x1) - y, + func (t0+h0, x1) - y, AbsTol, RelTol, normcontrol); + + if (max(d1, d2) <= 1.e-15) + h1 = max (1.e-6, h0*1.e-3); + else + h1 = (1.e-2 / max (d1, d2)) ^(1 / (order+1)); + endif + + h = min (100*h0, h1); + +endfunction + +## Local Variables: *** +## mode: octave *** +## End: ***