Mercurial > octave
changeset 22631:ddbc3f9f0136
maint: merge stable to default.
| author | Rik <rik@octave.org> |
|---|---|
| date | Sun, 16 Oct 2016 12:25:24 -0700 |
| parents | 5ab3c91fc4bb (current diff) a918e983a943 (diff) |
| children | 859d48b5648d |
| files | NEWS libinterp/corefcn/graphics.cc libinterp/corefcn/graphics.in.h scripts/ode/private/AbsRel_Norm.m scripts/ode/private/known_option_names.m scripts/plot/appearance/legend.m |
| diffstat | 29 files changed, 526 insertions(+), 540 deletions(-) [+] |
line wrap: on
line diff
--- a/NEWS Sun Oct 16 02:23:54 2016 -0500 +++ b/NEWS Sun Oct 16 12:25:24 2016 -0700 @@ -70,6 +70,10 @@ "luminance profile" and is also more similar to Matlab's new default colormap "parula". + ** The colormap function no longer supports the input argument "list" + to show built-in colormaps. Use "help colormap" to find the + built-in colormaps. + ** The graphics command "hold on" now ensures that each new plot added to an existing plot has a different color or linestyle according to the "ColorOrder" and/or "LineStyleOrder" properties. This is
--- a/libinterp/corefcn/graphics.cc Sun Oct 16 02:23:54 2016 -0500 +++ b/libinterp/corefcn/graphics.cc Sun Oct 16 12:25:24 2016 -0700 @@ -8809,6 +8809,14 @@ error ("set: cannot change the style of a uicontrol object after creation."); style = st; + + // if we know know what we are, can override value for listbox and popupmenu + if (style_is ("listbox") || style_is ("popupmenu")) + { + Matrix v = value.get ().matrix_value (); + if(v.numel () == 1 && v (0) == 0) + value.set (octave_value (1)); + } } Matrix
--- a/libinterp/corefcn/graphics.in.h Sun Oct 16 02:23:54 2016 -0500 +++ b/libinterp/corefcn/graphics.in.h Sun Oct 16 12:25:24 2016 -0700 @@ -5517,7 +5517,7 @@ radio_property style S , "{pushbutton}|togglebutton|radiobutton|checkbox|edit|text|slider|frame|listbox|popupmenu" string_property tooltipstring , "" radio_property units u , "normalized|inches|centimeters|points|{pixels}|characters" - row_vector_property value , Matrix (1, 1, 1.0) + row_vector_property value , Matrix (1, 1, 0.0) radio_property verticalalignment , "top|{middle}|bottom" END_PROPERTIES
--- a/scripts/general/randi.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/general/randi.m Sun Oct 16 12:25:24 2016 -0700 @@ -163,13 +163,13 @@ ## Test that no warning thrown if IMAX is exactly on the limits of the range %!function test_no_warning (func, varargin) -%! state = warning ("query"); -%! unwind_protect -%! warning ("error", "all"); -%! func (varargin{:}); -%! unwind_protect_cleanup -%! warning (state); -%! end_unwind_protect +%! state = warning ("query"); +%! unwind_protect +%! warning ("error", "all"); +%! func (varargin{:}); +%! unwind_protect_cleanup +%! warning (state); +%! end_unwind_protect %!endfunction %!test test_no_warning (@randi, max_int8, "int8"); %!test test_no_warning (@randi, max_uint8, "uint8");
--- a/scripts/general/structfun.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/general/structfun.m Sun Oct 16 12:25:24 2016 -0700 @@ -121,8 +121,8 @@ %! assert (o, l); %!function [a, b] = __twoouts (x) -%! a = x + x; -%! b = x * x; +%! a = x + x; +%! b = x * x; %!endfunction %!test
--- a/scripts/image/colormap.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/image/colormap.m Sun Oct 16 12:25:24 2016 -0700 @@ -45,7 +45,39 @@ ## For convenience, it is also possible to use this function with the ## command form, @code{colormap @var{map_name}}. ## -## @seealso{viridis} +## The list of built-in colormaps is: +## +## @c FIXME: It would be nice to display the actual colormap as an image +## @c in the PDF version of the documentation. +## @multitable @columnfractions 0.15 .85 +## @headitem Map @tab Description +## @item viridis @tab default +## @item jet @tab colormap traversing blue, cyan, green, yellow, red. +## @item cubehelix @tab colormap traversing black, blue, green, red, white with increasing intensity. +## @item hsv @tab cyclic colormap traversing Hue, Saturation, Value space. +## @item rainbow @tab colormap traversing red, yellow, blue, green, violet. +## @item ------------- @tab --------------------------------------------------------------------------------------------- +## @item hot @tab colormap traversing black, red, orange, yellow, white. +## @item cool @tab colormap traversing cyan, purple, magenta. +## @item spring @tab colormap traversing magenta to yellow. +## @item summer @tab colormap traversing green to yellow. +## @item autumn @tab colormap traversing red, orange, yellow. +## @item winter @tab colormap traversing blue to green. +## @item ------------- @tab --------------------------------------------------------------------------------------------- +## @item gray @tab colormap traversing black to white in shades of gray. +## @item bone @tab colormap traversing black, gray-blue, white. +## @item copper @tab colormap traversing black to light copper. +## @item pink @tab colormap traversing black, gray-pink, white. +## @item ocean @tab colormap traversing black, dark-blue, white. +## @item ------------- @tab --------------------------------------------------------------------------------------------- +## @item colorcube @tab equally spaced colors in RGB color space. +## @item flag @tab cyclic 4-color map of red, white, blue, black. +## @item lines @tab cyclic colormap with colors from axes @qcode{"ColorOrder"} property. +## @item prism @tab cyclic 6-color map of red, orange, yellow, green, blue, violet. +## @item ------------- @tab --------------------------------------------------------------------------------------------- +## @item white @tab all white colormap (no colors). +## @end multitable +## @seealso{viridis, jet, cubehelix, hsv, rainbow, hot, cool, spring, summer, autumn, winter, gray, bone, copper, pink, ocean, colorcube, flag, lines, prism, white} ## @end deftypefn ## Author: Tony Richardson <arichard@stark.cc.oh.us>
--- a/scripts/linear-algebra/normest1.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/linear-algebra/normest1.m Sun Oct 16 12:25:24 2016 -0700 @@ -248,28 +248,28 @@ endfunction %!function z = afun_A (flag, x, A, n) -%! switch flag -%! case {"dim"} -%! z = n; -%! case {"real"} -%! z = isreal (A); -%! case {"transp"} -%! z = A' * x; -%! case {"notransp"} -%! z = A * x; -%! endswitch +%! switch flag +%! case {"dim"} +%! z = n; +%! case {"real"} +%! z = isreal (A); +%! case {"transp"} +%! z = A' * x; +%! case {"notransp"} +%! z = A * x; +%! endswitch %!endfunction %!function z = afun_A_P (flag, x, A, m) -%! switch flag -%! case "dim" -%! z = length (A); -%! case "real" -%! z = isreal (A); -%! case "transp" -%! z = x; for i = 1:m, z = A' * z;, endfor -%! case "notransp" -%! z = x; for i = 1:m, z = A * z;, endfor -%! endswitch +%! switch flag +%! case "dim" +%! z = length (A); +%! case "real" +%! z = isreal (A); +%! case "transp" +%! z = x; for i = 1:m, z = A' * z;, endfor +%! case "notransp" +%! z = x; for i = 1:m, z = A * z;, endfor +%! endswitch %!endfunction %!test
--- a/scripts/ode/module.mk Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/module.mk Sun Oct 16 12:25:24 2016 -0700 @@ -3,10 +3,9 @@ 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/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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/ode23.m Sun Oct 16 12:25:24 2016 -0700 @@ -1,4 +1,5 @@ -## 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> @@ -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 @@ -85,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 @@ -110,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 = {}; @@ -128,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)) @@ -157,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", ... @@ -180,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; @@ -191,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 @@ -207,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 @@ -224,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 @@ -306,7 +303,6 @@ %!demo -%! %! ## Demonstrate convergence order for ode23 %! tol = 1e-5 ./ 10.^[0:8]; %! for i = 1 : numel (tol) @@ -335,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 @@ -486,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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/ode45.m Sun Oct 16 12:25:24 2016 -0700 @@ -1,4 +1,5 @@ -## 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> @@ -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 @@ -76,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) @@ -93,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 = {}; @@ -111,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)) @@ -141,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"}); @@ -162,9 +160,9 @@ "InitialSlope", "MaxOrder", "BDF"}); attributes = rmfield (attributes, {"Jacobian", "JPattern", "Vectorized", ... "MvPattern", "MassSingular", ... - "InitialSlope", "MaxOrder", "BDF"}); + "InitialSlope", "MaxOrder", "BDF"}); - odeopts = odemergeopts (odeopts, defaults, classes, attributes, 'ode45'); + odeopts = odemergeopts ("ode45", odeopts, defaults, classes, attributes); odeopts.funarguments = funarguments; odeopts.direction = direction; @@ -175,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 @@ -191,13 +189,11 @@ 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); - endif - + 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; @@ -208,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 @@ -290,7 +284,6 @@ %!demo -%! %! ## Demonstrate convergence order for ode45 %! tol = 1e-5 ./ 10.^[0:8]; %! for i = 1 : numel (tol) @@ -307,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)"); @@ -315,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 @@ -488,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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/odeget.m Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/odeplot.m Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/odeset.m Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,50 +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) - - n = length (x); - - if (nargin == 5) - y = zeros (size (x)); - endif - - 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); - 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 Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/integrate_adaptive.m Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/kahan.m Sun Oct 16 12:25:24 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/known_option_names.m Sun Oct 16 02:23:54 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/ode_event_handler.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/ode_event_handler.m Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/odedefaults.m Sun Oct 16 12:25:24 2016 -0700 @@ -1,4 +1,5 @@ -## 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. ## @@ -16,81 +17,86 @@ ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. +## -*- texinfo -*- +## @deftypefn {} {[@var{defaults}, @var{classes}, @var{attributes}] =} 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"); + 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"); -##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 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/odemergeopts.m Sun Oct 16 12:25:24 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); - + + options.(key) = useroptions.(key); + endif + endfor + endfunction +
--- a/scripts/ode/private/runge_kutta_23.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/runge_kutta_23.m Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/runge_kutta_45_dorpri.m Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/runge_kutta_interpolate.m Sun Oct 16 12:25:24 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 Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/ode/private/starting_stepsize.m Sun Oct 16 12:25:24 2016 -0700 @@ -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. @@ -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,15 +64,15 @@ 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); + 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
--- a/scripts/plot/appearance/legend.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/plot/appearance/legend.m Sun Oct 16 12:25:24 2016 -0700 @@ -1428,7 +1428,8 @@ %! xlabel ("Indices"); %! ylabel ("Random Values"); %! title ('Legend "off" deletes the legend'); -%! legend (cellstr (num2str ((1:10)")), "location", "northeastoutside"); +%! legend (cellstr (num2str ((0:10)')), "location", "northeastoutside"); +%! pause (1); %! legend off; %! axis ([0, 10, 0 1]);
--- a/scripts/plot/util/frame2im.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/plot/util/frame2im.m Sun Oct 16 12:25:24 2016 -0700 @@ -59,11 +59,11 @@ %!function f = make_rgb_f () -%! f = randi ([0 255], 10, 20, 3); +%! f = randi ([0 255], 10, 20, 3); %!endfunction %!function f = make_ind_f () -%! f = randi ([1 100], 10, 20, 3); +%! f = randi ([1 100], 10, 20, 3); %!endfunction %!test
--- a/scripts/plot/util/im2frame.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/plot/util/im2frame.m Sun Oct 16 12:25:24 2016 -0700 @@ -69,11 +69,11 @@ %!function f = make_rgb_f () -%! f = randi ([0 255], 10, 20, 3); +%! f = randi ([0 255], 10, 20, 3); %!endfunction %!function f = make_ind_f () -%! f = randi ([1 100], 10, 20, 3); +%! f = randi ([1 100], 10, 20, 3); %!endfunction %!test
--- a/scripts/testfun/test.m Sun Oct 16 02:23:54 2016 -0500 +++ b/scripts/testfun/test.m Sun Oct 16 12:25:24 2016 -0700 @@ -885,19 +885,19 @@ ## Test 'function' keyword %!function x = __test_a (y) -%! x = 2*y; +%! x = 2*y; %!endfunction %!assert (__test_a (2), 4) # Test a test function %!function __test_a (y) -%! x = 2*y; +%! x = 2*y; %!endfunction %!test %! __test_a (2); # Test a test function with no return value %!function [x,z] = __test_a (y) -%! x = 2*y; -%! z = 3*y; +%! x = 2*y; +%! z = 3*y; %!endfunction %!test %! [x,z] = __test_a (3); # Test a test function with multiple returns