Mercurial > octave
changeset 20634:80e630b37ba1
maint: Remove unnecessary 'v' prefix before variables in ODE m-files.
* ode_rk_interpolate.m: Deleted file.
* odepkg_event_handle.m: Deleted file.
* runge_kutta_interpolate.m: Renamed from ode_rk_interpolate.m. Remove 'v'
prefix on variables. Delete blank space at end of lines.
* ode_event_handler.m: Renamed from odepkg_event_handle.m. Remove 'v'
prefix on variables. Delete blank space at end of lines. Use 'evt' for
event rather than 'eve' in variable names. Use 'idx' rather than 'index'
in variable names.
* scripts/ode/module.mk: Add ode_event_handler.m and runge_kutta_interpolate.m
to build system.
* AbsRel_Norm.m, starting_stepsize.m, ode_struct_value_check.m, odeget.m,
odeset.m: Delete blank space at end of lines.
* integrate_adaptive, integrate_const.m, integrate_n_steps.m,
runge_kutta_45_dorpri.m:
Remove 'v' prefix on variable. Delete blank space at end of lines.
* ode45.m: Expand docstring to cover more of the inputs/outputs.
Remove 'v' prefix on variable. Use name of variable in input validation
warnings. Use name of function as prefix in warnings and error messages.
Delete long, unnecessary comments. Use faster 'isempty' rather than slow
'isequal' to check whether option has been set. Remove SubOpts variable.
Shorten lines < 80 chars.
| author | Rik <rik@octave.org> |
|---|---|
| date | Sun, 18 Oct 2015 09:55:41 -0700 |
| parents | 00caf63edcdf |
| children | 6e81f4b37e13 |
| files | 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/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_event_handler.m scripts/ode/private/ode_rk_interpolate.m scripts/ode/private/ode_struct_value_check.m scripts/ode/private/odepkg_event_handle.m scripts/ode/private/runge_kutta_45_dorpri.m scripts/ode/private/runge_kutta_interpolate.m scripts/ode/private/starting_stepsize.m |
| diffstat | 16 files changed, 868 insertions(+), 903 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ode/module.mk Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/module.mk Sun Oct 18 09:55:41 2015 -0700 @@ -8,10 +8,10 @@ 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/ode_rk_interpolate.m \ + scripts/ode/private/ode_event_handler.m \ scripts/ode/private/ode_struct_value_check.m \ scripts/ode/private/runge_kutta_45_dorpri.m \ + scripts/ode/private/runge_kutta_interpolate.m \ scripts/ode/private/starting_stepsize.m scripts_ode_FCN_FILES = \
--- a/scripts/ode/ode45.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/ode45.m Sun Oct 18 09:55:41 2015 -0700 @@ -20,388 +20,351 @@ ## -*- texinfo -*- ## @deftypefn {Function File} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init}) -## @deftypefnx {Function File} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init}, @var{opt}) +## @deftypefnx {Function File} {[@var{t}, @var{y}] =} ode45 (@var{fun}, @var{trange}, @var{init}, @var{ode_opt}) ## @deftypefnx {Function File} {[@var{t}, @var{y}] =} ode45 (@dots{}, @var{par1}, @var{par2}, @dots{}) -## @deftypefnx {Function File} {[@var{t}, @var{y}, @var{xe}, @var{ye}, @var{ie}] =} ode45 (@dots{}) -## @deftypefnx {Function File} {@var{solution} =} ode45 (@var{fun}, @var{trange}, @var{init}, @dots{}) +## @deftypefnx {Function File} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode45 (@dots{}) +## @deftypefnx {Function File} {@var{solution} =} ode45 (@dots{}) ## ## Solve a set of non-stiff Ordinary Differential Equations (non-stiff ODEs) ## with the well known explicit Dormand-Prince method of order 4. ## -## The first input argument must be a function handle or inline function that -## defines the ODE: @code{y' = f(t,y)}. The function must accept two inputs -## where the first is time @var{t} and the second is a column vector of -## unknowns @var{y}. +## @var{fun} is a function handle, inline function, or string containing the +## name of the function that defines the ODE: @code{y' = f(t,y)}. The function +## must accept two inputs where the first is time @var{t} and the second is a +## column vector of unknowns @var{y}. ## ## @var{trange} specifies the time interval over which the ODE will be -## evaluated. Usually, it is a two-element vector specifying the initial and +## evaluated. Typically, it is a two-element vector specifying the initial and ## final times (@code{[tinit, tfinal]}). If there are more than two elements ## then the solution will also be evaluated at these intermediate time -## instances unless the integrate function called is -## @command{integrate_n_steps}. If there is only one time value, then -## @code{ode45} will raise an error unless the options structure has -## non-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 @var{"TimeStepNumber"} is also -## specified then the integrate function @command{integrate_n_steps} will be -## used; otherwise, @command{integrate_adaptive} is used. For this last -## possibility the user can set the tolerance for the timestep computation by -## changing the option @var{"Tau"}, that has a default value of @math{1e-6}. +## instances unless the integrate function specified is +## @command{integrate_n_steps}. ## -## The third input argument @var{init} contains the initial value for the -## unknowns. If this is a row 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}. +## By default, @code{ode45} uses an adaptive timestep with the +## @code{integrate_adaptive} algorithm. The tolerance for the timestep +## computation may be changed by using the option @qcode{"Tau"}, that has a +## default value of @math{1e-6}. If the ODE option @qcode{"TimeStepSize"} is +## not empty, then the stepper called will be @code{integrate_const}. If, in +## addition, the option @qcode{"TimeStepNumber"} is also specified then the +## integrate function @code{integrate_n_steps} will be used. ## -## If present, the fourth input argument specifies options to the ODE solver. -## It is a structure typically generated by @code{odeset}. +## @var{init} contains the initial value for the unknowns. If it is a row +## 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 function usually produces just two outputs. Variable @var{t} is a +## The optional fourth argument @var{ode_opt} specifies non-default options to +## the ODE solver. It is a structure generated by @code{odeset}. +## +## The function typically returns two outputs. Variable @var{t} is a ## column vector and contains the times where the solution was found. The ## output @var{y} is a matrix in which each column refers to a different ## unknown of the problem and each row corresponds to a time in @var{t}. ## -## For example, solve an anonymous implementation of the Van der Pol equation +## The output can also be returned as a structure @var{solution} which +## has field @var{x} containing the time where the solution was evaluated and +## field @var{y} containing the solution matrix for the times in @var{x}. +## Use @code{fieldnames (@var{solution})} to see the other fields and additional +## information returned. +## +## If using the @qcode{"Events"} option then three additional outputs may +## be returned. @var{te} holds the time when an Event function returned a +## zero. @var{ye} holds the value of the solution at time @var{te}. @var{ie} +## contains an index indicating which Event function was triggered in the case +## of multiple Event functions. +## +## Example: Solve the Van der Pol equation ## ## @example ## @group -## fvdp = @@(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; -## [T,Y] = ode45 (fvdp, [0 20], [2 0]); +## fvdp = @@(@var{t},@var{y}) [@var{y}(2); (1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)]; +## [@var{t},@var{y}] = ode45 (fvdp, [0 20], [2 0]); ## @end group ## @end example ## @seealso{odeset, odeget} ## @end deftypefn -function varargout = ode45 (vfun, vtrange, vinit, varargin) +function varargout = ode45 (fun, trange, init, varargin) if (nargin < 3) print_usage (); endif - - vorder = 5; # runge_kutta_45_dorpri uses local extrapolation - vsolver = "ode45"; + + order = 5; # runge_kutta_45_dorpri uses local extrapolation + solver = "ode45"; if (nargin >= 4) if (! isstruct (varargin{1})) - ## varargin{1:len} are parameters for vfun - vodeoptions = odeset (); - vodeoptions.vfunarguments = varargin; + ## varargin{1:len} are parameters for fun + odeopts = odeset (); + odeopts.funarguments = varargin; elseif (length (varargin) > 1) ## varargin{1} is an ODE options structure vopt - vodeoptions = ode_struct_value_check ("ode45", varargin{1}, "ode45"); - vodeoptions.vfunarguments = {varargin{2:length(varargin)}}; + odeopts = ode_struct_value_check ("ode45", varargin{1}, "ode45"); + odeopts.funarguments = {varargin{2:length(varargin)}}; else # if (isstruct (varargin{1})) - vodeoptions = ode_struct_value_check ("ode45", varargin{1}, "ode45"); - vodeoptions.vfunarguments = {}; + odeopts = ode_struct_value_check ("ode45", varargin{1}, "ode45"); + odeopts.funarguments = {}; endif else # nargin == 3 - vodeoptions = odeset (); - vodeoptions.vfunarguments = {}; + odeopts = odeset (); + odeopts.funarguments = {}; endif - if (! isvector (vtrange) || ! isnumeric (vtrange)) + if (! isvector (trange) || ! isnumeric (trange)) error ("Octave:invalid-input-arg", - "second input argument must be a valid vector"); + "ode45: TRANGE must be a numeric vector"); endif - TimeStepNumber = odeget (vodeoptions, "TimeStepNumber", [], "fast"); - TimeStepSize = odeget (vodeoptions, "TimeStepSize", [], "fast"); - if (length (vtrange) < 2 + TimeStepNumber = odeget (odeopts, "TimeStepNumber", [], "fast"); + TimeStepSize = odeget (odeopts, "TimeStepSize", [], "fast"); + if (length (trange) < 2 && (isempty (TimeStepSize) || isempty (TimeStepNumber))) error ("Octave:invalid-input-arg", - "second input argument must be a valid vector"); - elseif (vtrange(2) == vtrange(1)) + "ode45: TRANGE must contain at least 2 elements"); + elseif (trange(1) == trange(2)) error ("Octave:invalid-input-arg", - "second input argument must be a valid vector"); + "ode45: invalid time span, TRANGE(1) == TRANGE(2)"); else - vodeoptions.vdirection = sign (vtrange(2) - vtrange(1)); + odeopts.direction = sign (trange(2) - trange(1)); endif - vtrange = vtrange(:); + trange = trange(:); - if (! isvector (vinit) || ! isnumeric (vinit)) + if (! isvector (init) || ! isnumeric (init)) error ("Octave:invalid-input-arg", - "third input argument must be a valid numerical value"); + "ode45: INIT must be a numeric vector"); endif - vinit = vinit(:); + init = init(:); - if (ischar (vfun)) + if (ischar (fun)) try - vfun = str2func (vfun); + fun = str2func (fun); catch warning (lasterr); end_try_catch endif - if (! isa (vfun, "function_handle")) + if (! isa (fun, "function_handle")) error ("Octave:invalid-input-arg", - "first input argument must be a valid function handle"); + "ode45: FUN must be a valid function handle"); endif - ## Start preprocessing, have a look which options are set in vodeoptions, + ## Start preprocessing, have a look which options are set in odeopts, ## check if an invalid or unused option is set if (isempty (TimeStepNumber) && isempty (TimeStepSize)) integrate_func = "adaptive"; - vodeoptions.vstepsizefixed = false; + odeopts.stepsizefixed = false; elseif (! isempty (TimeStepNumber) && ! isempty (TimeStepSize)) integrate_func = "n_steps"; - vodeoptions.vstepsizefixed = true; - if (sign (TimeStepSize) != vodeoptions.vdirection) + odeopts.stepsizefixed = true; + if (sign (TimeStepSize) != odeopts.direction) warning ("Octave:invalid-input-arg", - "option 'TimeStepSize' has a wrong sign", - "it will be corrected automatically"); + ["ode45: option 'TimeStepSize' has the wrong sign, ", ... + "but will be corrected automatically\n"]); TimeStepSize = -TimeStepSize; endif elseif (isempty (TimeStepNumber) && ! isempty (TimeStepSize)) integrate_func = "const"; - vodeoptions.vstepsizefixed = true; - if (sign (TimeStepSize) != vodeoptions.vdirection) + odeopts.stepsizefixed = true; + if (sign (TimeStepSize) != odeopts.direction) warning ("Octave:invalid-input-arg", - "option 'TimeStepSize' has a wrong sign", - "it will be corrected automatically"); + ["ode45: option 'TimeStepSize' has the wrong sign, ", + "but will be corrected automatically\n"]); TimeStepSize = -TimeStepSize; endif else warning ("Octave:invalid-input-arg", - "assuming an adaptive integrate function"); + "ode45: assuming an adaptive integrate function\n"); integrate_func = "adaptive"; endif - vodeoptions.TimeStepSize = TimeStepSize; - vodeoptions.TimeStepNumber = TimeStepNumber; - - ## 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 ("Octave:invalid-input-arg", - "Option 'RelTol' not set, new value %f is used", - vodeoptions.RelTol); - elseif (! isempty (vodeoptions.RelTol) && vodeoptions.vstepsizefixed) - warning ("Octave:invalid-input-arg", - "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; + if (isempty (odeopts.RelTol) && ! odeopts.stepsizefixed) + odeopts.RelTol = 1e-3; warning ("Octave:invalid-input-arg", - "Option 'AbsTol' not set, new value %f is used", - vodeoptions.AbsTol); - elseif (! isempty (vodeoptions.AbsTol) && vodeoptions.vstepsizefixed) + "ode45: option 'RelTol' not set, new value %f will be used\n", + odeopts.RelTol); + elseif (! isempty (odeopts.RelTol) && odeopts.stepsizefixed) warning ("Octave:invalid-input-arg", - "Option 'AbsTol' will be ignored if fixed time stamps are given"); - else - vodeoptions.AbsTol = vodeoptions.AbsTol(:); # Create column vector + ["ode45: option 'RelTol' is ignored", ... + " when fixed time stamps are given\n"]); 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 ("Octave:invalid-input-arg", - "Option 'NonNegative' will be ignored if mass matrix is set"); - endif - else - vodeoptions.vhavenonnegative = false; + if (isempty (odeopts.AbsTol) && ! odeopts.stepsizefixed) + odeopts.AbsTol = 1e-6; + warning ("Octave:invalid-input-arg", + "ode45: option 'AbsTol' not set, new value %f will be used\n", + odeopts.AbsTol); + elseif (! isempty (odeopts.AbsTol) && odeopts.stepsizefixed) + warning ("Octave:invalid-input-arg", + ["ode45: option 'AbsTol' is ignored", ... + " when fixed time stamps are given\n"]); + else + odeopts.AbsTol = odeopts.AbsTol(:); # Create column vector 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; + odeopts.normcontrol = strcmp (odeopts.NormControl, "on"); + + if (! isempty (odeopts.NonNegative)) + if (isempty (odeopts.Mass)) + odeopts.havenonnegative = true; + else + odeopts.havenonnegative = false; + warning ("Octave:invalid-input-arg", + ["ode45: option 'NonNegative' is ignored", ... + " when mass matrix is set\n"]); + endif + else + odeopts.havenonnegative = false; 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 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; + if (isempty (odeopts.OutputFcn) && nargout == 0) + odeopts.OutputFcn = @odeplot; + odeopts.haveoutputfunction = true; + else + odeopts.haveoutputfunction = ! isempty (odeopts.OutputFcn); endif - ## Implementation of the option Stats has been finished. - ## This option can be set by the user to another value than default value. + odeopts.haveoutputselection = ! isempty (odeopts.OutputSel); - ## 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, vtrange(1), - vinit, - vodeoptions.AbsTol, - vodeoptions.RelTol, - vodeoptions.vnormcontrol); - warning ("Octave:invalid-input-arg", - "option 'InitialStep' not set, estimated value %f is used", - vodeoptions.InitialStep); - elseif (isempty (vodeoptions.InitialStep)) - vodeoptions.InitialStep = TimeStepSize; + if (odeopts.Refine > 0) + odeopts.haverefine = true; + else + odeopts.haverefine = false; 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 (vtrange(end) - vtrange(1)) / 10; + if (isempty (odeopts.InitialStep) && strcmp (integrate_func, "adaptive")) + odeopts.InitialStep = ... + odeopts.direction * starting_stepsize (order, fun, trange(1), + init, + odeopts.AbsTol, + odeopts.RelTol, + odeopts.normcontrol); warning ("Octave:invalid-input-arg", - "Option 'MaxStep' not set, new value %f is used", - vodeoptions.MaxStep); + ["ode45: option 'InitialStep' not set,", ... + " estimated value %f will be used\n"], + odeopts.InitialStep); + elseif (isempty (odeopts.InitialStep)) + odeopts.InitialStep = TimeStepSize; 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; + if (isempty (odeopts.MaxStep) && ! odeopts.stepsizefixed) + odeopts.MaxStep = abs (trange(end) - trange(1)) / 10; + warning ("Octave:invalid-input-arg", + "ode45: option 'MaxStep' not set, new value %f will be used\n", + odeopts.MaxStep); endif + odeopts.haveeventfunction = ! isempty (odeopts.Events); + ## 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)) + if (! isempty (odeopts.Jacobian)) warning ("Octave:invalid-input-arg", - "option 'Jacobian' will be ignored by this solver"); + "ode45: option 'Jacobian' is ignored by this solver\n"); endif - if (! isequal (vodeoptions.JPattern, vodetemp.JPattern)) + if (! isempty (odeopts.JPattern)) warning ("Octave:invalid-input-arg", - "option 'JPattern' will be ignored by this solver"); + "ode45: option 'JPattern' is ignored by this solver\n"); endif - if (! isequal (vodeoptions.Vectorized, vodetemp.Vectorized)) + if (! isempty (odeopts.Vectorized)) warning ("Octave:invalid-input-arg", - "option 'Vectorized' will be ignored by this solver"); + "ode45: option 'Vectorized' is ignored by this solver\n"); 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 + if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass)) + havemasshandle = false; + mass = odeopts.Mass; # constant mass + elseif (isa (odeopts.Mass, "function_handle")) + havemasshandle = 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); + havemasshandle = false; # mass = diag (ones (length (init), 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; + massdependence = ! strcmp (odeopts.MStateDependence, "none"); + + ## Other options that are not used by this solver. + if (! isempty (odeopts.MvPattern)) + warning ("Octave:invalid-input-arg", + "ode45: option 'MvPattern' is ignored by this solver\n"); endif - ## Other options that are not used by this solver. - ## Print a warning message to tell the user that the option is ignored. - if (! isequal (vodeoptions.MvPattern, vodetemp.MvPattern)) + if (! isempty (odeopts.MassSingular)) warning ("Octave:invalid-input-arg", - "option 'MvPattern' will be ignored by this solver"); + "ode45: option 'MassSingular' is ignored by this solver\n"); endif - if (! isequal (vodeoptions.MassSingular, vodetemp.MassSingular)) + if (! isempty (odeopts.InitialSlope)) warning ("Octave:invalid-input-arg", - "option 'MassSingular' will be ignored by this solver"); + "ode45: option 'InitialSlope' is ignored by this solver\n"); endif - if (! isequal (vodeoptions.InitialSlope, vodetemp.InitialSlope)) + if (! isempty (odeopts.MaxOrder)) warning ("Octave:invalid-input-arg", - "option 'InitialSlope' will be ignored by this solver"); + "ode45: option 'MaxOrder' is ignored by this solver\n"); endif - if (! isequal (vodeoptions.MaxOrder, vodetemp.MaxOrder)) + if (! isempty (odeopts.BDF)) warning ("Octave:invalid-input-arg", - "option 'MaxOrder' will be ignored by this solver"); - endif - - if (! isequal (vodeoptions.BDF, vodetemp.BDF)) - warning ("Octave:invalid-input-arg", - "option 'BDF' will be ignored by this solver"); + "ode45: option 'BDF' is ignored by this solver\n"); 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{:}); + + if (havemasshandle) # Handle only the dynamic mass matrix, + if (massdependence) # constant mass matrices have already + mass = @(t,x) odeopts.Mass (t, x, odeopts.funarguments{:}); + fun = @(t,x) mass (t, x, odeopts.funarguments{:}) ... + \ fun (t, x, odeopts.funarguments{:}); + else # if (massdependence == false) + mass = @(t) odeopts.Mass (t, odeopts.funarguments{:}); + fun = @(t,x) mass (t, odeopts.funarguments{:}) ... + \ fun (t, x, odeopts.funarguments{:}); endif endif switch (integrate_func) case "adaptive" solution = integrate_adaptive (@runge_kutta_45_dorpri, - vorder, vfun, vtrange, vinit, SubOpts); + order, fun, trange, init, odeopts); case "n_steps" solution = integrate_n_steps (@runge_kutta_45_dorpri, - vfun, vtrange(1), vinit, - TimeStepSize, TimeStepNumber, SubOpts); + fun, trange(1), init, + TimeStepSize, TimeStepNumber, odeopts); case "const" solution = integrate_const (@runge_kutta_45_dorpri, - vfun, vtrange, vinit, - TimeStepSize, SubOpts); + fun, trange, init, + TimeStepSize, odeopts); 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{:}); + if (odeopts.haveoutputfunction) # Cleanup plotter + feval (odeopts.OutputFcn, solution.t(end), + solution.x(end,:)', "done", odeopts.funarguments{:}); endif - if (vodeoptions.vhaveeventfunction) # Cleanup event function handling - odepkg_event_handle (vodeoptions.Events, solution.t(end), + if (odeopts.haveeventfunction) # Cleanup event function handling + ode_event_handler (odeopts.Events, solution.t(end), solution.x(end,:)', "done", - vodeoptions.vfunarguments{:}); + odeopts.funarguments{:}); 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)-(vnsteps)+1; # vcntcycl from 1..end - vnfevals = 6 * (solution.vcntcycles-1) + 1; # number of ode evaluations - vndecomps = 0; # number of LU decompositions - vnpds = 0; # number of partial derivatives - vnlinsols = 0; # no. of linear systems solutions + if (strcmp (odeopts.Stats, "on")) + havestats = true; + nsteps = solution.cntloop-2; # cntloop from 2..end + nfailed = (solution.cntcycles-1)-(nsteps)+1; # cntcycl from 1..end + nfevals = 6 * (solution.cntcycles-1) + 1; # number of ode evaluations + ndecomps = 0; # number of LU decompositions + npds = 0; # number of partial derivatives + nlinsols = 0; # no. of linear systems solutions ## Print cost statistics if no output argument is given if (nargout == 0) - printf ("Number of successful steps: %d\n", vnsteps); - printf ("Number of failed attempts: %d\n", vnfailed); - printf ("Number of function calls: %d\n", vnfevals); + printf ("Number of successful steps: %d\n", nsteps); + printf ("Number of failed attempts: %d\n", nfailed); + printf ("Number of function calls: %d\n", nfevals); endif else - vhavestats = false; + havestats = false; endif if (nargout == 2) @@ -410,29 +373,29 @@ elseif (nargout == 1) 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 occurred - varargout{1}.xe = solution.vevent{3}; # Time info when an event occurred - varargout{1}.ye = solution.vevent{4}; # Results when an event occurred + varargout{1}.solver = solver; # Solver name is saved in field solver + if (odeopts.haveeventfunction) + varargout{1}.ie = solution.event{2}; # Index info which event occurred + varargout{1}.xe = solution.event{3}; # Time info when an event occurred + varargout{1}.ye = solution.event{4}; # Results when an event occurred endif - if (vhavestats) + if (havestats) 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; + varargout{1}.stats.nsteps = nsteps; + varargout{1}.stats.nfailed = nfailed; + varargout{1}.stats.nfevals = nfevals; + varargout{1}.stats.npds = npds; + varargout{1}.stats.ndecomps = ndecomps; + varargout{1}.stats.nlinsols = nlinsols; endif elseif (nargout == 5) varargout = cell (1,5); varargout{1} = solution.t; varargout{2} = solution.x; - if (vodeoptions.vhaveeventfunction) - varargout{3} = solution.vevent{3}; # Time info when an event occurred - varargout{4} = solution.vevent{4}; # Results when an event occurred - varargout{5} = solution.vevent{2}; # Index info which event occurred + if (odeopts.haveeventfunction) + varargout{3} = solution.event{3}; # Time info when an event occurred + varargout{4} = solution.event{4}; # Results when an event occurred + varargout{5} = solution.event{2}; # Index info which event occurred endif endif @@ -443,210 +406,211 @@ ## 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)]; +%!function ydot = fpol (t, y) # The Van der Pol +%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; %!endfunction -%!function [vref] = fref () # The computed reference solution -%! vref = [0.32331666704577, -1.83297456798624]; +%!function ref = fref () # The computed reference solution +%! ref = [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 jac = fjac (t, y, varargin) # its Jacobian +%! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]; %!endfunction -%!function [vjac] = fjcc (vt, vy, varargin) # sparse type -%! vjac = sparse ([0, 1; -1 - 2 * vy(1) * vy(2), 1 - vy(1)^2]) +%!function jac = fjcc (t, y, varargin) # sparse type +%! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]) %!endfunction -%!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 [val, trm, dir] = feve (t, y, varargin) +%! val = fpol (t, y, varargin); # We use the derivatives +%! trm = zeros (2,1); # that's why component 2 +%! dir = ones (2,1); # seems to not be exact %!endfunction -%!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 [val, trm, dir] = fevn (t, y, varargin) +%! val = fpol (t, y, varargin); # We use the derivatives +%! trm = ones (2,1); # that's why component 2 +%! dir = ones (2,1); # seems to not be exact %!endfunction -%!function [vmas] = fmas (vt, vy, varargin) -%! vmas = [1, 0; 0, 1]; # Dummy mass matrix for tests +%!function mas = fmas (t, y, varargin) +%! mas = [1, 0; 0, 1]; # Dummy mass matrix for tests %!endfunction -%!function [vmas] = fmsa (vt, vy, varargin) -%! vmas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix +%!function mas = fmsa (t, y, varargin) +%! mas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix %!endfunction -%!function [vout] = fout (vt, vy, vflag, varargin) -%! if (regexp (char (vflag), 'init') == 1) -%! if (any (size (vt) != [2, 1])) error ('"fout" step "init"'); endif -%! elseif (isempty (vflag)) -%! if (any (size (vt) != [1, 1])) error ('"fout" step "calc"'); endif -%! vout = false; -%! elseif (regexp (char (vflag), 'done') == 1) -%! if (any (size (vt) != [1, 1])) error ('"fout" step "done"'); endif +%!function out = fout (t, y, flag, varargin) +%! if (regexp (char (flag), 'init') == 1) +%! if (any (size (t) != [2, 1])) error ('"fout" step "init"'); endif +%! elseif (isempty (flag)) +%! if (any (size (t) != [1, 1])) error ('"fout" step "calc"'); endif +%! out = false; +%! elseif (regexp (char (flag), 'done') == 1) +%! if (any (size (t) != [1, 1])) error ('"fout" step "done"'); endif %! else -%! error ('"fout" invalid vflag'); +%! error ('"fout" invalid flag'); %! endif %!endfunction %! -%! ## Turn off output of warning messages for all tests, turn them on -%! ## again if the last test is called -%!error # ouput argument -%! warning ("off", "Octave:invalid-input-arg"); -%! B = ode45 (1, [0 25], [3 15 1]); -%!error # input argument number one -%! [vt, vy] = ode45 (1, [0 25], [3 15 1]); +%! ## Turn off output of warning messages for all tests, +%! ## turn them on again after the last test is called +%!test +%! warning ("off", "Octave:invalid-input-arg", "local"); +%!error # first input arg is not a function +%! [t, y] = ode45 (1, [0 25], [3 15 1]); %!error # input argument number one as name of non existing function -%! [vt, vy] = ode45 ("non-existing-function", [0 25], [3 15 1]); +%! [t, y] = ode45 ("non-existing-function", [0 25], [3 15 1]); %!error # input argument number two -%! [vt, vy] = ode45 (@fpol, 1, [3 15 1]); +%! [t, y] = 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); +%! [t, y] = ode45 (@fpol, [0 2], [2 0]); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); %!test # not too many steps -%! [vt, vy] = ode45 (@fpol, [0 2], [2 0]); -%! assert (size (vt) < 20); +%! [t, y] = ode45 (@fpol, [0 2], [2 0]); +%! assert (size (t) < 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); +%! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%! [t, y] = ode45 (fvdb, [0 2], [2 0]); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); %!test # string instead of function -%! [vt, vy] = ode45 ("fpol", [0 2], [2 0]); -%! assert ([vt(end), vy(end,:)], [2, fref], 1e-2); +%! [t, y] = ode45 ("fpol", [0 2], [2 0]); +%! assert ([t(end), y(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); +%! [t, y] = ode45 (@fpol, [0 2], [2 0], 12, 13, "KL"); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); %!test # empty ODEOPT 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); +%! opt = odeset; +%! [t, y] = ode45 (@fpol, [0 2], [2 0], opt, 12, 13, "KL"); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); %!error # strange ODEOPT structure -%! vopt = struct ("foo", 1); -%! [vt, vy] = ode45 (@fpol, [0 2], [2 0], vopt); +%! opt = struct ("foo", 1); +%! [t, y] = ode45 (@fpol, [0 2], [2 0], opt); %!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); +%! opt = odeset("TimeStepSize", 0.1); +%! [t, y] = ode45 (@fpol, [0,2], [2 0], opt); +%! assert (t(:), [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); +%! [t, y] = ode45 (@(t,y) y, [-2 0], 2); +%! assert ([t(end), y(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); +%! opt = odeset ("InitialStep", 1e-8); +%! [t, y] = ode45 (@fpol, [0 0.2], [2 0], opt); +%! assert ([t(2)-t(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); +%! opt = odeset ("MaxStep", 1e-3); +%! sol = ode45 (@fpol, [0 0.2], [2 0], opt); +%! assert ([sol.x(5)-sol.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); +%! sol = ode45 (@fpol, [0 1 2], [2 0]); +%! assert (any((sol.x-1) == 0)); +%! assert ([sol.x(end), sol.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); +%! sol = ode45 (@fpol, [0 -2], [2 0]); +%! assert ([sol.x(end), sol.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); +%! sol = ode45 (@fpol, [2 -2], fref); +%! assert ([sol.x(end), sol.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); +%! sol = ode45 (@fpol, [2 0 -2], fref); +%! idx = find(sol.x < 0, 1, "first") - 1; +%! assert ([sol.x(idx), sol.y(idx,:)], [0 2 0], 1e-2); +%! assert ([sol.x(end), sol.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); +%! sol = ode45 (@(t,y) y, [0 -1], 1); +%! assert ([sol.x(end), sol.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); +%! sol = ode45 (@(t,y) y, [-2 0], 2); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("MaxStep", 1e-3); +%! sol = ode45 (@fpol, [0 -2], [2 0], opt); +%! assert ([abs(sol.x(8)-sol.x(7))], [1e-3], 1e-3); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("AbsTol", 1e-5); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("RelTol", 1e-8, "NormControl", "on"); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("NonNegative", 2); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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 # 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")); +%! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%!test # Stats must add further elements in sol +%! opt = odeset ("Stats", "on"); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert (isfield (sol, "stats")); +%! assert (isfield (sol.stats, "nsteps")); +%!test # Events option add further elements in sol +%! opt = odeset ("Events", @feve); +%! sol = ode45 (@fpol, [0 10], [2 0], opt); +%! assert (isfield (sol, "ie")); +%! assert (sol.ie(1), 2); +%! assert (isfield (sol, "xe")); +%! assert (isfield (sol, "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], +%! opt = odeset ("Events", @fevn, "NormControl", "on"); +%! sol = ode45 (@fpol, [0 10], [2 0], opt); +%! assert ([sol.ie, sol.xe, sol.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], +%! opt = odeset ("Events", @fevn, "NormControl", "on"); +%! [t, y, vxe, ye, vie] = ode45 (@fpol, [0 10], [2 0], opt); +%! assert ([vie, vxe, ye], %! [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); +%! opt = odeset ("Jacobian", @fjac); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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 -%! +%! opt = odeset ("Jacobian", @fjcc); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); + +## test for JPattern option is missing +## test for Vectorized 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); +%! opt = odeset ("Mass", @fmas); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("Mass", eye (2,2)); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("Mass", sparse (eye (2,2))); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("Mass", @fmsa); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.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); +%! opt = odeset ("Mass", @fmas, "MStateDependence", "strong"); +%! sol = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%!test # Set BDF option to something other than default +%! opt = odeset ("BDF", "on"); +%! [t, y] = ode45 (@fpol, [0 2], [2 0], opt); +%! assert ([t(end), y(end,:)], [2, fref], 1e-3); + +## test for MvPattern option is missing +## test for InitialSlope option is missing +## test for MaxOrder option is missing + %!test -%!## test for MvPattern option is missing -%!## test for InitialSlope option is missing -%!## test for MaxOrder option is missing %! -%! warning ("on", "Octave:invalid-input-arg"); +%! #warning ("on", "Octave:invalid-input-arg");
--- a/scripts/ode/odeget.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/odeget.m Sun Oct 18 09:55:41 2015 -0700 @@ -31,7 +31,7 @@ ## ## If called called with an optional third input argument, and @var{field} is ## not set in the structure @var{ode_opt}, then return the default value -## @var{default} instead. +## @var{default} instead. ## @seealso{odeset} ## @end deftypefn @@ -90,7 +90,7 @@ "MaxStep"; "MStateDependence"; "MvPattern"; "NonNegative"; "NormControl"; "OutputFcn"; "OutputSel"; "Refine"; "RelTol"; "Stats"; "Vectorized"}; - + exactmatch = true; match = find (strcmpi (field, options)); if (isempty (match))
--- a/scripts/ode/odeset.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/odeset.m Sun Oct 18 09:55:41 2015 -0700 @@ -218,7 +218,7 @@ ## function to print all possible options function print_options () - + disp ("List of all possible ODE solver options."); disp ("Default values are in square brackets."); disp (""); @@ -255,7 +255,7 @@ %! odeoptA = odeset (); %!demo -%! # A new ODE options structure with manually set options +%! # A new ODE options structure with manually set options %! # for "AbsTol" and "RelTol" is created. %! %! odeoptB = odeset ("AbsTol", 1e-2, "RelTol", 1e-1); @@ -272,7 +272,7 @@ %!test %! odeoptA = odeset (); %! assert (isstruct (odeoptA)); -%! assert (numel (fieldnames (odeoptA)), 23); +%! assert (numel (fieldnames (odeoptA)), 23); %! assert (all (structfun ("isempty", odeoptA))); %!shared odeoptB, odeoptC @@ -288,7 +288,7 @@ %!test %! odeoptD = odeset (odeoptB, odeoptC); -%! assert (odeoptD, odeoptC); +%! assert (odeoptD, odeoptC); ## Test custom user-defined option %!test
--- a/scripts/ode/private/AbsRel_Norm.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/private/AbsRel_Norm.m Sun Oct 18 09:55:41 2015 -0700 @@ -28,12 +28,12 @@ if (length (x_old) != n || length (y) != n) error ("Octave:invalid-input-arg", "invalid dimensions of input arguments"); endif - + if ((length (AbsTol) != 1 && length (AbsTol) != n) || (length (RelTol) != 1 && length (RelTol) != n)) error ("Octave:invalid-input-arg", "invalid dimensions of input arguments"); endif - + sc = AbsTol + max (abs (x), abs (x_old)) .* RelTol; if (normcontrol) retval = max (abs (x - y) ./ sc);
--- a/scripts/ode/private/integrate_adaptive.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/private/integrate_adaptive.m Sun Oct 18 09:55:41 2015 -0700 @@ -73,10 +73,10 @@ dt = odeget (options, "InitialStep", [], "fast"); if (isempty (dt)) dt = starting_stepsize (order, func, t, x, options.AbsTol, options.RelTol, - options.vnormcontrol); + options.normcontrol); endif - dir = odeget (options, "vdirection", [], "fast"); + dir = odeget (options, "direction", [], "fast"); dt = dir * min (abs (dt), options.MaxStep); options.comp = 0.0; @@ -87,30 +87,30 @@ fac = 0.38^(1/(order+1)); # formula taken from Hairer ## Initialize the OutputFcn - if (options.vhaveoutputfunction) - if (options.vhaveoutputselection) - solution.vretout = x(options.OutputSel,end); + if (options.haveoutputfunction) + if (options.haveoutputselection) + solution.retout = x(options.OutputSel,end); else - solution.vretout = x; + solution.retout = x; endif - feval (options.OutputFcn, tspan, solution.vretout, - "init", options.vfunarguments{:}); + feval (options.OutputFcn, tspan, solution.retout, + "init", options.funarguments{:}); endif ## Initialize the EventFcn - if (options.vhaveeventfunction) - odepkg_event_handle (options.Events, tspan(end), x, - "init", options.vfunarguments{:}); + if (options.haveeventfunction) + ode_event_handler (options.Events, tspan(end), x, + "init", options.funarguments{:}); endif - if (options.vhavenonnegative) + if (options.havenonnegative) nn = options.NonNegative; endif - solution.vcntloop = 2; - solution.vcntcycles = 1; - solution.vcntsave = 2; - solution.vunhandledtermination = true; + solution.cntloop = 2; + solution.cntcycles = 1; + solution.cntsave = 2; + solution.unhandledtermination = true; ireject = 0; k_vals = []; @@ -122,20 +122,20 @@ [t_new, x_new, x_est, new_k_vals] = ... stepper (func, t_old, x_old, dt, options, k_vals, t_new); - solution.vcntcycles++; + solution.cntcycles++; - if (options.vhavenonnegative) - x_new(nn, end) = abs (x_new(nn, end)); + if (options.havenonnegative) + x_new(nn, end) = abs (x_new(nn, end)); x_est(nn, end) = abs (x_est(nn, end)); endif err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol, - options.vnormcontrol, x_est); + options.normcontrol, x_est); ## Accept solution only if err <= 1.0 if (err <= 1) - solution.vcntloop++; + solution.cntloop++; ireject = 0; ## if output time steps are fixed @@ -149,28 +149,28 @@ t(t_caught) = tspan(t_caught); iout = max (t_caught); x(:, t_caught) = ... - ode_rk_interpolate (order, [t_old t_new], [x_old x_new], + runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], tspan(t_caught), new_k_vals, dt, - options.vfunarguments{:}); + options.funarguments{:}); istep++; ## Call Events function only if a valid result has been found. - ## Stop integration if veventbreak is true. - if (options.vhaveeventfunction) + ## Stop integration if eventbreak is true. + if (options.haveeventfunction) break_loop = false; for idenseout = 1:numel (t_caught) id = t_caught(idenseout); td = t(id); - solution.vevent = ... - odepkg_event_handle (options.Events, t(id), x(:, id), [], - options.vfunarguments{:}); - if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) - t(id) = solution.vevent{3}(end); + solution.event = ... + ode_event_handler (options.Events, t(id), x(:, id), [], + options.funarguments{:}); + if (! isempty (solution.event{1}) && solution.event{1} == 1) + t(id) = solution.event{3}(end); t = t(1:id); - x(:, id) = solution.vevent{4}(end, :).'; + x(:, id) = solution.event{4}(end, :).'; x = x(:,1:id); - solution.vunhandledtermination = false; + solution.unhandledtermination = false; break_loop = true; break; endif @@ -182,22 +182,22 @@ ## Call OutputFcn only if a valid result has been found. ## Stop integration if function returns false. - if (options.vhaveoutputfunction) - vcnt = options.Refine + 1; - vapproxtime = linspace (t_old, t_new, vcnt); - vapproxvals = interp1 ([t_old, t(t_caught), t_new], - [x_old, x(:, t_caught), x_new] .', - vapproxtime, 'linear') .'; - if (options.vhaveoutputselection) - vapproxvals = vapproxvals(options.OutputSel, :); + if (options.haveoutputfunction) + cnt = options.Refine + 1; + approxtime = linspace (t_old, t_new, cnt); + approxvals = interp1 ([t_old, t(t_caught), t_new], + [x_old, x(:, t_caught), x_new] .', + approxtime, 'linear') .'; + if (options.haveoutputselection) + approxvals = approxvals(options.OutputSel, :); endif - for ii = 1:numel (vapproxtime) - vpltret = feval (options.OutputFcn, vapproxtime(ii), - vapproxvals(:, ii), [], - options.vfunarguments{:}); + for ii = 1:numel (approxtime) + pltret = feval (options.OutputFcn, approxtime(ii), + approxvals(:, ii), [], + options.funarguments{:}); endfor - if (vpltret) # Leave main loop - solution.vunhandledtermination = false; + if (pltret) # Leave main loop + solution.unhandledtermination = false; break; endif endif @@ -211,36 +211,36 @@ iout = istep; ## Call Events function only if a valid result has been found. - ## Stop integration if veventbreak is true. - if (options.vhaveeventfunction) - solution.vevent = ... - odepkg_event_handle (options.Events, t(istep), x(:, istep), [], - options.vfunarguments{:}); - if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) - t(istep) = solution.vevent{3}(end); - x(:, istep) = solution.vevent{4}(end, :).'; - solution.vunhandledtermination = false; + ## Stop integration if eventbreak is true. + if (options.haveeventfunction) + solution.event = ... + ode_event_handler (options.Events, t(istep), x(:, istep), [], + options.funarguments{:}); + if (! isempty (solution.event{1}) && solution.event{1} == 1) + t(istep) = solution.event{3}(end); + x(:, istep) = solution.event{4}(end, :).'; + solution.unhandledtermination = false; break; endif endif ## Call OutputFcn only if a valid result has been found. ## Stop integration if function returns false. - if (options.vhaveoutputfunction) - vcnt = options.Refine + 1; - vapproxtime = linspace (t_old, t_new, vcnt); - vapproxvals = interp1 ([t_old, t_new], - [x_old, x_new] .', - vapproxtime, 'linear') .'; - if (options.vhaveoutputselection) - vapproxvals = vapproxvals(options.OutputSel, :); + if (options.haveoutputfunction) + cnt = options.Refine + 1; + approxtime = linspace (t_old, t_new, cnt); + approxvals = interp1 ([t_old, t_new], + [x_old, x_new] .', + approxtime, 'linear') .'; + if (options.haveoutputselection) + approxvals = approxvals(options.OutputSel, :); endif - for ii = 1:numel (vapproxtime) - vpltret = feval (options.OutputFcn, vapproxtime(ii), - vapproxvals(:, ii), [], options.vfunarguments{:}); + for ii = 1:numel (approxtime) + pltret = feval (options.OutputFcn, approxtime(ii), + approxvals(:, ii), [], options.funarguments{:}); endfor - if (vpltret) # Leave main loop - solution.vunhandledtermination = false; + if (pltret) # Leave main loop + solution.unhandledtermination = false; break; endif endif @@ -252,8 +252,7 @@ x_old = x_new; k_vals = new_k_vals; - solution.vcntloop += 1; - vcntiter = 0; + solution.cntloop += 1; else @@ -286,22 +285,22 @@ ## Check if integration of the ode has been successful if (dir * t(end) < dir * tspan(end)) - if (solution.vunhandledtermination == true) + if (solution.unhandledtermination == true) error ("integrate_adaptive:unexpected_termination", [" Solving was not successful. ", ... " The iterative integration loop exited at time", ... " t = %f before the endpoint at tend = %f was reached. ", ... " This may happen if the stepsize becomes too small. ", ... " Try to reduce the value of 'InitialStep'", ... - " and/or 'MaxStep' with the command 'odeset'.\n"], + " and/or 'MaxStep' with the command 'odeset'."], t(end), tspan(end)); else warning ("integrate_adaptive:unexpected_termination", ["Solver was stopped by a call of 'break'", ... - " in the main iteration loop at time ", ... + " in the main iteration loop at time", ... " t = %f before the endpoint at tend = %f was reached. ", ... - " This may happen because the @odeplot function ", ... - " returned 'true' or the @event function returned 'true'.\n"], + " This may happen because the @odeplot function", ... + " returned 'true' or the @event function returned 'true'."], t(end), tspan(end)); endif endif
--- a/scripts/ode/private/integrate_const.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/private/integrate_const.m Sun Oct 18 09:55:41 2015 -0700 @@ -69,8 +69,8 @@ t = tspan(1); x = x0(:); - vdirection = odeget (options, "vdirection", [], "fast"); - if (sign (dt) != vdirection) + direction = odeget (options, "direction", [], "fast"); + if (sign (dt) != direction) error ("Octave:invalid-input-arg", "option 'InitialStep' has a wrong sign"); endif @@ -81,33 +81,33 @@ 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; + if (options.haveoutputfunction) + if (options.haveoutputselection) + solution.retout = x(options.OutputSel,end); + else + solution.retout = x; endif - feval (options.OutputFcn, tspan, solution.vretout, "init", - options.vfunarguments{:}); + feval (options.OutputFcn, tspan, solution.retout, "init", + options.funarguments{:}); endif ## Initialize the EventFcn - if (options.vhaveeventfunction) - odepkg_event_handle (options.Events, t(end), x, "init", - options.vfunarguments{:}); + if (options.haveeventfunction) + ode_event_handler (options.Events, t(end), x, "init", + options.funarguments{:}); endif - - solution.vcntloop = 2; - solution.vcntcycles = 1; - vcntiter = 0; - solution.vunhandledtermination = true; - solution.vcntsave = 2; - + + solution.cntloop = 2; + solution.cntcycles = 1; + cntiter = 0; + solution.unhandledtermination = true; + solution.cntsave = 2; + z = t; u = x; - k_vals = feval (func, t , x, options.vfunarguments{:}); + k_vals = feval (func, t , x, options.funarguments{:}); while (counter <= k) ## computing the integration step from t to t+dt @@ -116,15 +116,15 @@ [tk, comp] = kahan (tk,comp, dt); options.comp = comp; s(end) = tk; - - if (options.vhavenonnegative) + + if (options.havenonnegative) 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); + + if (options.haveoutputfunction && options.haverefine) + SaveVUForRefine = u(:,end); endif ## values on this interval for time and solution @@ -139,7 +139,7 @@ ## 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) ) + elseif (direction * z(end) > direction * tspan(counter) ) ## initializing counter for the following cycle i = 2; while (i <= length (z)) @@ -153,7 +153,7 @@ endif ## else, loop until there are requested values inside this subinterval while ((counter <= k) - && vdirection * z(i) > vdirection * tspan(counter) ) + && direction * z(i) > direction * 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)]; @@ -166,7 +166,7 @@ endwhile ## if new time requested is not out of this interval - if (counter <= k && vdirection * z(end) > vdirection * tspan(counter)) + if (counter <= k && direction * z(end) > direction * tspan(counter)) ## update the counter i++; else @@ -179,65 +179,65 @@ x = [x,u(:,2:end)]; t = [t;z(2:end)]; - solution.vcntsave += 1; - solution.vcntloop += 1; - vcntiter = 0; - + solution.cntsave += 1; + solution.cntloop += 1; + cntiter = 0; + ## Call OutputFcn only if a valid result has been found. ## Stop integration if 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; + if (options.haveoutputfunction) + for cnt = 0:options.Refine # Approximation between told and t + if (options.haverefine) # Do interpolation + approxtime = (cnt + 1) / (options.Refine + 2); + approxvals = (1 - approxtime) * SaveVUForRefine ... + + (approxtime) * y(:,end); + approxtime = s(end) + approxtime*dt; else - vapproxvals = x(:,end); - vapproxtime = t(end); + approxvals = x(:,end); + approxtime = t(end); endif - if (options.vhaveoutputselection) - vapproxvals = vapproxvals(options.OutputSel); + if (options.haveoutputselection) + approxvals = approxvals(options.OutputSel); endif - vpltret = feval (options.OutputFcn, vapproxtime, vapproxvals, [], - options.vfunarguments{:}); - if (vpltret) # Leave refinement loop + pltret = feval (options.OutputFcn, approxtime, approxvals, [], + options.funarguments{:}); + if (pltret) # Leave refinement loop break; endif endfor - if (vpltret) # Leave main loop - solution.vunhandledtermination = false; + if (pltret) # Leave main loop + solution.unhandledtermination = false; break; endif endif - + ## Call Events function only if a valid result has been found. - ## 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; + ## Stop integration if eventbreak is true. + if (options.haveeventfunction) + solution.event = ode_event_handler (options.Events, t(end), x(:,end), + [], options.funarguments{:}); + if (! isempty (solution.event{1}) && solution.event{1} == 1) + t(solution.cntloop-1,:) = solution.event{3}(end,:); + x(:,solution.cntloop-1) = solution.event{4}(end,:)'; + solution.unhandledtermination = false; break; endif endif - + ## Update counters that count the number of iteration cycles - solution.vcntcycles += 1; # Needed for cost statistics - vcntiter += 1; # Needed to find iteration problems + solution.cntcycles += 1; # Needed for cost statistics + cntiter += 1; # Needed to find iteration problems ## Stop solving because, in the last 5,000 steps, no successful valid ## value has been found - if (vcntiter >= 5_000) + if (cntiter >= 5_000) error (["integrate_const: Solving was not successful. ", ... " The iterative integration loop exited at time", ... " t = %f before the endpoint at tend = %f was reached. ", ... " This happened because the iterative integration loop", ... " did not find a valid solution at this time stamp. ", ... " Try to reduce the value of 'InitialStep' and/or", ... - " 'MaxStep' with the command 'odeset'.\n"], + " 'MaxStep' with the command 'odeset'."], s(end), tspan(end)); endif @@ -247,36 +247,36 @@ endif endwhile - + ## Check if integration of the ode has been successful - if (vdirection * z(end) < vdirection * tspan(end)) - if (solution.vunhandledtermination == true) + if (direction * z(end) < direction * tspan(end)) + if (solution.unhandledtermination == true) error ("integrate_const:unexpected_termination", [" Solving was not successful. ", ... " The iterative integration loop exited at time", ... " t = %f before the endpoint at tend = %f was reached. ", ... " This may happen if the stepsize becomes too small. ", ... " Try to reduce the value of 'InitialStep'", ... - " and/or 'MaxStep' with the command 'odeset'.\n"], + " and/or 'MaxStep' with the command 'odeset'."], z(end), tspan(end)); else warning ("integrate_const:unexpected_termination", ["Solver was stopped by a call of 'break'", ... - " in the main iteration loop at time ", ... + " in the main iteration loop at time", ... " t = %f before the endpoint at tend = %f was reached. ", ... - " This may happen because the @odeplot function ", ... - " returned 'true' or the @event function returned 'true'.\n"], + " This may happen because the @odeplot function", ... + " returned 'true' or the @event function returned 'true'."], 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); + d = direction * t((end-(j-1)):end) > direction * 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
--- a/scripts/ode/private/integrate_n_steps.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/private/integrate_n_steps.m Sun Oct 18 09:55:41 2015 -0700 @@ -67,60 +67,60 @@ solution = struct (); ## first values for time and solution - x = x0(:); + x = x0(:); t = t0; - vdirection = odeget (options, "vdirection", [], "fast"); - if (sign (dt) != vdirection) + direction = odeget (options, "direction", [], "fast"); + if (sign (dt) != direction) error ("Octave:invalid-input-arg", "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; + if (options.haveoutputfunction) + if (options.haveoutputselection) + solution.retout = x(options.OutputSel,end); + else + solution.retout = x; endif - feval (options.OutputFcn, tspan, solution.vretout, "init", - options.vfunarguments{:}); + feval (options.OutputFcn, tspan, solution.retout, "init", + options.funarguments{:}); endif ## Initialize the EventFcn - if (options.vhaveeventfunction) - odepkg_event_handle (options.Events, t(end), x, "init", - options.vfunarguments{:}); + if (options.haveeventfunction) + ode_event_handler (options.Events, t(end), x, "init", + options.funarguments{:}); endif - - solution.vcntloop = 2; - solution.vcntcycles = 1; - vcntiter = 0; - solution.vunhandledtermination = true; - solution.vcntsave = 2; - + + solution.cntloop = 2; + solution.cntcycles = 1; + cntiter = 0; + solution.unhandledtermination = true; + solution.cntsave = 2; + z = t; u = x; - k_vals = feval (func, t , x, options.vfunarguments{:}); + k_vals = feval (func, t , x, options.funarguments{:}); 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) + + if (options.havenonnegative) 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); + + if (options.haveoutputfunction && options.haverefine) + SaveVUForRefine = u(:,end); endif ## values on this interval for time and solution @@ -129,92 +129,93 @@ x = [x,u(:,2:end)]; t = [t;z(2:end)]; - solution.vcntsave += 1; - solution.vcntloop += 1; - vcntiter = 0; - + solution.cntsave += 1; + solution.cntloop += 1; + cntiter = 0; + ## Call OutputFcn only if a valid result has been found. ## Stop integration if 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; + if (options.haveoutputfunction) + for cnt = 0:options.Refine # Approximation between told and t + if (options.haverefine) # Do interpolation + approxtime = (cnt + 1) / (options.Refine + 2); + approxvals = (1 - approxtime) * SaveVUForRefine ... + + (approxtime) * y(:,end); + approxtime = s(end) + approxtime*dt; else - vapproxvals = x(:,end); - vapproxtime = t(end); + approxvals = x(:,end); + approxtime = t(end); endif - if (options.vhaveoutputselection) - vapproxvals = vapproxvals(options.OutputSel); + if (options.haveoutputselection) + approxvals = approxvals(options.OutputSel); endif - vpltret = feval (options.OutputFcn, vapproxtime, vapproxvals, [], - options.vfunarguments{:}); - if (vpltret) # Leave refinement loop + pltret = feval (options.OutputFcn, approxtime, approxvals, [], + options.funarguments{:}); + if (pltret) # Leave refinement loop break; endif endfor - if (vpltret) # Leave main loop - solution.vunhandledtermination = false; + if (pltret) # Leave main loop + solution.unhandledtermination = false; break; endif endif - + ## Call Events function only if a valid result has been found. - ## 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; + ## Stop integration if eventbreak is true. + if (options.haveeventfunction) + solution.event = ode_event_handler (options.Events, t(end), x(:,end), + [], options.funarguments{:}); + if (! isempty (solution.event{1}) && solution.event{1} == 1) + t(solution.cntloop-1,:) = solution.event{3}(end,:); + x(:,solution.cntloop-1) = solution.event{4}(end,:)'; + solution.unhandledtermination = false; break; endif endif - + ## Update counters that count the number of iteration cycles - solution.vcntcycles += 1; # Needed for cost statistics - vcntiter += 1; # Needed to find iteration problems + solution.cntcycles += 1; # Needed for cost statistics + cntiter += 1; # Needed to find iteration problems ## Stop solving because, in the last 5,000 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"], + if (cntiter >= 5_000) + error (["integrate_n_steps: Solving was not successful. ", ... + " The iterative integration loop exited at time", ... + " t = %f before the endpoint at tend = %f was reached. ", ... + " This happened because the iterative integration loop", ... + " did not find a valid solution at this time stamp. ", ... + " Try to reduce the value of 'InitialStep' and/or", ... + " 'MaxStep' with the command 'odeset'.], 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) + #if (direction * z(end) < direction * tspan(end)) + # if (solution.unhandledtermination == true) # error ("integrate_n_steps:unexpected_termination", # [" Solving was not successful. ", ... # " The iterative integration loop exited at time", ... # " t = %f before the endpoint at tend = %f was reached. ", ... # " This may happen if the stepsize becomes too small. ", ... # " Try to reduce the value of 'InitialStep'", ... - # " and/or 'MaxStep' with the command 'odeset'.\n"], + # " and/or 'MaxStep' with the command 'odeset'."], # z(end), tspan(end)); # else # warning ("integrate_n_steps:unexpected_termination", # ["Solver was stopped by a call of 'break'", ... - # " in the main iteration loop at time ", ... + # " in the main iteration loop at time", ... # " t = %f before the endpoint at tend = %f was reached. ", ... - # " This may happen because the @odeplot function ", ... - # " returned 'true' or the @event function returned 'true'.\n"], + # " This may happen because the @odeplot function", ... + # " returned 'true' or the @event function returned 'true'."], # z(end), tspan(end)); # endif #endif solution.t = t; solution.x = x'; - + endfunction
--- a/scripts/ode/private/kahan.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/private/kahan.m Sun Oct 18 09:55:41 2015 -0700 @@ -27,8 +27,8 @@ ## obtained by adding a sequence of finite precision floating point numbers, ## compared to the straightforward approach. For more details ## see @url{http://en.wikipedia.org/wiki/Kahan_summation_algorithm}. -## This function is called by @command{integrate_adaptive} and -## @command{integrate_const} to better catch equality comparisons. +## This function is called by @code{integrate_adaptive} and +## @code{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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/ode_event_handler.m Sun Oct 18 09:55:41 2015 -0700 @@ -0,0 +1,159 @@ +## 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{sol} =} ode_event_handler (@var{@@evtfun}, @var{t}, @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{@@evtfun} in the form of a function handle. +## +## The second input argument @var{t} 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 @qcode{"init"} +## then initialize internal persistent variables of the function +## @code{ode_event_handler} and return an empty cell array of size 4, +## +## @item @qcode{"calc"} +## then do the evaluation of the event function and return the solution +## @var{sol} as type cell array of size 4, +## +## @item @qcode{"done"} +## then cleanup internal variables of the function +## @code{ode_event_handler} 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 +## @code{ode_event_handler} to the event function. +## +## This function is an ODE 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 retval = ode_event_handler (evtfun, t, y, flag = "", varargin) + + ## No error handling has been implemented in this function to achieve + ## the highest performance available. + + ## retval{1} is true (to terminate) or false (to continue) + ## retval{2} is the index information for which an event occurred + ## retval{3} is the time information column vector + ## retval{4} is the line by line result information matrix + + ## These persistent variables store the results and time value from the + ## processing in the previous time stamp. + ## evtold the results from the event function + ## told the time stamp + ## yold the ODE result + ## retcell the return values cell array + ## evtcnt the counter for how often this function has been called + ## has been called + persistent evtold told yold retcell evtcnt; + + ## Call the event function if an event function has been defined to + ## initialize the internal variables of the event function and to get + ## a value for evtold. + if (strcmp (flag, "init")) + + if (! iscell (y)) + inpargs = {evtfun, t, y}; + else + inpargs = {evtfun, t, y{1}, y{2}}; + y = y{1}; # Delete cell element 2 + endif + if (nargin > 4) + inpargs = {inpargs{:}, varargin{:}}; + endif + [evtold, term, dir] = feval (inpargs{:}); + + ## FIXME: This actually seems to assume that everything must be row vectors + ## We assume that all return values must be column vectors + evtold = evtold(:)'; term = term(:)'; dir = dir(:)'; + told = t; yold = y; evtcnt = 1; retcell = cell (1,4); + + ## Process the event, i.e., + ## find the zero crossings for either a rising or falling edge + elseif (isempty (flag)) + + if (! iscell (y)) + inpargs = {evtfun, t, y}; + else + inpargs = {evtfun, t, y{1}, y{2}}; + y = y{1}; # Delete cell element 2 + endif + if (nargin > 4) + inpargs = {inpargs{:}, varargin{:}}; + endif + [evt, term, dir] = feval (inpargs{:}); + + ## FIXME: This actually seems to assume that everything must be row vectors + ## We assume that all return values must be column vectors + evt = evt(:)'; term = term(:)'; dir = dir(:)'; + + ## Check if one or more signs of the event has changed + signum = (sign (evtold) != sign (evt)); + if (any (signum)) # One or more values have changed + idx = find (signum); # Get the index of the changed values + + if (any (dir(idx) == 0)) + ## Rising or falling (both are possible) + ## Don't change anything, keep the index + elseif (any (dir(idx) == sign (evt(idx)))) + ## Detected rising or falling, need a new index + idx = find (dir == sign (evt)); + else + ## Found a zero crossing but must not be notified + idx = []; + endif + + ## Create new output values if a valid index has been found + if (! isempty (idx)) + ## Change the persistent result cell array + retcell{1} = any (term(idx)); # Stop integration or not + retcell{2}(evtcnt,1) = idx(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 + tnew = t - evt(1,idx) * (t - told) / (evt(1,idx) - evtold(1,idx)); + ynew = (y - (t - tnew) * (y - yold) / (t - told))'; + retcell{3}(evtcnt,1) = tnew; + retcell{4}(evtcnt,:) = ynew; + evtcnt += 1; + endif + + endif # Check for one or more signs ... + evtold = evt; told = t; retval = retcell; yold = y; + + elseif (strcmp (flag, "done")) # Clear this event handling function + clear ("evtold", "told", "retcell", "yold", "evtcnt"); + retcell = cell (1,4); + + endif + +endfunction +
--- a/scripts/ode/private/ode_rk_interpolate.m Fri Oct 16 10:31:02 2015 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,105 +0,0 @@ -## Copyright (C) 2015 Carlo de Falco -## Copyright (C) 2015 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/>. - -function u_interp = ode_rk_interpolate (order, z, u, t, k_vals, dt, args) - - switch (order) - - ## FIXME: Can interpolations for orders 1-4 simply be deleted? 2015-10-14. - #{ - case 1 - u_interp = linear_interpolation (z, u, t); - case 2 - if (! isempty (k_vals)) - der = k_vals(:,1); - else - der = feval (func, z(1) , u(:,1), args); - endif - u_interp = quadratic_interpolation (z, u, der, t); - case 3 - u_interp = ... - hermite_cubic_interpolation (z, u, k_vals, t); - 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: - u_interp = hermite_quartic_interpolation (z, u, k_vals, t); - - ## 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 interpolation instead"]); - der(:,1) = feval (func, z(1), u(:,1), args); - der(:,2) = feval (func, z(2), u(:,2), args); - u_interp = hermite_cubic_interpolation (z, u, der, t); - - endswitch - -endfunction - - -## The function below can be used in an ODE solver to interpolate the solution -## at the time t_out using 4th order hermite interpolation. -function x_out = hermite_quartic_interpolation (t, x, der, t_out) - - persistent coefs_u_half = ... - [(6025192743/30085553152), 0, (51252292925/65400821598), ... - (-2691868925/45128329728), (187940372067/1594534317056), ... - (-1776094331/19743644256), (11237099/235043384)].'; - - ## 4th order approximation of y in t+dt/2 as proposed by Shampine in - ## Lawrence, Shampine, "Some Practical Runge-Kutta Formulas", 1986. - dt = t(2) - t(1); - u_half = x(:,1) + (1/2) * dt * (der(:,1:7) * coefs_u_half); - - ## Rescale time on [0,1] - s = (t_out - t(1)) / dt; - - ## 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 (rows (x), length (t_out)); - x_out = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4) .* x(:,1) + ... - ( s - 4*s.^2 + 5*s.^3 - 2*s.^4) .* (dt * der(:,1)) + ... - ( 16*s.^2 - 32*s.^3 + 16*s.^4) .* u_half + ... - ( - 5*s.^2 + 14*s.^3 - 8*s.^4) .* x(:,2) + ... - ( s.^2 - 3*s.^3 + 2*s.^4) .* (dt * der(:,end)); - -endfunction -
--- a/scripts/ode/private/ode_struct_value_check.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/private/ode_struct_value_check.m Sun Oct 18 09:55:41 2015 -0700 @@ -267,7 +267,7 @@ %! ode_struct_value_check (odeset); %!demo -%! # Create the ODE options structure A with odeset and check it +%! # Create the ODE options structure A with odeset and check it %! # with ode_struct_value_check. This actually is unnecessary %! # because odeset automatically calls ode_struct_value_check before %! # returning.
--- a/scripts/ode/private/odepkg_event_handle.m Fri Oct 16 10:31:02 2015 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,158 +0,0 @@ -## 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{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 the 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 @qcode{"init"} -## then initialize internal persistent variables of the function -## @command{odepkg_event_handle} and return an empty cell array of size 4, -## -## @item @qcode{"calc"} -## then do the evaluation of the event function and return the solution -## @var{sol} as type cell array of size 4, -## -## @item @qcode{"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 ODE 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 - endif - if (nargin > 4) - vinpargs = {vinpargs{:}, varargin{:}}; - endif - [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 - endif - if (nargin > 4) - vinpargs = {vinpargs{:}, varargin{:}}; - endif - [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 = []; - endif - - ## 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; - endif - - endif # 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); - - endif - -endfunction -
--- a/scripts/ode/private/runge_kutta_45_dorpri.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/private/runge_kutta_45_dorpri.m Sun Oct 18 09:55:41 2015 -0700 @@ -83,17 +83,17 @@ k = zeros (rows (x), 7); if (! isempty (options)) # extra arguments for function evaluator - args = options.vfunarguments; + args = options.funarguments; else args = {}; endif if (! isempty (k_vals)) # k values from previous step are passed k(:,1) = k_vals(:,end); # FSAL property - else + else k(:,1) = feval (f, t, x, args{:}); endif - + k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).' , args{:}); k(:,3) = feval (f, s(3), x + k(:,1:2) * aa(3, 1:2).', args{:}); k(:,4) = feval (f, s(4), x + k(:,1:3) * aa(4, 1:3).', args{:}); @@ -111,6 +111,6 @@ cc_prime = dt * c_prime; x_est = x + k * cc_prime(:); endif - + endfunction
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/runge_kutta_interpolate.m Sun Oct 18 09:55:41 2015 -0700 @@ -0,0 +1,105 @@ +## Copyright (C) 2015 Carlo de Falco +## Copyright (C) 2015 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/>. + +function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, args) + + switch (order) + + ## FIXME: Can interpolations for orders 1-4 simply be deleted? 2015-10-14. + #{ + case 1 + u_interp = linear_interpolation (z, u, t); + case 2 + if (! isempty (k_vals)) + der = k_vals(:,1); + else + der = feval (func, z(1) , u(:,1), args); + endif + u_interp = quadratic_interpolation (z, u, der, t); + case 3 + u_interp = ... + hermite_cubic_interpolation (z, u, k_vals, t); + 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: + u_interp = hermite_quartic_interpolation (z, u, k_vals, t); + + ## 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.funarguments{:}); + ## 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 interpolation instead"]); + der(:,1) = feval (func, z(1), u(:,1), args); + der(:,2) = feval (func, z(2), u(:,2), args); + u_interp = hermite_cubic_interpolation (z, u, der, t); + + endswitch + +endfunction + + +## The function below can be used in an ODE solver to interpolate the solution +## at the time t_out using 4th order hermite interpolation. +function x_out = hermite_quartic_interpolation (t, x, der, t_out) + + persistent coefs_u_half = ... + [(6025192743/30085553152), 0, (51252292925/65400821598), ... + (-2691868925/45128329728), (187940372067/1594534317056), ... + (-1776094331/19743644256), (11237099/235043384)].'; + + ## 4th order approximation of y in t+dt/2 as proposed by Shampine in + ## Lawrence, Shampine, "Some Practical Runge-Kutta Formulas", 1986. + dt = t(2) - t(1); + u_half = x(:,1) + (1/2) * dt * (der(:,1:7) * coefs_u_half); + + ## Rescale time on [0,1] + s = (t_out - t(1)) / dt; + + ## 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 (rows (x), length (t_out)); + x_out = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4) .* x(:,1) + ... + ( s - 4*s.^2 + 5*s.^3 - 2*s.^4) .* (dt * der(:,1)) + ... + ( 16*s.^2 - 32*s.^3 + 16*s.^4) .* u_half + ... + ( - 5*s.^2 + 14*s.^3 - 8*s.^4) .* x(:,2) + ... + ( s.^2 - 3*s.^3 + 2*s.^4) .* (dt * der(:,end)); + +endfunction +
--- a/scripts/ode/private/starting_stepsize.m Fri Oct 16 10:31:02 2015 -0700 +++ b/scripts/ode/private/starting_stepsize.m Sun Oct 18 09:55:41 2015 -0700 @@ -26,7 +26,7 @@ ## equations, @var{t0} is the initial time, and @var{x0} is the initial ## condition. @var{AbsTol} and @var{RelTol} are the absolute and relative ## tolerance on the ODE integration taken from an ode options structure. -## +## ## References: ## [1] E. Hairer, S.P. Norsett and G. Wanner, ## @cite{Solving Ordinary Differential Equations I: Nonstiff Problems},