Mercurial > octave
changeset 22626:869c02fde46c stable
Further clean-up of ode functions.
* scripts/ode/private/AbsRel_Norm.m: Renamed to AbsRel_norm.m
* scripts/ode/module.mk: Add AbsRel_norm.m to build system.
* ode23.m, ode45.m: Remove extra comma from Copyright statement.
Add ode45 to @seealso links in docstring.
Add FIXME notes for questionable code.
Use numel instead of length for clarity.
Use space after function name and opening parenthesis.
Wrap long lines to less than 80 characters.
Change odemergeopts function call to match new prototype.
Use single quotes to simplify strings that contain double quotes.
Use 2-space indent in %!function blocks.
* odeget.m: Remove extra comma from Copyright statement.
Use double quote in preference to single quote.
Delete whitespace at end of lines.
* odeplot.m: Re-write docstring. Include @seealso links in docstring.
Declare all persistent variables in a single declaration.
Use in-place += operator for efficiency.
Add FIXME notes for questionable code.
* odeset.m: Remove extra comma from Copyright statement.
Add additional calling form with 1 output and 0 inputs to docstring.
Add FIXME notes for questionable code.
Delete whitespace at end of lines.
* odeset.m(print_options): Use single quotes to simplify strings with double
quotes. Put default value of option first in list.
* integrate_adaptive.m: Wrap long lines < 80 characters.
Delete whitespace at end of lines.
Correct indentation of declared values after '='.
* kahan.m: Reise docstring.
* ode_event_handler.m: Use retval in docstring to match functin prototype.
Delete unnecessary comments.
* odedefaults.m: Remove extra comma from Copyright statement.
Match variable names in docstring to function prototype.
Use space after function name and before opening parenthesis.
Delete whitespace at end of lines.
* odemergeopts.m: Remove extra comma from Copyright statement.
Delete extra space in function prototype.
Change function prototype to have "caller" as first argument to match rest of
Octave.
* runge_kutta_23.m: Clean up declaration of persistent variables.
Correct misspellings in comments.
* runge_kutta_45_dorpri.m: Put description of input arguments before output
arguments in docstring.
Clean up declaration of persistent variables.
* runge_kutta_interpolate.m: Use double quotes in preference to single quotes.
Eliminate line continuations for code that would fit on a single line.
Remove obsolete code that calls non-existent functions.
Capitalize Hermite in comments.
Cleanup declaration of persistent variables.
* starting_stepsize.m: Replace calls to AbsRel_Norm with AbsRel_norm.
| author | Rik <rik@octave.org> |
|---|---|
| date | Sat, 15 Oct 2016 22:39:29 -0700 |
| parents | 081a201b77c7 |
| children | 7b190a2f11cb |
| files | scripts/ode/module.mk scripts/ode/ode23.m scripts/ode/ode45.m scripts/ode/odeget.m scripts/ode/odeplot.m scripts/ode/odeset.m scripts/ode/private/AbsRel_Norm.m scripts/ode/private/AbsRel_norm.m scripts/ode/private/integrate_adaptive.m scripts/ode/private/kahan.m scripts/ode/private/ode_event_handler.m scripts/ode/private/odedefaults.m scripts/ode/private/odemergeopts.m scripts/ode/private/runge_kutta_23.m scripts/ode/private/runge_kutta_45_dorpri.m scripts/ode/private/runge_kutta_interpolate.m scripts/ode/private/starting_stepsize.m |
| diffstat | 17 files changed, 358 insertions(+), 375 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ode/module.mk Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/module.mk Sat Oct 15 22:39:29 2016 -0700 @@ -3,7 +3,7 @@ scripts/ode/private scripts_ode_PRIVATE_FCN_FILES = \ - scripts/ode/private/AbsRel_Norm.m \ + scripts/ode/private/AbsRel_norm.m \ scripts/ode/private/integrate_adaptive.m \ scripts/ode/private/kahan.m \ scripts/ode/private/odedefaults.m \
--- a/scripts/ode/ode23.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/ode23.m Sat Oct 15 22:39:29 2016 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2016, Carlo de Falco -## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## Copyright (C) 2016 Carlo de Falco +## Copyright (C) 2016 Francesco Faccio <francesco.faccio@mail.polimi.it> ## Copyright (C) 2014-2016 Jacopo Corno <jacopo.corno@gmail.com> ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it> ## Copyright (C) 2006-2016 Thomas Treichl <treichl@users.sourceforge.net> @@ -86,9 +86,10 @@ ## [@var{t},@var{y}] = ode23 (fvdp, [0, 20], [2, 0]); ## @end group ## @end example -## @seealso{odeset, odeget} +## @seealso{odeset, odeget, ode45} ## @end deftypefn +## FIXME: We store ChangeLog information in Mercurial. Can this be deleted? ## ChangeLog: ## 20010703 the function file "ode23.m" was written by Marc Compere ## under the GPL for the use with this software. This function has been @@ -111,10 +112,10 @@ ## varargin{1:len} are parameters for fun odeopts = odeset (); funarguments = varargin; - elseif (length (varargin) > 1) + elseif (numel (varargin) > 1) ## varargin{1} is an ODE options structure opt odeopts = varargin{1}; - funarguments = {varargin{2:length(varargin)}}; + funarguments = {varargin{2:numel (varargin)}}; else # if (isstruct (varargin{1})) odeopts = varargin{1}; funarguments = {}; @@ -129,7 +130,7 @@ "ode23: TRANGE must be a numeric vector"); endif - if (length (trange) < 2) + if (numel (trange) < 2) error ("Octave:invalid-input-arg", "ode23: TRANGE must contain at least 2 elements"); elseif (trange(2) == trange(1)) @@ -158,18 +159,16 @@ "ode23: FUN must be a valid function handle"); endif - ## Start preprocessing, have a look which options are set in odeopts, - ## check if an invalid or unused option is set - - + ## check if an invalid or unused option is set. + ## FIXME: Why persistent? Won't these have different values for every + ## run of ode23? persistent defaults = []; persistent classes = []; persistent attributes = []; - - [defaults, classes, attributes] = odedefaults (numel (init), trange(1), - trange(end)); + [defaults, classes, attributes] = odedefaults (numel (init), + trange(1), trange(end)); defaults = rmfield (defaults, {"Jacobian", "JPattern", "Vectorized", ... "MvPattern", "MassSingular", ... @@ -181,7 +180,7 @@ "MvPattern", "MassSingular", ... "InitialSlope", "MaxOrder", "BDF"}); - odeopts = odemergeopts (odeopts, defaults, classes, attributes, 'ode23'); + odeopts = odemergeopts ("ode23", odeopts, defaults, classes, attributes); odeopts.funarguments = funarguments; odeopts.direction = direction; @@ -192,7 +191,7 @@ else odeopts.havenonnegative = false; warning ("Octave:invalid-input-arg", - ["ode23: option \"NonNegative\" is ignored", ... + ['ode23: option "NonNegative" is ignored', ... " when mass matrix is set\n"]); endif else @@ -208,14 +207,12 @@ if (isempty (odeopts.InitialStep)) odeopts.InitialStep = odeopts.direction * ... - starting_stepsize (order, fun, trange(1), - init, odeopts.AbsTol, - odeopts.RelTol, - strcmp (odeopts.NormControl, - "on"), odeopts.funarguments); + starting_stepsize (order, fun, trange(1), init, + odeopts.AbsTol, odeopts.RelTol, + strcmp (odeopts.NormControl, "on"), + odeopts.funarguments); endif - if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass)) havemasshandle = false; mass = odeopts.Mass; # constant mass @@ -225,33 +222,32 @@ havemasshandle = false; # mass = diag (ones (length (init), 1), 0); endif - ## Starting the initialization of the core solver ode23 if (havemasshandle) # Handle only the dynamic mass matrix, - if (! strcmp (odeopts.MStateDependence, "none")) # constant mass matrices have already + if (! strcmp (odeopts.MStateDependence, "none")) + ## FIXME: How is this comment supposed to end? + ## 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 ((! strcmp (odeopts.MStateDependence, "none")) == false) + \ fun (t, x, odeopts.funarguments{:}); + else mass = @(t) odeopts.Mass (t, odeopts.funarguments{:}); fun = @(t,x) mass (t, odeopts.funarguments{:}) ... - \ fun (t, x, odeopts.funarguments{:}); + \ fun (t, x, odeopts.funarguments{:}); endif endif - - solution = integrate_adaptive (@runge_kutta_23, ... + solution = integrate_adaptive (@runge_kutta_23, order, fun, trange, init, odeopts); - ## Postprocessing, do whatever when terminating integration algorithm if (odeopts.haveoutputfunction) # Cleanup plotter - feval (odeopts.OutputFcn, solution.t(end), ... + feval (odeopts.OutputFcn, solution.t(end), solution.x(end,:)', "done", odeopts.funarguments{:}); endif if (! isempty (odeopts.Events)) # Cleanup event function handling - ode_event_handler (odeopts.Events, solution.t(end), ... + ode_event_handler (odeopts.Events, solution.t(end), solution.x(end,:)', "done", odeopts.funarguments{:}); endif @@ -307,7 +303,6 @@ %!demo -%! %! ## Demonstrate convergence order for ode23 %! tol = 1e-5 ./ 10.^[0:8]; %! for i = 1 : numel (tol) @@ -336,45 +331,45 @@ ## for this function. ## For further tests we also define a reference solution (computed at high ## accuracy) -%!function ydot = fpol (t, y) # The Van der Pol -%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%!function ydot = fpol (t, y) # The Van der Pol ODE +%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; %!endfunction %!function ref = fref () # The computed reference sol -%! ref = [0.32331666704577, -1.83297456798624]; +%! ref = [0.32331666704577, -1.83297456798624]; %!endfunction %!function jac = fjac (t, y, varargin) # its Jacobian -%! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]; +%! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]; %!endfunction %!function jac = fjcc (t, y, varargin) # sparse type -%! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]); +%! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]); %!endfunction %!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 +%! val = fpol (t, y, varargin); # We use the derivatives +%! trm = zeros (2,1); # that's why component 2 +%! dir = ones (2,1); # does not seem to be exact %!endfunction %!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 +%! val = fpol (t, y, varargin); # We use the derivatives +%! trm = ones (2,1); # that's why component 2 +%! dir = ones (2,1); # does not seem to be exact %!endfunction %!function mas = fmas (t, y, varargin) -%! mas = [1, 0; 0, 1]; # Dummy mass matrix for tests +%! mas = [1, 0; 0, 1]; # Dummy mass matrix for tests %!endfunction %!function mas = fmsa (t, y, varargin) -%! mas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix +%! mas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix %!endfunction %!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 flag"); -%! endif +%! 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 flag'); +%! endif %!endfunction %! %!test # two output arguments @@ -487,27 +482,21 @@ %! opt = odeset ("Mass", @fmas, "MStateDependence", "strong"); %! sol = ode23 (@fpol, [0 2], [2 0], opt); %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); -%! -%! ## test for MvPattern option is missing -%! ## test for InitialSlope option is missing -%! ## test for MaxOrder option is missing + +## FIXME: Missing tests. +## test for MvPattern option is missing +## test for InitialSlope option is missing +## test for MaxOrder option is missing ## Test input validation %!error ode23 () %!error ode23 (1) %!error ode23 (1,2) -%!error <TRANGE must be a numeric> -%! ode23 (@fpol, {[0 25]}, [3 15 1]); -%!error <TRANGE must be a .* vector> -%! ode23 (@fpol, [0 25; 25 0], [3 15 1]); -%!error <TRANGE must contain at least 2 elements> -%! ode23 (@fpol, [1], [3 15 1]); -%!error <invalid time span> -%! ode23 (@fpol, [1 1], [3 15 1]); -%!error <INIT must be a numeric> -%! ode23 (@fpol, [0 25], {[3 15 1]}); -%!error <INIT must be a .* vector> -%! ode23 (@fpol, [0 25], [3 15 1; 3 15 1]); -%!error <FUN must be a valid function handle> -%! ode23 (1, [0 25], [3 15 1]); +%!error <TRANGE must be a numeric> ode23 (@fpol, {[0 25]}, [3 15 1]) +%!error <TRANGE must be a .* vector> ode23 (@fpol, [0 25; 25 0], [3 15 1]) +%!error <TRANGE must contain at least 2 elements> ode23 (@fpol, [1], [3 15 1]) +%!error <invalid time span> ode23 (@fpol, [1 1], [3 15 1]) +%!error <INIT must be a numeric> ode23 (@fpol, [0 25], {[3 15 1]}) +%!error <INIT must be a .* vector> ode23 (@fpol, [0 25], [3 15 1; 3 15 1]) +%!error <FUN must be a valid function handle> ode23 (1, [0 25], [3 15 1])
--- a/scripts/ode/ode45.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/ode45.m Sat Oct 15 22:39:29 2016 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2016, Carlo de Falco -## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## Copyright (C) 2016 Carlo de Falco +## Copyright (C) 2016 Francesco Faccio <francesco.faccio@mail.polimi.it> ## Copyright (C) 2014-2016 Jacopo Corno <jacopo.corno@gmail.com> ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it> ## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net> @@ -77,7 +77,7 @@ ## [@var{t},@var{y}] = ode45 (fvdp, [0, 20], [2, 0]); ## @end group ## @end example -## @seealso{odeset, odeget} +## @seealso{odeset, odeget, ode23} ## @end deftypefn function varargout = ode45 (fun, trange, init, varargin) @@ -94,10 +94,10 @@ ## varargin{1:len} are parameters for fun odeopts = odeset (); funarguments = varargin; - elseif (length (varargin) > 1) + elseif (numel (varargin) > 1) ## varargin{1} is an ODE options structure opt odeopts = varargin{1}; - funarguments = {varargin{2:length(varargin)}}; + funarguments = {varargin{2:numel (varargin)}}; else # if (isstruct (varargin{1})) odeopts = varargin{1}; funarguments = {}; @@ -112,8 +112,7 @@ "ode45: TRANGE must be a numeric vector"); endif - - if (length (trange) < 2) + if (numel (trange) < 2) error ("Octave:invalid-input-arg", "ode45: TRANGE must contain at least 2 elements"); elseif (trange(1) == trange(2)) @@ -142,19 +141,17 @@ "ode45: FUN must be a valid function handle"); endif - ## Start preprocessing, have a look which options are set in odeopts, ## check if an invalid or unused option is set - + ## FIXME: Why persistent when it is changed with every run of ode45? persistent defaults = []; persistent classes = []; persistent attributes = []; + [defaults, classes, attributes] = odedefaults (numel (init), + trange(1), trange(end)); - [defaults, classes, attributes] = odedefaults (numel (init), trange(1), - trange(end)); - - defaults = odeset (defaults, 'Refine', 4); + defaults = odeset (defaults, "Refine", 4); defaults = rmfield (defaults, {"Jacobian", "JPattern", "Vectorized", ... "MvPattern", "MassSingular", ... "InitialSlope", "MaxOrder", "BDF"}); @@ -165,7 +162,7 @@ "MvPattern", "MassSingular", ... "InitialSlope", "MaxOrder", "BDF"}); - odeopts = odemergeopts (odeopts, defaults, classes, attributes, 'ode45'); + odeopts = odemergeopts ("ode45", odeopts, defaults, classes, attributes); odeopts.funarguments = funarguments; odeopts.direction = direction; @@ -176,7 +173,7 @@ else odeopts.havenonnegative = false; warning ("Octave:invalid-input-arg", - ["ode45: option 'NonNegative' is ignored", ... + ['ode45: option "NonNegative" is ignored', ... " when mass matrix is set\n"]); endif else @@ -192,14 +189,12 @@ if (isempty (odeopts.InitialStep)) odeopts.InitialStep = odeopts.direction * ... - starting_stepsize (order, fun, trange(1), - init, odeopts.AbsTol, - odeopts.RelTol, - strcmp (odeopts.NormControl, - "on"), odeopts.funarguments); + starting_stepsize (order, fun, trange(1), init, + odeopts.AbsTol, odeopts.RelTol, + strcmp (odeopts.NormControl, "on"), + odeopts.funarguments); endif - if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass)) havemasshandle = false; mass = odeopts.Mass; # constant mass @@ -209,33 +204,31 @@ havemasshandle = false; # mass = diag (ones (length (init), 1), 0); endif - ## Starting the initialization of the core solver ode45 if (havemasshandle) # Handle only the dynamic mass matrix, - if (! strcmp (odeopts.MStateDependence, "none")) # constant mass matrices have already + if (! strcmp (odeopts.MStateDependence, "none")) + ### 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 ((! strcmp (odeopts.MStateDependence, "none")) == false) + \ fun (t, x, odeopts.funarguments{:}); + else mass = @(t) odeopts.Mass (t, odeopts.funarguments{:}); fun = @(t,x) mass (t, odeopts.funarguments{:}) ... - \ fun (t, x, odeopts.funarguments{:}); + \ fun (t, x, odeopts.funarguments{:}); endif endif - solution = integrate_adaptive (@runge_kutta_45_dorpri, order, fun, trange, init, odeopts); - ## Postprocessing, do whatever when terminating integration algorithm if (odeopts.haveoutputfunction) # Cleanup plotter - feval (odeopts.OutputFcn, solution.t(end), ... + feval (odeopts.OutputFcn, solution.t(end), solution.x(end,:)', "done", odeopts.funarguments{:}); endif if (! isempty (odeopts.Events)) # Cleanup event function handling - ode_event_handler (odeopts.Events, solution.t(end), ... + ode_event_handler (odeopts.Events, solution.t(end), solution.x(end,:)', "done", odeopts.funarguments{:}); endif @@ -291,7 +284,6 @@ %!demo -%! %! ## Demonstrate convergence order for ode45 %! tol = 1e-5 ./ 10.^[0:8]; %! for i = 1 : numel (tol) @@ -308,7 +300,7 @@ %! loglog (h, tol, "-ob", %! h, err, "-b", %! h, (h/h(end)) .^ 4 .* tol(end), "k--", -%! h, (h/h(end)) .^ 5 .* tol(end), "k-") +%! h, (h/h(end)) .^ 5 .* tol(end), "k-"); %! axis tight %! xlabel ("h"); %! ylabel ("err(h)"); @@ -316,49 +308,49 @@ %! legend ("imposed tolerance", "ode45 (relative) error", %! "order 4", "order 5", "location", "northwest"); -## We are using the "Van der Pol" implementation for all tests that are done +## We are using the Van der Pol equation for all tests that are done ## for this function. ## For further tests we also define a reference solution (computed at high ## accuracy) -%!function ydot = fpol (t, y) # The Van der Pol -%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%!function ydot = fpol (t, y) # The Van der Pol ODE +%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; %!endfunction %!function ref = fref () # The computed reference solution -%! ref = [0.32331666704577, -1.83297456798624]; +%! ref = [0.32331666704577, -1.83297456798624]; %!endfunction %!function jac = fjac (t, y, varargin) # its Jacobian -%! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]; +%! jac = [0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]; %!endfunction %!function jac = fjcc (t, y, varargin) # sparse type -%! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]); +%! jac = sparse ([0, 1; -1 - 2 * y(1) * y(2), 1 - y(1)^2]); %!endfunction %!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 +%! val = fpol (t, y, varargin); # We use the derivatives +%! trm = zeros (2,1); # that's why component 2 +%! dir = ones (2,1); # does not seem to be exact %!endfunction %!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 +%! val = fpol (t, y, varargin); # We use the derivatives +%! trm = ones (2,1); # that's why component 2 +%! dir = ones (2,1); # does not seem to be exact %!endfunction %!function mas = fmas (t, y, varargin) -%! mas = [1, 0; 0, 1]; # Dummy mass matrix for tests +%! mas = [1, 0; 0, 1]; # Dummy mass matrix for tests %!endfunction %!function mas = fmsa (t, y, varargin) -%! mas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix +%! mas = sparse ([1, 0; 0, 1]); # A sparse dummy matrix %!endfunction %!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 flag'); -%! endif +%! 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 flag'); +%! endif %!endfunction %! %!test # two output arguments @@ -489,21 +481,19 @@ %! sol = ode45 (@fpol, [0 2], [2 0], opt); %! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +## FIXME: Missing tests. +## test for MvPattern option is missing +## test for InitialSlope option is missing +## test for MaxOrder option is missing + %!error ode45 () %!error ode45 (1) %!error ode45 (1,2) -%!error <TRANGE must be a numeric> -%! ode45 (@fpol, {[0 25]}, [3 15 1]); -%!error <TRANGE must be a .* vector> -%! ode45 (@fpol, [0 25; 25 0], [3 15 1]); -%!error <TRANGE must contain at least 2 elements> -%! ode45 (@fpol, [1], [3 15 1]); -%!error <invalid time span> -%! ode45 (@fpol, [1 1], [3 15 1]); -%!error <INIT must be a numeric> -%! ode45 (@fpol, [0 25], {[3 15 1]}); -%!error <INIT must be a .* vector> -%! ode45 (@fpol, [0 25], [3 15 1; 3 15 1]); -%!error <FUN must be a valid function handle> -%! ode45 (1, [0 25], [3 15 1]); +%!error <TRANGE must be a numeric> ode45 (@fpol, {[0 25]}, [3 15 1]) +%!error <TRANGE must be a .* vector> ode45 (@fpol, [0 25; 25 0], [3 15 1]) +%!error <TRANGE must contain at least 2 elements> ode45 (@fpol, [1], [3 15 1]) +%!error <invalid time span> ode45 (@fpol, [1 1], [3 15 1]) +%!error <INIT must be a numeric> ode45 (@fpol, [0 25], {[3 15 1]}) +%!error <INIT must be a .* vector> ode45 (@fpol, [0 25], [3 15 1; 3 15 1]) +%!error <FUN must be a valid function handle> ode45 (1, [0 25], [3 15 1])
--- a/scripts/ode/odeget.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/odeget.m Sat Oct 15 22:39:29 2016 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2016, Carlo de Falco -## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## Copyright (C) 2016 Carlo de Falco +## Copyright (C) 2016 Francesco Faccio <francesco.faccio@mail.polimi.it> ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it> ## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net> ## @@ -41,19 +41,19 @@ function val = odeget (ode_opt, field, default = [], opt = "") - validateattributes (ode_opt, {'struct'}, {'nonempty'}); - validateattributes (field, {'char'}, {'nonempty'}); - + validateattributes (ode_opt, {"struct"}, {"nonempty"}); + validateattributes (field, {"char"}, {"nonempty"}); + if (! isfield (ode_opt, field)) - error ('Octave:odeget:InvalidPropName', + error ("Octave:odeget:InvalidPropName", 'odeget: Unrecognized property name "%s".', field); else val = ode_opt.(field); - if (isempty (val)) + if (isempty (val)) val = default; endif endif - + endfunction @@ -73,12 +73,12 @@ %!assert (odeget (odeset (), "Mass"), []) %!assert (odeget (odeset (), "AbsTol", 1e-9), 1e-9) %!assert (odeget (odeset ("AbsTol", 1e-9), "AbsTol", []), 1e-9) -%!assert (odeget (odeset ('foo', 42), 'foo'), 42) +%!assert (odeget (odeset ("foo", 42), "foo"), 42) %!error odeget () %!error odeget (1) %!error odeget (1,2,3,4,5) %!error odeget (1, "opt1") %!error odeget (struct ("opt1", 1), 1) -%!error odeget (struct ("opt1", 1), "foo"); +%!error odeget (struct ("opt1", 1), "foo")
--- a/scripts/ode/odeplot.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/odeplot.m Sat Oct 15 22:39:29 2016 -0700 @@ -19,33 +19,34 @@ ## Author: Thomas Treichl <treichl@users.sourceforge.net> ## -*- texinfo -*- -## @deftypefn {} {[@var{ret}] =} odeplot (@var{t}, @var{y}, @var{flag}) +## @deftypefn {} {@var{ret} =} odeplot (@var{t}, @var{y}, @var{flag}) ## -## Open a new figure window and plot the results from the variable @var{y} of -## type column vector over time while solving. The types and the values of -## the input parameter @var{t} and the output parameter @var{ret} depend on -## the input value @var{flag} that is of type string. If @var{flag} is +## Open a new figure window and plot input @var{y} over time during the +## solving of an ode problem. +## +## The input @var{y} is a column vector. The types and values of the input +## parameter @var{t} and output parameter @var{ret} depend on the input +## @var{flag} that is of type string. If @var{flag} is ## ## @table @option -## @item @qcode{"init"} -## then @var{t} must be a double column vector of length 2 with the first and -## the last time step and nothing is returned from this function, +## @item @qcode{"init"} +## then @var{t} must be a column vector of length 2 with the first and +## the last time step; Nothing is returned from this function. ## -## @item @qcode{""} -## then @var{t} must be a double scalar specifying the actual time step and -## the return value is false (resp. value 0) for @qcode{"not stop solving"}, +## @item @qcode{""} +## then @var{t} must be a scalar double specifying the actual time step; +## The return value is false (resp. value 0) for @qcode{"not stop solving"}. ## -## @item @qcode{"done"} -## then @var{t} must be a double scalar specifying the last time step and -## nothing is returned from this function. +## @item @qcode{"done"} +## then @var{t} must be a scalar double specifying the last time step; +## Nothing is returned from this function. ## @end table ## ## This function is called by an ode solver function if it was specified in -## an options structure with the @code{odeset}. This function is an -## 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. +## an options structure with @code{odeset}. This function is an internal +## helper function; It should never be necessary for this function to be +## directly called by a user. There is only minimal error detection +## implemented in order to to achieve the highest performance. ## ## For example, solve an anonymous implementation of the ## @nospell{@qcode{"Van der Pol"}} equation and display the results while @@ -59,34 +60,38 @@ ## sol = ode45 (fvdb, [0 20], [2 0], opt); ## @end group ## @end example -## @end deftypefn ## ## @seealso{odeset, odeget} +## @end deftypefn function ret = odeplot (t, y, flag, varargin) ## No input argument check is done for a higher processing speed - persistent fig; persistent told; - persistent yold; persistent counter; + persistent fig told yold counter; if (strcmp (flag, "init")) ## Nothing to return, t is either the time slot [tstart tstop] - ## or [t0, t1, ..., tn], y is the inital value vector "init" - fig = figure; told = t(1,1); yold = y(:,1); + ## or [t0, t1, ..., tn], y is the initial value vector "init" counter = 1; + fig = figure (); + told = t(1,1); + yold = y(:,1); elseif (isempty (flag)) ## Return something, either false for "not stopping ## the integration" or true for "stopping the integration" - counter = counter + 1; figure (fig); + counter += 1; + figure (fig); told(counter,1) = t(1,1); yold(:,counter) = y(:,1); + ## FIXME: Why not use '.' rather than 'o' and skip the markersize? + ## FIXME: Why not just update the xdata, ydata properties? + ## Calling plot involves a lot of overhead. plot (told, yold, "-o", "markersize", 1); drawnow; ret = false; elseif (strcmp (flag, "done")) - ## Cleanup has to be done, clear the persistent variables because - ## we don't need them anymore + ## Cleanup after ode solver has finished. clear ("figure", "told", "yold", "counter"); endif @@ -95,8 +100,9 @@ %!demo -%! % solve an anonymous implementation of the "Van der Pol" equation -%! % and display the results while solving +%! # Solve an anonymous implementation of the Van der Pol equation +%! # and display the results while solving %! fvdb = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; %! opt = odeset ("OutputFcn", @odeplot, "RelTol", 1e-6); %! sol = ode45 (fvdb, [0 20], [2 0], opt); +
--- a/scripts/ode/odeset.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/odeset.m Sat Oct 15 22:39:29 2016 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2016, Carlo de Falco -## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## Copyright (C) 2016 Carlo de Falco +## Copyright (C) 2016 Francesco Faccio <francesco.faccio@mail.polimi.it> ## Copyright (C) 2013-2016 Roberto Porcu' <roberto.porcu@polimi.it> ## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net> ## @@ -20,15 +20,18 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {} {} odeset () +## @deftypefn {} {@var{odestruct} =} odeset () ## @deftypefnx {} {@var{odestruct} =} odeset (@var{"field1"}, @var{value1}, @var{"field2"}, @var{value2}, @dots{}) ## @deftypefnx {} {@var{odestruct} =} odeset (@var{oldstruct}, @var{"field1"}, @var{value1}, @var{"field2"}, @var{value2}, @dots{}) ## @deftypefnx {} {@var{odestruct} =} odeset (@var{oldstruct}, @var{newstruct}) +## @deftypefnx {} {} odeset () ## ## Create or modify an ODE options structure. ## -## When called without an input argument, return a new ODE options structure -## that contains all possible fields initialized to their default values. +## When called with no input argument and one output argument, return a new ODE +## options structure that contains all possible fields initialized to their +## default values. If no output argument is requested, display a list of +## the common ODE solver options along with their default value. ## ## If called with string input arguments @var{"field1"}, @var{"field2"}, ## @dots{} identifying valid ODE options then return a new ODE options @@ -44,6 +47,7 @@ ## @var{newstruct} overwrite all values from the structure @var{oldstruct} with ## new values from the structure @var{newstruct}. Empty values in ## @var{newstruct} will not overwrite values in @var{oldstruct}. +## ## @seealso{odeget} ## @end deftypefn @@ -52,7 +56,6 @@ persistent p; if (isempty (p)) - p = inputParser (); p.addParameter ("AbsTol", []); p.addParameter ("BDF", []); @@ -63,6 +66,7 @@ p.addParameter ("JConstant", []); p.addParameter ("JPattern", []); p.addParameter ("Mass", []); + ## FIXME: MassConstant does not appear in Matlab documentation for odeset p.addParameter ("MassConstant", []); p.addParameter ("MassSingular", []); p.addParameter ("MaxOrder", []); @@ -78,7 +82,6 @@ p.addParameter ("Stats", []); p.addParameter ("Vectorized", []); p.KeepUnmatched = true; - endif if (nargin == 0 && nargout == 0) @@ -88,20 +91,24 @@ odestruct = p.Results; odestruct_extra = p.Unmatched; - s1 = cellfun (@(x) ifelse (iscell(x), {x}, x), - struct2cell(odestruct), - 'UniformOutput', false); + ## FIXME: For speed, shouldn't this merge of structures only occur + ## when there is something in odestruct_extra? + ## FIXME: Should alphabetical order of fieldnames be maintained + ## by using sort? + s1 = cellfun (@(x) ifelse (iscell (x), {x}, x), + struct2cell (odestruct), + "UniformOutput", false); - s2 = cellfun (@(x) ifelse (iscell(x), {x}, x), - struct2cell(odestruct_extra), - 'UniformOutput', false); - + s2 = cellfun (@(x) ifelse (iscell (x), {x}, x), + struct2cell (odestruct_extra), + "UniformOutput", false); + C = [fieldnames(odestruct) s1; fieldnames(odestruct_extra) s2]; - + odestruct = struct (C'{:}); endif - + endfunction ## function to print all possible options @@ -110,29 +117,29 @@ disp ("List of the most common ODE solver options."); disp ("Default values are in square brackets."); disp (""); - disp (" AbsTol: scalar or vector, >0, [1e-6]"); - disp (" BDF: binary, {'on', ['off']}"); - disp (" Events: function_handle, []"); - disp (" InitialSlope: vector, []"); - disp (" InitialStep: scalar, >0, []"); - disp (" Jacobian: matrix or function_handle, []"); - disp (" JConstant: binary, {'on', ['off']}"); - disp (" JPattern: sparse matrix, []"); - disp (" Mass: matrix or function_handle, []"); - disp (" MassConstant: binary, {'on', ['off']}"); - disp (" MassSingular: switch, {'yes', ['maybe'], 'no'}"); - disp (" MaxOrder: switch, {1, 2, 3, 4, [5]}"); - disp (" MaxStep: scalar, >0, []"); - disp (" MStateDependence: switch, {'none', ['weak'], 'strong'}"); - disp (" MvPattern: sparse matrix, []"); - disp (" NonNegative: vector of integers, []"); - disp (" NormControl: binary, {'on', ['off']}"); - disp (" OutputFcn: function_handle, []"); - disp (" OutputSel: scalar or vector, []"); - disp (" Refine: scalar, integer, >0, []"); - disp (" RelTol: scalar, >0, [1e-3]"); - disp (" Stats: binary, {'on', ['off']}"); - disp (" Vectorized: binary, {'on', ['off']}"); + disp (' AbsTol: scalar or vector, >0, [1e-6]'); + disp (' BDF: binary, {["off"], "on"}'); + disp (' Events: function_handle, []'); + disp (' InitialSlope: vector, []'); + disp (' InitialStep: scalar, >0, []'); + disp (' Jacobian: matrix or function_handle, []'); + disp (' JConstant: binary, {["off"], "on"}'); + disp (' JPattern: sparse matrix, []'); + disp (' Mass: matrix or function_handle, []'); + disp (' MassConstant: binary, {["off"], "on"}'); + disp (' MassSingular: switch, {["maybe"], "no", "yes"}'); + disp (' MaxOrder: switch, {[5], 1, 2, 3, 4, }'); + disp (' MaxStep: scalar, >0, []'); + disp (' MStateDependence: switch, {["weak"], "none", "strong"}'); + disp (' MvPattern: sparse matrix, []'); + disp (' NonNegative: vector of integers, []'); + disp (' NormControl: binary, {["off"], "on"}'); + disp (' OutputFcn: function_handle, []'); + disp (' OutputSel: scalar or vector, []'); + disp (' Refine: scalar, integer, >0, []'); + disp (' RelTol: scalar, >0, [1e-3]'); + disp (' Stats: binary, {["off"], "on"}'); + disp (' Vectorized: binary, {["off"], "on"}'); endfunction @@ -192,8 +199,7 @@ %!error <argument 'OPT1' is not a valid parameter> odeset (odeset (), "opt1") %!error odeset (odeset (), 1, 1) -##FIXME: Add not exact match option +## FIXME: Add inexact match option ## %!warning <no exact match for 'Rel'. Assuming 'RelTol'> odeset ("Rel", 1); ## %!error <Possible fields found: InitialSlope, InitialStep> odeset ("Initial", 1) -
--- a/scripts/ode/private/AbsRel_Norm.m Sat Oct 15 10:30:48 2016 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,37 +0,0 @@ -## Copyright (C) 2014-2016 Jacopo Corno <jacopo.corno@gmail.com> -## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or (at -## your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -## General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## <http://www.gnu.org/licenses/>. - -## -*- texinfo -*- -## @deftypefn {} {retval =} AbsRel_Norm (@var{x}, @var{x_old}, @var{AbsTol}, @var{RelTol}, @var{normcoontrol}, @var{y}) -## Undocumented internal function. -## @end deftypefn - -function retval = AbsRel_Norm (x, x_old, AbsTol, RelTol, normcontrol, y = zeros (size (x))) - - n = numel (x); - - sc = AbsTol + max (abs (x), abs (x_old)) .* RelTol; - if (normcontrol) - retval = max (abs (x - y) ./ sc); - else - retval = sqrt ((1 / n) * sumsq ((x - y) ./ sc)); - endif - -endfunction -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/AbsRel_norm.m Sat Oct 15 22:39:29 2016 -0700 @@ -0,0 +1,37 @@ +## Copyright (C) 2014-2016 Jacopo Corno <jacopo.corno@gmail.com> +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or (at +## your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## <http://www.gnu.org/licenses/>. + +## -*- texinfo -*- +## @deftypefn {} {retval =} AbsRel_norm (@var{x}, @var{x_old}, @var{AbsTol}, @var{RelTol}, @var{normcoontrol}, @var{y}) +## Undocumented internal function. +## @end deftypefn + +function retval = AbsRel_norm (x, x_old, AbsTol, RelTol, normcontrol, y = zeros (size (x))) + + n = numel (x); + + sc = AbsTol + max (abs (x), abs (x_old)) .* RelTol; + if (normcontrol) + retval = max (abs (x - y) ./ sc); + else + retval = sqrt ((1 / n) * sumsq ((x - y) ./ sc)); + endif + +endfunction +
--- a/scripts/ode/private/integrate_adaptive.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/integrate_adaptive.m Sat Oct 15 22:39:29 2016 -0700 @@ -21,10 +21,10 @@ ## @deftypefn {} {@var{solution} =} integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@func}, @var{tspan}, @var{x0}, @var{options}) ## ## This function file can be called by an ODE solver function in order to -## integrate the set of ODEs on the interval @var{[t0, t1]} with an -## adaptive timestep. +## integrate the set of ODEs on the interval @var{[t0, t1]} with an adaptive +## timestep. ## -## The function returns a structure @var{solution} with two fieldss: @var{t} +## The function returns a structure @var{solution} with two fields: @var{t} ## and @var{y}. @var{t} is a column vector and contains the time stamps. ## @var{y} is a matrix in which each column refers to a different unknown ## of the problem and the row number is the same as the @var{t} row number. @@ -61,8 +61,6 @@ ## stepper. ## ## @end deftypefn -## -## @seealso{ode45, ode23} function solution = integrate_adaptive (stepper, order, func, tspan, x0, options) @@ -75,8 +73,10 @@ ## Get first initial timestep dt = options.InitialStep; if (isempty (dt)) - dt = starting_stepsize (order, func, t, x, options.AbsTol, options.RelTol, - strcmp (options.NormControl, "on"), options.funarguments); + dt = starting_stepsize (order, func, t, x, + options.AbsTol, options.RelTol, + strcmp (options.NormControl, "on"), + options.funarguments); endif dir = options.direction; @@ -118,13 +118,13 @@ k_vals = []; iout = istep = 1; - + while (dir * t_old < dir * tspan(end)) ## Compute integration step from t_old to t_new = t_old + dt [t_new, options.comp] = kahan (t_old, options.comp, dt); [t_new, x_new, x_est, new_k_vals] = ... - stepper (func, t_old, x_old, dt, options, k_vals, t_new); + stepper (func, t_old, x_old, dt, options, k_vals, t_new); solution.cntcycles += 1; @@ -133,7 +133,8 @@ x_est(nn, end) = abs (x_est(nn, end)); endif - err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol, + ## FIXME: Take strcmp out of while loop and calculate just once + err = AbsRel_norm (x_new, x_old, options.AbsTol, options.RelTol, strcmp (options.NormControl, "on"), x_est); ## Accept solution only if err <= 1.0 @@ -153,9 +154,9 @@ t(t_caught) = tspan(t_caught); iout = max (t_caught); x(:, t_caught) = ... - runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], - t(t_caught), new_k_vals, dt, func, - options.funarguments); + runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], + t(t_caught), new_k_vals, dt, func, + options.funarguments); istep += 1; @@ -167,8 +168,8 @@ id = t_caught(idenseout); td = t(id); solution.event = ... - ode_event_handler (options.Events, t(id), x(:, id), [], - options.funarguments{:}); + 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); @@ -191,7 +192,7 @@ 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') .'; + approxtime, "linear") .'; if (! isempty (options.OutputSel)) approxvals = approxvals(options.OutputSel, :); endif @@ -208,7 +209,7 @@ endif - else + else # not fixed times t(++istep) = t_new; x(:, istep) = x_new; @@ -218,8 +219,8 @@ ## Stop integration if eventbreak is true. if (! isempty (options.Events)) solution.event = ... - ode_event_handler (options.Events, t(istep), x(:, istep), [], - options.funarguments{:}); + 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, :).'; @@ -235,7 +236,7 @@ approxtime = linspace (t_old, t_new, cnt); approxvals = interp1 ([t_old, t_new], [x_old, x_new] .', - approxtime, 'linear') .'; + approxtime, "linear") .'; if (! isempty (options.OutputSel)) approxvals = approxvals(options.OutputSel, :); endif @@ -256,12 +257,12 @@ x_old = x_new; k_vals = new_k_vals; - else + else # error condition ireject += 1; - ## Stop solving because, in the last 5,000 steps, no successful valid - ## value has been found + ## Stop solving if, in the last 5,000 steps, no successful valid + ## value has been found. if (ireject >= 5_000) error (["integrate_adaptive: Solving was not successful. ", ... " The iterative integration loop exited at time", ... @@ -280,10 +281,10 @@ dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1)))); dt = dir * min (abs (dt), options.MaxStep); if (! (abs (dt) > eps (t (end)))) - break + break; endif - - ## make sure we don't go past tpan(end) + + ## Make sure we don't go past tpan(end) dt = dir * min (abs (dt), abs (tspan(end) - t_old)); endwhile
--- a/scripts/ode/private/kahan.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/kahan.m Sat Oct 15 22:39:29 2016 -0700 @@ -27,14 +27,14 @@ ## 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 @code{integrate_adaptive} and -## @code{integrate_const} to better catch equality comparisons. +## This function is called by @code{integrate_adaptive} to better catch +## equality comparisons. ## -## The first input argument is the variable that will contain the summation, -## so that is also returned as first output argument in order to reuse it in -## next calls to @code{kahan} function. +## The first input argument is the variable that will contain the summation. +## This variable is also returned as the first output argument in order to +## reuse it in subsequent calls to @code{kahan} function. ## -## The second input argument contains the compensation term and it is returned +## The second input argument contains the compensation term and is returned ## as the second output argument so that it can be reused in future calls of ## the same summation. ##
--- a/scripts/ode/private/ode_event_handler.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/ode_event_handler.m Sat Oct 15 22:39:29 2016 -0700 @@ -17,7 +17,7 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {} {@var{sol} =} ode_event_handler (@var{@@evtfun}, @var{t}, @var{y}, @var{flag}, @var{par1}, @var{par2}, @dots{}) +## @deftypefn {} {@var{retval} =} 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. @@ -36,7 +36,7 @@ ## ## @item @qcode{"calc"} ## then do the evaluation of the event function and return the solution -## @var{sol} as type cell array of size 4, +## @var{retval} as type cell array of size 4, ## ## @item @qcode{"done"} ## then cleanup internal variables of the function @@ -52,8 +52,6 @@ ## 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) @@ -72,7 +70,6 @@ ## 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
--- a/scripts/ode/private/odedefaults.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/odedefaults.m Sat Oct 15 22:39:29 2016 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2016, Carlo de Falco -## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## Copyright (C) 2016 Carlo de Falco +## Copyright (C) 2016 Francesco Faccio <francesco.faccio@mail.polimi.it> ## ## This file is part of Octave. ## @@ -18,7 +18,7 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {} {[@var{d}, @var{c}, @var{a}] =} odedefaults (@var{n}, @var{t0}, @var{tf}) +## @deftypefn {} {[@var{defaults}, @var{classes}, @var{attributes}] =} odedefaults (@var{n}, @var{t0}, @var{tf}) ## Undocumented internal function. ## @end deftypefn @@ -27,7 +27,7 @@ persistent defaults = struct ("AbsTol", 1e-6, "BDF", "off", "Events", [], - "InitialSlope", zeros(n,1), + "InitialSlope", zeros (n,1), "InitialStep", [], "Jacobian", [], "JConstant", "off", @@ -43,13 +43,13 @@ "NormControl", "off", "OutputFcn", [], "OutputSel", [], - "Refine", 1, + "Refine", 1, "RelTol", 1e-3, "Stats", "off", "Vectorized", "off"); - + defaults.MaxStep = (0.1 * abs (t0-tf)); - + persistent classes = struct ("AbsTol", {{"float"}}, "BDF", "char", "Events", {{"function_handle"}}, @@ -73,7 +73,7 @@ "RelTol", {{"float"}}, "Stats", "char", "Vectorized", "char"); - + persistent attributes = struct ("AbsTol", {{"real", "vector", "positive"}}, "BDF", {{"on", "off"}}, "Events", {{}}, @@ -99,3 +99,4 @@ "Stats", {{"on", "off"}}, "Vectorized", {{"on", "off"}}); endfunction +
--- a/scripts/ode/private/odemergeopts.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/odemergeopts.m Sat Oct 15 22:39:29 2016 -0700 @@ -1,4 +1,4 @@ -## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> +## Copyright (C) 2016 Francesco Faccio <francesco.faccio@mail.polimi.it> ## ## This file is part of Octave. ## @@ -16,8 +16,8 @@ ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. -function options = odemergeopts (useroptions, options, classes, - attributes, fun_name); +function options = odemergeopts (caller, useroptions, options, classes, + attributes); for [value, key] = options; @@ -25,18 +25,21 @@ if (! strcmp (classes.(key), "char")) validateattributes (useroptions.(key), classes.(key), - attributes.(key), fun_name, key); + attributes.(key), caller, key); elseif (ischar (useroptions.(key))) - validatestring (useroptions.(key), attributes.(key), fun_name, key); + validatestring (useroptions.(key), attributes.(key), caller, key); else error ("Octave:invalid-input-arg", - [fun_name ": invalid value assigned to field '%s'"], key); + [caller ": invalid value assigned to field '%s'"], key); endif options.(key) = useroptions.(key); endif + endfor + endfunction +
--- a/scripts/ode/private/runge_kutta_23.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/runge_kutta_23.m Sat Oct 15 22:39:29 2016 -0700 @@ -65,12 +65,12 @@ k_vals = [], t_next = t + dt) - persistent a = [0 0 0; - 1/2 0 0; - 0 3/4 0]; - persistent b = [0 1/2 3/4 1]; - persistent c = [(2/9) (1/3) (4/9)]; - persistent c_prime = [(7/24) (1/4) (1/3) (1/8)]; + persistent a = [0 0 0; + 1/2 0 0; + 0 3/4 0]; + persistent b = [0, 1/2, 3/4, 1]; + persistent c = [2/9, 1/3, 4/9]; + persistent c_prime = [7/24, 1/4, 1/3, 1/8]; s = t + dt * b; cc = dt * c; @@ -92,9 +92,9 @@ k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).', args{:}); k(:,3) = feval (f, s(3), x + k(:,2) * aa(3, 2).', args{:}); - ## compute new time and new values for the unkwnowns + ## compute new time and new values for the unknowns ## t_next = t + dt; - x_next = x + k(:,1:3) * cc(:); # 3rd order approximation + x_next = x + k(:,1:3) * cc(:); # 3rd order approximation ## if the estimation of the error is required if (nargout >= 3) @@ -105,3 +105,4 @@ endif endfunction +
--- a/scripts/ode/private/runge_kutta_45_dorpri.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/runge_kutta_45_dorpri.m Sat Oct 15 22:39:29 2016 -0700 @@ -27,17 +27,6 @@ ## For the definition of this method see ## @url{http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method}. ## -## First output argument is the final integration time value. -## -## Second output parameter is the higher order computed solution at time -## @var{t_next} (local extrapolation). -## -## Third output parameter is a lower order solution for the estimation of the -## error. -## -## Fourth output parameter is matrix containing the Runge-Kutta evaluations -## to use in an FSAL scheme or for dense output. -## ## First input argument is the function describing the system of ODEs to be ## integrated. ## @@ -53,9 +42,18 @@ ## ## Sixth input parameter is optional and describes the Runge-Kutta evaluations ## of the previous step to use in an FSAL scheme. -## @end deftypefn +## +## First output argument is the final integration time value. +## +## Second output parameter is the higher order computed solution at time +## @var{t_next} (local extrapolation). ## -## @seealso{odepkg} +## Third output parameter is a lower order solution for the estimation of the +## error. +## +## Fourth output parameter is matrix containing the Runge-Kutta evaluations +## to use in an FSAL scheme or for dense output. +## @end deftypefn function [t_next, x_next, x_est, k] = runge_kutta_45_dorpri (f, t, x, dt, options = [], @@ -68,10 +66,12 @@ 44/45 -56/15 32/9 0 0 0; 19372/6561 -25360/2187 64448/6561 -212/729 0 0; 9017/3168 -355/33 46732/5247 49/176 -5103/18656 0]; - persistent b = [0 1/5 3/10 4/5 8/9 1 1]; - persistent c = [(35/384) 0 (500/1113) (125/192) (-2187/6784) (11/84)]; - persistent c_prime = [(5179/57600) 0 (7571/16695) (393/640), ... - (-92097/339200) (187/2100) (1/40)]; + persistent b = [0, 1/5, 3/10, 4/5, 8/9, 1, 1]; + persistent c = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]; + persistent c_prime = [5179/57600, 0, 7571/16695, 393/640, ... + -92097/339200, 187/2100, 1/40]; + ## FIXME: Which source is c_prime derived from? + ## Can't the Shampine clause be deleted if it will never be used? ## According to Shampine 1986: ## persistent c_prime = [(1951/21600) 0 (22642/50085) (451/720), ... ## (-12231/42400) (649/6300) (1/60)];
--- a/scripts/ode/private/runge_kutta_interpolate.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/runge_kutta_interpolate.m Sat Oct 15 22:39:29 2016 -0700 @@ -21,10 +21,8 @@ switch (order) - ## FIXME: Can interpolations for orders 1-4 simply be deleted? 2015-10-14. - case 1 - u_interp = interp1 (z, u', t, 'linear'); + u_interp = interp1 (z, u', t, "linear"); case 2 if (! isempty (k_vals)) @@ -33,9 +31,11 @@ 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); + u_interp = hermite_cubic_interpolation (z, u, k_vals, t); + + ## FIXME: Do we need an algorithm for order = 4? #{ case 4 ## if ode45 is used without local extrapolation this function @@ -48,17 +48,6 @@ ## 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"]); @@ -71,7 +60,7 @@ endfunction ## The function below can be used in an ODE solver to interpolate the solution -## at the time t_out using 2th order hermite interpolation. +## at the time t_out using 2nd order Hermite interpolation. function x_out = quadratic_interpolation (t, x, der, t_out) # coefficients of the parabola @@ -79,19 +68,32 @@ b = der(:) - 2*t(1).*a; c = x(:,1) - a*t(1)^2 - b*t(1); - # evauate in t_out + # evaluate in t_out x_out = a*t_out.^2 + b*t_out + c; 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. +## solution at the time t_out using 3rd order Hermite interpolation. +function x_out = hermite_cubic_interpolation (t, x, der, t_out) + + dt = (t(2) - t(1)); + s = (t_out - t(1)) / dt; + x_out = ((1 + 2*s) .* (1-s).^2) .* x(:,1) + ... + (s .* (1-s).^2 * dt ) .* der(:,1) + ... + ((3-2*s) .* s.^2 ) .* x(:,2) + ... + ((s-1) .* s.^2 * dt ) .* der(:,2); + +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)].'; + [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. @@ -116,16 +118,3 @@ endfunction -## The function below can be used in an ODE solver to interpolate the -## solution at the time t_out using 3rd order hermite interpolation. -function x_out = hermite_cubic_interpolation (t, x, der, t_out) - - dt = (t(2) - t(1)); - s = (t_out - t(1)) / dt; - x_out = ((1 + 2*s) .* (1-s).^2) .* x(:,1) + ... - (s .* (1-s).^2 * dt ) .* der(:,1) + ... - ((3-2*s) .* s.^2 ) .* x(:,2) + ... - ((s-1) .* s.^2 * dt ) .* der(:,2); - -endfunction -
--- a/scripts/ode/private/starting_stepsize.m Sat Oct 15 10:30:48 2016 +0200 +++ b/scripts/ode/private/starting_stepsize.m Sat Oct 15 22:39:29 2016 -0700 @@ -40,14 +40,14 @@ args = {}) ## compute norm of initial conditions - d0 = AbsRel_Norm (x0, x0, AbsTol, RelTol, normcontrol); + d0 = AbsRel_norm (x0, x0, AbsTol, RelTol, normcontrol); ## compute norm of the function evaluated at initial conditions y = func (t0, x0, args{:}); if (iscell (y)) y = y{1}; endif - d1 = AbsRel_Norm (y, y, AbsTol, RelTol, normcontrol); + d1 = AbsRel_norm (y, y, AbsTol, RelTol, normcontrol); if (d0 < 1e-5 || d1 < 1e-5) h0 = 1e-6; @@ -64,7 +64,7 @@ yh = yh{1}; endif d2 = (1 / h0) * ... - AbsRel_Norm (yh - y, yh - y, AbsTol, RelTol, normcontrol); + AbsRel_norm (yh - y, yh - y, AbsTol, RelTol, normcontrol); if (max (d1, d2) <= 1e-15) h1 = max (1e-6, h0 * 1e-3);