Mercurial > octave-libtiff
changeset 22625:081a201b77c7 stable
Clean up ode options implementation to follow Octave coding standards.
* known_option_names.m: Delete file
* scripts/ode/module.mk: Remove known_option_names from build system.
* ode23.m, ode45.m: Fix typo in docstring. Correct indentation.
Remove trailing whitespace.
* AbsRel_Norm.m: Use default for input argument to simplify function.
Remove input validation for private, internal function.
* odedefaults.m: Add docstring. Use persistent variables for performance.
* odemergeopts.m: Fix indentation.
* starting_stepsize.m: Show input func as '@func' in docstring.
author | Carlo de Falco <carlo.defalco@polimi.it> |
---|---|
date | Sat, 15 Oct 2016 10:30:48 +0200 |
parents | 37b7b86f62f2 |
children | 869c02fde46c |
files | scripts/ode/module.mk scripts/ode/ode23.m scripts/ode/ode45.m scripts/ode/private/AbsRel_Norm.m scripts/ode/private/known_option_names.m scripts/ode/private/odedefaults.m scripts/ode/private/odemergeopts.m scripts/ode/private/starting_stepsize.m |
diffstat | 8 files changed, 108 insertions(+), 150 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ode/module.mk Sat Oct 15 02:51:30 2016 -0500 +++ b/scripts/ode/module.mk Sat Oct 15 10:30:48 2016 +0200 @@ -6,7 +6,6 @@ scripts/ode/private/AbsRel_Norm.m \ scripts/ode/private/integrate_adaptive.m \ scripts/ode/private/kahan.m \ - scripts/ode/private/known_option_names.m \ scripts/ode/private/odedefaults.m \ scripts/ode/private/odemergeopts.m \ scripts/ode/private/ode_event_handler.m \
--- a/scripts/ode/ode23.m Sat Oct 15 02:51:30 2016 -0500 +++ b/scripts/ode/ode23.m Sat Oct 15 10:30:48 2016 +0200 @@ -1,3 +1,4 @@ +## 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> @@ -43,8 +44,8 @@ ## ## By default, @code{ode23} uses an adaptive timestep with the ## @code{integrate_adaptive} algorithm. The tolerance for the timestep -## computation may be changed by using the options @qcode{"RelTol"}, -## and @qcode{"AbsTol"},. +## computation may be changed by using the options @qcode{"RelTol"} +## and @qcode{"AbsTol"}. ## ## @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 @@ -211,7 +212,7 @@ init, odeopts.AbsTol, odeopts.RelTol, strcmp (odeopts.NormControl, - "on"), odeopts.funarguments); + "on"), odeopts.funarguments); endif @@ -231,18 +232,18 @@ 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{:}); + \ fun (t, x, odeopts.funarguments{:}); else # if ((! strcmp (odeopts.MStateDependence, "none")) == false) 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, ... order, fun, trange, init, odeopts); - + ## Postprocessing, do whatever when terminating integration algorithm if (odeopts.haveoutputfunction) # Cleanup plotter
--- a/scripts/ode/ode45.m Sat Oct 15 02:51:30 2016 -0500 +++ b/scripts/ode/ode45.m Sat Oct 15 10:30:48 2016 +0200 @@ -1,3 +1,4 @@ +## 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> @@ -41,7 +42,7 @@ ## ## 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 options @qcode{"RelTol"}, +## computation may be changed by using the options @qcode{"RelTol"} ## and @qcode{"AbsTol"}. ## ## @var{init} contains the initial value for the unknowns. If it is a row @@ -162,7 +163,7 @@ "InitialSlope", "MaxOrder", "BDF"}); attributes = rmfield (attributes, {"Jacobian", "JPattern", "Vectorized", ... "MvPattern", "MassSingular", ... - "InitialSlope", "MaxOrder", "BDF"}); + "InitialSlope", "MaxOrder", "BDF"}); odeopts = odemergeopts (odeopts, defaults, classes, attributes, 'ode45'); @@ -195,8 +196,8 @@ init, odeopts.AbsTol, odeopts.RelTol, strcmp (odeopts.NormControl, - "on"), odeopts.funarguments); - endif + "on"), odeopts.funarguments); + endif if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass)) @@ -215,18 +216,18 @@ 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{:}); + \ fun (t, x, odeopts.funarguments{:}); else # if ((! strcmp (odeopts.MStateDependence, "none")) == false) 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
--- a/scripts/ode/private/AbsRel_Norm.m Sat Oct 15 02:51:30 2016 -0500 +++ b/scripts/ode/private/AbsRel_Norm.m Sat Oct 15 10:30:48 2016 +0200 @@ -22,23 +22,10 @@ ## Undocumented internal function. ## @end deftypefn -function retval = AbsRel_Norm (x, x_old, AbsTol, RelTol, normcontrol, y) - - n = length (x); - - if (nargin == 5) - y = zeros (size (x)); - endif +function retval = AbsRel_Norm (x, x_old, AbsTol, RelTol, normcontrol, y = zeros (size (x))) - 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 - + n = numel (x); + sc = AbsTol + max (abs (x), abs (x_old)) .* RelTol; if (normcontrol) retval = max (abs (x - y) ./ sc);
--- a/scripts/ode/private/known_option_names.m Sat Oct 15 02:51:30 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,35 +0,0 @@ -## Copyright (C) 2016 Carlo de Falco -## -## 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 {} {@var{val} =} known_option_names () -## Return a list of known names for ode options. -## @seealso{odeset, odeget} -## @end deftypefn - -function ret = known_option_names () - -ret = {"AbsTol"; "BDF"; "Events"; "InitialSlope"; - "InitialStep"; "Jacobian"; "JConstant"; "JPattern"; - "Mass"; "MassConstant"; "MassSingular"; "MaxOrder"; - "MaxStep"; "MStateDependence"; "MvPattern"; - "NonNegative"; "NormControl"; "OutputFcn"; "OutputSel"; - "Refine"; "RelTol"; "Stats"; "Vectorized"; - "TimeStepSize"; "TimeStepNumber"}; - -endfunction
--- a/scripts/ode/private/odedefaults.m Sat Oct 15 02:51:30 2016 -0500 +++ b/scripts/ode/private/odedefaults.m Sat Oct 15 10:30:48 2016 +0200 @@ -1,3 +1,4 @@ +## Copyright (C) 2016, Carlo de Falco ## Copyright (C) 2016, Francesco Faccio <francesco.faccio@mail.polimi.it> ## ## This file is part of Octave. @@ -16,81 +17,85 @@ ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. +## -*- texinfo -*- +## @deftypefn {} {[@var{d}, @var{c}, @var{a}] =} odedefaults (@var{n}, @var{t0}, @var{tf}) +## Undocumented internal function. +## @end deftypefn + function [defaults, classes, attributes] = odedefaults (n, t0, tf) -defaults = odeset ('AbsTol', 1e-6, - 'BDF', 'off', - 'Events', [], - 'InitialSlope', zeros(n,1), - 'InitialStep', [], - 'Jacobian', [], - 'JConstant', 'off', - 'JPattern', [], - 'Mass', [], - 'MassConstant', 'off', - 'MassSingular', 'maybe', - 'MaxOrder', 5, - 'MaxStep', 0.1*abs(t0-tf), - 'MStateDependence', 'weak', - 'MvPattern', [], - 'NonNegative', [], - 'NormControl', 'off', - 'OutputFcn', [], - 'OutputSel', [], - 'Refine', 1, - 'RelTol', 1e-3, - 'Stats', 'off', - 'Vectorized', 'off'); - -classes = odeset ('Abstol', {"float"}, - 'BDF', "char", - 'Events', {"function_handle"}, - 'InitialSlope', {"float"}, - 'InitialStep', {"float"}, - 'Jacobian', {"float", "function_handle", "cell"}, - 'JConstant', "char", - 'JPattern', {"float"}, - 'Mass', {"float", "function_handle"}, - 'MassConstant', "char", - 'MassSingular', "char", - 'MaxOrder', {"float"}, - 'MaxStep', {"float"}, - 'MStateDependence', "char", - 'MvPattern', {"float"}, - 'NonNegative', {"float"}, - 'NormControl', "char", - 'OutputFcn', {"function_handle"}, - 'OutputSel', {"float"}, - 'Refine', {"float"}, - 'RelTol', {"float"}, - 'Stats', "char", - 'Vectorized', "char"); - - -##FIXME: How can I check Jacobian where it's a cell????? Maybe it's better to check it inside the solver -##FIXME: Vectorized can be a cell of stings -attributes = odeset ('AbsTol', {"real", "vector", "positive"}, - 'BDF', {"on", "off"}, - 'Events', {}, - 'InitialSlope', {"real", "vector", "numel", n}, - 'InitialStep', {"positive", "scalar"}, - 'Jacobian', {}, - 'JConstant', {"on", "off"}, - 'JPattern', {"vector"}, - 'Mass', {}, - 'MassConstant', {"on", "off"}, - 'MassSingular', {"no", "maybe", "yes"}, - 'MaxOrder', {">=", 0, "<=", 5, "integer"}, - 'MaxStep', {"positive", "scalar", "real"}, - 'MStateDependence', {"weak", "strong", "none"}, - 'MvPattern', {"vector"}, - 'NonNegative', {"vector", "integer", "positive"}, - 'NormControl', {"on", "off"}, - 'OutputFcn', {}, - 'OutputSel', {"vector", "integer", "positive",... - ">", 0, "<=", n}, - 'Refine', {"scalar", ">", 0, "integer"}, - 'RelTol', {"scalar", "positive", "real"}, - 'Stats', {"on", "off"}, - 'Vectorized', {"on", "off"}); + persistent defaults = struct ("AbsTol", 1e-6, + "BDF", "off", + "Events", [], + "InitialSlope", zeros(n,1), + "InitialStep", [], + "Jacobian", [], + "JConstant", "off", + "JPattern", [], + "Mass", [], + "MassConstant", "off", + "MassSingular", "maybe", + "MaxOrder", 5, + "MaxStep", 0.1 * abs (t0-tf), + "MStateDependence", "weak", + "MvPattern", [], + "NonNegative", [], + "NormControl", "off", + "OutputFcn", [], + "OutputSel", [], + "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"}}, + "InitialSlope", {{"float"}}, + "InitialStep", {{"float"}}, + "Jacobian", {{"float", "function_handle", "cell"}}, + "JConstant", "char", + "JPattern", {{"float"}}, + "Mass", {{"float", "function_handle"}}, + "MassConstant", "char", + "MassSingular", "char", + "MaxOrder", {{"float"}}, + "MaxStep", {{"float"}}, + "MStateDependence", "char", + "MvPattern", {{"float"}}, + "NonNegative", {{"float"}}, + "NormControl", "char", + "OutputFcn", {{"function_handle"}}, + "OutputSel", {{"float"}}, + "Refine", {{"float"}}, + "RelTol", {{"float"}}, + "Stats", "char", + "Vectorized", "char"); + + persistent attributes = struct ("AbsTol", {{"real", "vector", "positive"}}, + "BDF", {{"on", "off"}}, + "Events", {{}}, + "InitialSlope", {{"real", "vector", "numel", n}}, + "InitialStep", {{"positive", "scalar"}}, + "Jacobian", {{}}, + "JConstant", {{"on", "off"}}, + "JPattern", {{"vector"}}, + "Mass", {{}}, + "MassConstant", {{"on", "off"}}, + "MassSingular", {{"no", "maybe", "yes"}}, + "MaxOrder", {{">=", 0, "<=", 5, "integer"}}, + "MaxStep", {{"positive", "scalar", "real"}}, + "MStateDependence", {{"weak", "strong", "none"}}, + "MvPattern", {{"vector"}}, + "NonNegative", {{"vector", "integer", "positive"}}, + "NormControl", {{"on", "off"}}, + "OutputFcn", {{}}, + "OutputSel", {{"vector", "integer", "positive",... + ">", 0, "<=", n}}, + "Refine", {{"scalar", ">", 0, "integer"}}, + "RelTol", {{"scalar", "positive", "real"}}, + "Stats", {{"on", "off"}}, + "Vectorized", {{"on", "off"}}); endfunction
--- a/scripts/ode/private/odemergeopts.m Sat Oct 15 02:51:30 2016 -0500 +++ b/scripts/ode/private/odemergeopts.m Sat Oct 15 10:30:48 2016 +0200 @@ -32,11 +32,11 @@ else error ("Octave:invalid-input-arg", - [fun_name ": invalid value assigned to field '%s'"], key); + [fun_name ": invalid value assigned to field '%s'"], key); endif - - options.(key) = useroptions.(key); - + + options.(key) = useroptions.(key); + endif endfor endfunction
--- a/scripts/ode/private/starting_stepsize.m Sat Oct 15 02:51:30 2016 -0500 +++ b/scripts/ode/private/starting_stepsize.m Sat Oct 15 10:30:48 2016 +0200 @@ -17,12 +17,12 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {} {@var{h} =} starting_stepsize (@var{order}, @var{@@func}, @var{t0}, @var{x0}, @var{AbsTol}, @var{RelTol}, @var{normcontrol}) +## @deftypefn {} {@var{h} =} starting_stepsize (@var{order}, @var{func}, @var{t0}, @var{x0}, @var{AbsTol}, @var{RelTol}, @var{normcontrol}) ## ## Determine a good initial timestep for an ODE solver of order @var{order} ## using the algorithm described in reference [1]. ## -## The input argument @var{@@func}, is the function describing the differential +## The input argument @var{func}, is the function describing the differential ## 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. @@ -67,12 +67,12 @@ AbsRel_Norm (yh - y, yh - y, AbsTol, RelTol, normcontrol); if (max (d1, d2) <= 1e-15) - h1 = max (1e-6, h0*1e-3); + h1 = max (1e-6, h0 * 1e-3); else h1 = (1e-2 / max (d1, d2)) ^(1 / (order+1)); endif - h = min (100*h0, h1); + h = min (100 * h0, h1); endfunction