Mercurial > octave
changeset 22663:9a939479308f
maint: Periodic merge of stable to default.
| author | Rik <rik@octave.org> |
|---|---|
| date | Mon, 24 Oct 2016 09:04:30 -0700 |
| parents | a93887d7f0da (current diff) 655157b34a9f (diff) |
| children | f1bb2f0bcfec |
| files | |
| diffstat | 9 files changed, 173 insertions(+), 130 deletions(-) [+] |
line wrap: on
line diff
--- a/libinterp/corefcn/fft.cc Sun Oct 23 13:57:00 2016 -0700 +++ b/libinterp/corefcn/fft.cc Mon Oct 24 09:04:30 2016 -0700 @@ -49,8 +49,9 @@ octave_value retval; octave_value arg = args(0); + octave_idx_type n_points = -1; dim_vector dims = arg.dims (); - octave_idx_type n_points = -1; + int ndims = dims.ndims (); int dim = -1; if (nargin > 1) @@ -72,7 +73,7 @@ double dval = args(2).double_value (); if (octave::math::isnan (dval)) error ("%s: DIM cannot be NaN", fcn); - else if (dval < 1 || dval > dims.ndims ()) + else if (dval < 1 || dval > ndims) error ("%s: DIM must be a valid dimension along which to perform FFT", fcn); else @@ -80,21 +81,19 @@ dim = octave::math::nint (dval) - 1; } - for (octave_idx_type i = 0; i < dims.ndims (); i++) + // FIXME: This seems strange and unnecessary (10/21/16). + // How would you ever arrive at an octave_value object without correct dims? + // We certainly don't make this check every other place in Octave. + for (octave_idx_type i = 0; i < ndims; i++) if (dims(i) < 0) return retval; if (dim < 0) { - for (octave_idx_type i = 0; i < dims.ndims (); i++) - if (dims(i) > 1) - { - dim = i; - break; - } + dim = dims.first_non_singleton (); // And if the first argument is scalar? - if (dim < 0) + if (dim == ndims) dim = 1; } @@ -103,7 +102,7 @@ else dims(dim) = n_points; - if (dims.any_zero () || n_points == 0) + if (n_points == 0 || dims.any_zero ()) { if (arg.is_single_type ()) return octave_value (FloatNDArray (dims)); @@ -111,6 +110,16 @@ return octave_value (NDArray (dims)); } + if (n_points == 1) + { + octave_value_list idx (ndims); + for (octave_idx_type i = 0; i < ndims; i++) + idx(i) = idx_vector::colon; + idx(dim) = idx_vector (0); + + return arg.do_index_op (idx); + } + if (arg.is_single_type ()) { if (arg.is_real_type ()) @@ -230,9 +239,9 @@ } /* -%% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) -%% Comalco Research and Technology -%% 02 May 2000 +## Author: David Billinghurst (David.Billinghurst@riotinto.com.au) +## Comalco Research and Technology +## 02 May 2000 %!test %! N = 64; %! n = 4; @@ -246,9 +255,9 @@ %! %! assert (S, answer, 4*N*eps); -%% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) -%% Comalco Research and Technology -%% 02 May 2000 +## Author: David Billinghurst (David.Billinghurst@riotinto.com.au) +## Comalco Research and Technology +## 02 May 2000 %!test %! N = 64; %! n = 7; @@ -261,9 +270,9 @@ %! %! assert (ifft (S), s, 4*N*eps); -%% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) -%% Comalco Research and Technology -%% 02 May 2000 +## Author: David Billinghurst (David.Billinghurst@riotinto.com.au) +## Comalco Research and Technology +## 02 May 2000 %!test %! N = 64; %! n = 4; @@ -277,9 +286,9 @@ %! %! assert (S, answer, 4*N*eps ("single")); -%% Author: David Billinghurst (David.Billinghurst@riotinto.com.au) -%% Comalco Research and Technology -%% 02 May 2000 +## Author: David Billinghurst (David.Billinghurst@riotinto.com.au) +## Comalco Research and Technology +## 02 May 2000 %!test %! N = 64; %! n = 7;
--- a/scripts/ode/ode23.m Sun Oct 23 13:57:00 2016 -0700 +++ b/scripts/ode/ode23.m Mon Oct 24 09:04:30 2016 -0700 @@ -60,8 +60,9 @@ ## unknown of the problem and each row corresponds to a time in @var{t}. ## ## The output can also be returned as a structure @var{solution} which -## has field @var{x} containing the time where the solution was evaluated and -## field @var{y} containing the solution matrix for the times in @var{x}. +## has a field @var{x} containing a row vector of times where the solution +## was evaluated and a field @var{y} containing the solution matrix such +## that each column corresponds to a time in @var{x}. ## Use @code{fieldnames (@var{solution})} to see the other fields and ## additional information returned. ## @@ -155,15 +156,13 @@ [defaults, classes, attributes] = odedefaults (numel (init), trange(1), trange(end)); - defaults = rmfield (defaults, {"Jacobian", "JPattern", "Vectorized", ... - "MvPattern", "MassSingular", ... - "InitialSlope", "MaxOrder", "BDF"}); - classes = rmfield (classes, {"Jacobian", "JPattern", "Vectorized", ... - "MvPattern", "MassSingular", ... - "InitialSlope", "MaxOrder", "BDF"}); - attributes = rmfield (attributes, {"Jacobian", "JPattern", "Vectorized", ... - "MvPattern", "MassSingular", ... - "InitialSlope", "MaxOrder", "BDF"}); + persistent ode23_ignore_options = ... + {"BDF", "InitialSlope", "Jacobian", "JPattern", + "MassSingular", "MaxOrder", "MvPattern", "Vectorized"}; + + defaults = rmfield (defaults, ode23_ignore_options); + classes = rmfield (classes, ode23_ignore_options); + attributes = rmfield (attributes, ode23_ignore_options); odeopts = odemergeopts ("ode23", odeopts, defaults, classes, attributes); @@ -211,8 +210,7 @@ if (havemasshandle) # Handle only the dynamic mass matrix, if (! strcmp (odeopts.MStateDependence, "none")) - ## FIXME: How is this comment supposed to end? - ## constant mass matrices have already + ## constant mass matrices have already been handled mass = @(t,x) odeopts.Mass (t, x, odeopts.funarguments{:}); fun = @(t,x) mass (t, x, odeopts.funarguments{:}) ... \ fun (t, x, odeopts.funarguments{:}); @@ -223,6 +221,15 @@ endif endif + if (nargout == 1) + ## Single output requires auto-selected intermediate times, + ## which is obtained by NOT specifying specific solution times. + trange = [trange(1); trange(end)]; + odeopts.Refine = []; # disable Refine when single output requested + elseif (numel (trange) > 2) + odeopts.Refine = []; # disable Refine when specific times requested + endif + solution = integrate_adaptive (@runge_kutta_23, order, fun, trange, init, odeopts); @@ -232,7 +239,7 @@ endif if (! isempty (odeopts.Events)) # Cleanup event function handling ode_event_handler (odeopts.Events, solution.t(end), - solution.x(end,:)', "done", odeopts.funarguments{:}); + solution.x(end,:).', "done", odeopts.funarguments{:}); endif ## Print additional information if option Stats is set @@ -255,8 +262,8 @@ varargout{1} = solution.t; # Time stamps are first output argument varargout{2} = solution.x; # Results are second output argument elseif (nargout == 1) - varargout{1}.x = solution.t; # Time stamps are saved in field x - varargout{1}.y = solution.x; # Results are saved in field y + varargout{1}.x = solution.t.'; # Time stamps are saved in field x (row vector) + varargout{1}.y = solution.x.'; # Results are saved in field y (row vector) varargout{1}.solver = solver; # Solver name is saved in field solver if (! isempty (odeopts.Events)) varargout{1}.ie = solution.event{2}; # Index info which event occurred @@ -321,12 +328,6 @@ %!function ref = fref () # The computed reference sol %! 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]; -%!endfunction -%!function jac = fjcc (t, y, varargin) # sparse type -%! 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 @@ -391,41 +392,41 @@ %!test # Solve in backward direction starting at t=0 %! ref = [-1.205364552835178, 0.951542399860817]; %! sol = ode23 (@fpol, [0 -2], [2 0]); -%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 5e-3); +%! assert ([sol.x(end); sol.y(:,end)], [-2; ref'], 5e-3); %!test # Solve in backward direction starting at t=2 %! ref = [-1.205364552835178, 0.951542399860817]; %! sol = ode23 (@fpol, [2 0 -2], fref); -%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 2e-2); +%! assert ([sol.x(end); sol.y(:,end)], [-2; ref'], 2e-2); %!test # Solve another anonymous function in backward direction %! ref = [-1, 0.367879437558975]; %! sol = ode23 (@(t,y) y, [0 -1], 1); -%! assert ([sol.x(end), sol.y(end,:)], ref, 1e-2); +%! assert ([sol.x(end); sol.y(:,end)], ref', 1e-2); %!test # Solve another anonymous function below zero %! ref = [0, 14.77810590694212]; %! sol = ode23 (@(t,y) y, [-2 0], 2); -%! assert ([sol.x(end), sol.y(end,:)], ref, 1e-2); +%! assert ([sol.x(end); sol.y(:,end)], ref', 1e-2); %!test # Solve in backward direction starting at t=0 with MaxStep option %! ref = [-1.205364552835178, 0.951542399860817]; %! opt = odeset ("MaxStep", 1e-3); %! sol = ode23 (@fpol, [0 -2], [2 0], opt); %! assert ([abs(sol.x(8)-sol.x(7))], [1e-3], 1e-3); -%! assert ([sol.x(end), sol.y(end,:)], [-2, ref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [-2; ref'], 1e-3); %!test # AbsTol option %! opt = odeset ("AbsTol", 1e-5); %! sol = ode23 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # AbsTol and RelTol option %! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8); %! sol = ode23 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # RelTol and NormControl option -- higher accuracy %! opt = odeset ("RelTol", 1e-8, "NormControl", "on"); %! sol = ode23 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-4); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-4); %!test # Keeps initial values while integrating %! opt = odeset ("NonNegative", 2); %! sol = ode23 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, 2, 0], 1e-1); +%! assert ([sol.x(end); sol.y(:,end)], [2; 2; 0], 1e-1); %!test # Details of OutputSel and Refine can't be tested %! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); %! sol = ode23 (@fpol, [0 2], [2 0], opt); @@ -455,28 +456,41 @@ %!test # Mass option as function %! opt = odeset ("Mass", @fmas); %! sol = ode23 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # Mass option as matrix %! opt = odeset ("Mass", eye (2,2)); %! sol = ode23 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # Mass option as sparse matrix %! opt = odeset ("Mass", sparse (eye (2,2))); %! sol = ode23 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # Mass option as function and sparse matrix %! opt = odeset ("Mass", @fmsa); %! sol = ode23 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # Mass option as function and MStateDependence %! 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); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); -## FIXME: Missing tests. -## test for InitialSlope option is missing -## test for MaxOrder option is missing -## test for MvPattern option is missing +## Note: The following options have no effect on this solver +## therefore it makes no sense to test them here: +## +## "BDF" +## "InitialSlope" +## "JPattern" +## "Jacobian" +## "MassSingular" +## "MaxOrder" +## "MvPattern" +## "Vectorized" + +%!test # Check that imaginary part of solution does not get inverted +%! sol = ode23 (@(x,y) 1, [0 1], 1i); +%! assert (imag (sol.y), ones (size (sol.y))) +%! [x, y] = ode23 (@(x,y) 1, [0 1], 1i); +%! assert (imag (y), ones (size (y))) ## Test input validation %!error ode23 ()
--- a/scripts/ode/ode45.m Sun Oct 23 13:57:00 2016 -0700 +++ b/scripts/ode/ode45.m Mon Oct 24 09:04:30 2016 -0700 @@ -58,8 +58,9 @@ ## unknown of the problem and each row corresponds to a time in @var{t}. ## ## The output can also be returned as a structure @var{solution} which -## has field @var{x} containing the time where the solution was evaluated and -## field @var{y} containing the solution matrix for the times in @var{x}. +## has a field @var{x} containing a row vector of times where the solution +## was evaluated and a field @var{y} containing the solution matrix such +## that each column corresponds to a time in @var{x}. ## Use @code{fieldnames (@var{solution})} to see the other fields and ## additional information returned. ## @@ -146,16 +147,16 @@ [defaults, classes, attributes] = odedefaults (numel (init), trange(1), trange(end)); - defaults = odeset (defaults, "Refine", 4); - defaults = rmfield (defaults, {"Jacobian", "JPattern", "Vectorized", ... - "MvPattern", "MassSingular", ... - "InitialSlope", "MaxOrder", "BDF"}); - classes = rmfield (classes, {"Jacobian", "JPattern", "Vectorized", ... - "MvPattern", "MassSingular", ... - "InitialSlope", "MaxOrder", "BDF"}); - attributes = rmfield (attributes, {"Jacobian", "JPattern", "Vectorized", ... - "MvPattern", "MassSingular", ... - "InitialSlope", "MaxOrder", "BDF"}); + ## FIXME: Refine is not correctly implemented yet + defaults = odeset (defaults, "Refine", 4); + + persistent ode45_ignore_options = ... + {"BDF", "InitialSlope", "Jacobian", "JPattern", + "MassSingular", "MaxOrder", "MvPattern", "Vectorized"}; + + defaults = rmfield (defaults, ode45_ignore_options); + classes = rmfield (classes, ode45_ignore_options); + attributes = rmfield (attributes, ode45_ignore_options); odeopts = odemergeopts ("ode45", odeopts, defaults, classes, attributes); @@ -214,6 +215,15 @@ endif endif + if (nargout == 1) + ## Single output requires auto-selected intermediate times, + ## which is obtained by NOT specifying specific solution times. + trange = [trange(1); trange(end)]; + odeopts.Refine = []; # disable Refine when single output requested + elseif (numel (trange) > 2) + odeopts.Refine = []; # disable Refine when specific times requested + endif + solution = integrate_adaptive (@runge_kutta_45_dorpri, order, fun, trange, init, odeopts); @@ -223,7 +233,7 @@ endif if (! isempty (odeopts.Events)) # Cleanup event function handling ode_event_handler (odeopts.Events, solution.t(end), - solution.x(end,:)', "done", odeopts.funarguments{:}); + solution.x(end,:).', "done", odeopts.funarguments{:}); endif ## Print additional information if option Stats is set @@ -246,8 +256,8 @@ varargout{1} = solution.t; # Time stamps are first output argument varargout{2} = solution.x; # Results are second output argument elseif (nargout == 1) - varargout{1}.x = solution.t; # Time stamps are saved in field x - varargout{1}.y = solution.x; # Results are saved in field y + varargout{1}.x = solution.t.'; # Time stamps are saved in field x (row vector) + varargout{1}.y = solution.x.'; # Results are saved in field y (row vector) varargout{1}.solver = solver; # Solver name is saved in field solver if (! isempty (odeopts.Events)) varargout{1}.ie = solution.event{2}; # Index info which event occurred @@ -312,12 +322,6 @@ %!function ref = fref () # The computed reference solution %! 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]; -%!endfunction -%!function jac = fjcc (t, y, varargin) # sparse type -%! 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 @@ -385,54 +389,54 @@ %! opt = odeset ("MaxStep", 1e-3); %! sol = ode45 (@fpol, [0 0.2], [2 0], opt); %! assert ([sol.x(5)-sol.x(4)], [1e-3], 1e-3); -%!test # Solve with intermidiate step -%! sol = ode45 (@fpol, [0 1 2], [2 0]); -%! assert (any((sol.x-1) == 0)); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%!test # Solve with intermediate step +%! [t, y] = ode45 (@fpol, [0 1 2], [2 0]); +%! assert (any((t-1) == 0)); +%! assert ([t(end), y(end,:)], [2, fref], 1e-3); %!test # Solve in backward direction starting at t=0 %! vref = [-1.205364552835178, 0.951542399860817]; %! sol = ode45 (@fpol, [0 -2], [2 0]); -%! assert ([sol.x(end), sol.y(end,:)], [-2, vref], 1e-2); +%! assert ([sol.x(end); sol.y(:,end)], [-2; vref'], 1e-2); %!test # Solve in backward direction starting at t=2 %! vref = [-1.205364552835178, 0.951542399860817]; %! sol = ode45 (@fpol, [2 -2], fref); -%! assert ([sol.x(end), sol.y(end,:)], [-2, vref], 1e-2); -%!test # Solve in backward direction starting at t=2, with intermidiate step +%! assert ([sol.x(end); sol.y(:,end)], [-2; vref'], 1e-2); +%!test # Solve in backward direction starting at t=2, with intermediate step %! vref = [-1.205364552835178, 0.951542399860817]; -%! sol = ode45 (@fpol, [2 0 -2], fref); -%! idx = find(sol.x < 0, 1, "first") - 1; -%! assert ([sol.x(idx), sol.y(idx,:)], [0 2 0], 1e-2); -%! assert ([sol.x(end), sol.y(end,:)], [-2, vref], 1e-2); +%! [t, y] = ode45 (@fpol, [2 0 -2], fref); +%! idx = find(y < 0, 1, "first") - 1; +%! assert ([t(idx), y(idx,:)], [0,2,0], 1e-2); +%! assert ([t(end), y(end,:)], [-2, vref], 1e-2); %!test # Solve another anonymous function in backward direction %! vref = [-1, 0.367879437558975]; %! sol = ode45 (@(t,y) y, [0 -1], 1); -%! assert ([sol.x(end), sol.y(end,:)], vref, 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], vref', 1e-3); %!test # Solve another anonymous function below zero %! vref = [0, 14.77810590694212]; %! sol = ode45 (@(t,y) y, [-2 0], 2); -%! assert ([sol.x(end), sol.y(end,:)], vref, 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], vref', 1e-3); %!test # Solve in backward direction starting at t=0 with MaxStep option %! vref = [-1.205364552835178, 0.951542399860817]; %! opt = odeset ("MaxStep", 1e-3); %! sol = ode45 (@fpol, [0 -2], [2 0], opt); %! assert ([abs(sol.x(8)-sol.x(7))], [1e-3], 1e-3); -%! assert ([sol.x(end), sol.y(end,:)], [-2, vref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [-2; vref'], 1e-3); %!test # AbsTol option %! opt = odeset ("AbsTol", 1e-5); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # AbsTol and RelTol option %! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # RelTol and NormControl option -- higher accuracy %! opt = odeset ("RelTol", 1e-8, "NormControl", "on"); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-5); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-5); %!test # Keeps initial values while integrating %! opt = odeset ("NonNegative", 2); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, 2, 0], 0.5); +%! assert ([sol.x(end); sol.y(:,end)], [2; 2; 0], 0.5); %!test # Details of OutputSel and Refine can't be tested %! opt = odeset ("OutputFcn", @fout, "OutputSel", 1, "Refine", 5); %! sol = ode45 (@fpol, [0 2], [2 0], opt); @@ -463,28 +467,41 @@ %!test # Mass option as function %! opt = odeset ("Mass", @fmas); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # Mass option as matrix %! opt = odeset ("Mass", eye (2,2)); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # Mass option as sparse matrix %! opt = odeset ("Mass", sparse (eye (2,2))); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # Mass option as function and sparse matrix %! opt = odeset ("Mass", @fmsa); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); %!test # Mass option as function and MStateDependence %! opt = odeset ("Mass", @fmas, "MStateDependence", "strong"); %! sol = ode45 (@fpol, [0 2], [2 0], opt); -%! assert ([sol.x(end), sol.y(end,:)], [2, fref], 1e-3); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-3); -## FIXME: Missing tests. -## test for InitialSlope option is missing -## test for MaxOrder option is missing -## test for MvPattern option is missing +## Note: The following options have no effect on this solver +## therefore it makes no sense to test them here: +## +## "BDF" +## "InitialSlope" +## "JPattern" +## "Jacobian" +## "MassSingular" +## "MaxOrder" +## "MvPattern" +## "Vectorized" + +%!test # Check that imaginary part of solution does not get inverted +%! sol = ode45 (@(x,y) 1, [0 1], 1i); +%! assert (imag (sol.y), ones (size (sol.y))) +%! [x, y] = ode45 (@(x,y) 1, [0 1], 1i); +%! assert (imag (y), ones (size (y))) %!error ode45 () %!error ode45 (1)
--- a/scripts/ode/odeset.m Sun Oct 23 13:57:00 2016 -0700 +++ b/scripts/ode/odeset.m Mon Oct 24 09:04:30 2016 -0700 @@ -60,6 +60,7 @@ ## ## @item BDF ## Use BDF formulas in implicit multistep methods. +## @strong{Note:} This option is not yet implemented. ## ## @item Events ## Event function. An event function must have the form @@ -103,6 +104,7 @@ ## @item MvPattern ## If the mass matrix is sparse and non-constant but maintains a ## constant sparsity pattern, specify the sparsity pattern. +## @strong{Note:} This option is not yet implemented. ## ## @item NonNegative ## Specify elements of the state vector that are expected to remain @@ -125,6 +127,7 @@ ## time step or also at intermediate time instances. The value should be ## a scalar indicating the number of equally spaced time points to use ## within each timestep at which to return output. +## @strong{Note:} This option is not yet implemented. ## ## @item RelTol ## Relative error tolerance.
--- a/scripts/ode/private/integrate_adaptive.m Sun Oct 23 13:57:00 2016 -0700 +++ b/scripts/ode/private/integrate_adaptive.m Mon Oct 24 09:04:30 2016 -0700 @@ -193,6 +193,9 @@ approxvals = interp1 ([t_old, t(t_caught), t_new], [x_old, x(:, t_caught), x_new] .', approxtime, "linear") .'; + if (isvector (approxvals)) + approxvals = approxvals.'; + endif if (! isempty (options.OutputSel)) approxvals = approxvals(options.OutputSel, :); endif @@ -241,6 +244,9 @@ approxvals = interp1 ([t_old, t_new], [x_old, x_new] .', approxtime, "linear") .'; + if (isvector (approxvals)) + approxvals = approxvals.'; + endif if (! isempty (options.OutputSel)) approxvals = approxvals(options.OutputSel, :); endif
--- a/scripts/ode/private/ode_event_handler.m Sun Oct 23 13:57:00 2016 -0700 +++ b/scripts/ode/private/ode_event_handler.m Mon Oct 24 09:04:30 2016 -0700 @@ -88,7 +88,7 @@ [evt, term, dir] = feval (inpargs{:}); ## We require that all return values be row vectors - evt = evt(:)'; term = term(:)'; dir = dir(:)'; + evt = evt(:).'; term = term(:).'; dir = dir(:).'; ## Check if one or more signs of the event has changed signum = (sign (evtold) != sign (evt)); @@ -115,7 +115,7 @@ ## calculate new values for the integration results, we do both by ## a linear interpolation tnew = t - evt(1,idx) * (t - told) / (evt(1,idx) - evtold(1,idx)); - ynew = (y - (t - tnew) * (y - yold) / (t - told))'; + ynew = (y - (t - tnew) * (y - yold) / (t - told)).'; retcell{3}(evtcnt,1) = tnew; retcell{4}(evtcnt,:) = ynew; evtcnt += 1; @@ -142,7 +142,7 @@ [evtold, ~, ~] = feval (inpargs{:}); ## We require that all return values be row vectors - evtold = evtold(:)'; told = t; yold = y; + evtold = evtold(:).'; told = t; yold = y; evtcnt = 1; retval = retcell = cell (1,4);
--- a/scripts/ode/private/odedefaults.m Sun Oct 23 13:57:00 2016 -0700 +++ b/scripts/ode/private/odedefaults.m Mon Oct 24 09:04:30 2016 -0700 @@ -47,6 +47,7 @@ "Stats", "off", "Vectorized", "off"); + defaults.InitialSlope = zeros (n,1); defaults.MaxStep = 0.1 * abs (tf -t0); persistent classes = struct ("AbsTol", {{"float"}}, @@ -95,5 +96,9 @@ "RelTol", {{"scalar", "positive", "real"}}, "Stats", {{"on", "off"}}, "Vectorized", {{"on", "off"}}); + + attributes.InitialSlope = {"real", "vector", "numel", n}; + attributes.OutputSel = {"vector", "integer", "positive", ">", 0, "<=", n}; + endfunction
--- a/scripts/ode/private/runge_kutta_45_dorpri.m Sun Oct 23 13:57:00 2016 -0700 +++ b/scripts/ode/private/runge_kutta_45_dorpri.m Mon Oct 24 09:04:30 2016 -0700 @@ -60,6 +60,9 @@ k_vals = [], t_next = t + dt) + ## Reference: Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (2008), + ## Solving ordinary differential equations I: Nonstiff problems, + ## Berlin, New York: Springer-Verlag, ISBN 978-3-540-56670-0 persistent a = [0 0 0 0 0 0; 1/5 0 0 0 0 0; 3/40 9/40 0 0 0 0; @@ -70,18 +73,13 @@ 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)]; s = t + dt * b; cc = dt * c; aa = dt * a; k = zeros (rows (x), 7); - if (! isempty (options)) # extra arguments for function evaluator + if (! isempty (options)) # extra arguments for function evaluator args = options.funarguments; else args = {};
--- a/scripts/ode/private/runge_kutta_interpolate.m Sun Oct 23 13:57:00 2016 -0700 +++ b/scripts/ode/private/runge_kutta_interpolate.m Mon Oct 24 09:04:30 2016 -0700 @@ -22,7 +22,7 @@ switch (order) case 1 - u_interp = interp1 (z, u', t, "linear"); + u_interp = interp1 (z, u.', t, "linear"); case 2 if (! isempty (k_vals)) @@ -35,15 +35,6 @@ case 3 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 - ## doesn't require a new function evaluation. - u_interp = dorpri_interpolation ([z(i-1) z(i)], - [u(:,i-1) u(:,i)], - k_vals, tspan(counter)); - #} case 5 ## ode45 with Dormand-Prince scheme: u_interp = hermite_quartic_interpolation (z, u, k_vals, t);