# HG changeset patch # User Fabio Cassini # Date 1501686868 -7200 # Wed Aug 02 17:14:28 2017 +0200 # Node ID 776c5839997ab54fb096442ac0ef8e8899c1d1f4 # Parent d9ca3f15f739bd2f6c38abec074ce882f950bb24 Added ODE solver ode23tb, modified ode23, ode45 and common files to follow the new standard *scripts/ode/ode23tb.m: main file of ode23tb *scripts/ode/private/runge_kutta_23tb.m: stepper of ode23tb *scripts/ode/private/AbsRel_norm.m: modified to follow new standard (x_est substituted with err_est *scripts/ode/private/integrate_adaptive.m: modified to follow new standard (x_est substituted with err_est) *scripts/ode/private/starting_stepsize.m: modified to follow new standard (x_est substituded with err_est) *scripts/ode/private/runge_kutta_23.m: modified to follow new standard (err computed inside the stepper) *scripts/ode/private/runge_kutta_45_dorpri.m: modified to follow new standard (err computed inside the stepper) diff -r d9ca3f15f739 -r 776c5839997a scripts/ode/ode23tb.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/ode23tb.m Wed Aug 02 17:14:28 2017 +0200 @@ -0,0 +1,514 @@ +## Copyright (C) 2017 Fabio Cassini +## +## 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 +## . + +## -*- texinfo -*- +## @deftypefn {} {[@var{t}, @var{y}] =} ode23tb (@var{fun}, @var{trange}, @var{init}) +## @deftypefnx {} {[@var{t}, @var{y}] =} ode23tb (@var{fun}, @var{trange}, @var{init}, @var{ode_opt}) +## @deftypefnx {} {[@var{t}, @var{y}, @var{te}, @var{ye}, @var{ie}] =} ode23tb (@dots{}) +## @deftypefnx {} {@var{solution} =} ode23tb (@dots{}) +## @deftypefnx {} {} ode23tb (@dots{}) +## +## Solve a set of stiff Ordinary Differential Equations (stiff ODEs) +## with the implicit @nospell{TR-BDF2} method of order 2. +## For the definition of this method see +## @url{https://math.la.asu.edu/~gardner/TRBDF.pdf}. +## +## @var{fun} is a function handle, inline function, or string containing the +## name of the function that defines the ODE: @code{y' = f(t,y)}. The function +## must accept two inputs where the first is time @var{t} and the second is a +## column vector of unknowns @var{y}. +## +## @var{trange} specifies the time interval over which the ODE will be +## evaluated. Typically, it is a two-element vector specifying the initial and +## final times (@code{[tinit, tfinal]}). If there are more than two elements +## then the solution will also be evaluated at these intermediate time +## instances. +## +## By default, @code{ode23tb} 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"}. +## +## @var{init} contains the initial value for the unknowns. If it is a row +## vector then the solution @var{y} will be a matrix in which each column is +## the solution for the corresponding initial value in @var{init}. +## +## The optional fourth argument @var{ode_opt} specifies non-default options to +## the ODE solver. It is a structure generated by @code{odeset}. +## +## An analytically computed Jacobian is needed to integrate the ODE. +## +## The function typically returns two outputs. Variable @var{t} is a +## column vector and contains the times where the solution was found. The +## output @var{y} is a matrix in which each column refers to a different +## unknown of the problem and each row corresponds to a time in @var{t}. +## +## The output can also be returned as a structure @var{solution} which +## 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. +## +## If no output arguments are requested, and no @code{OutputFcn} is +## specified in @var{ode_opt}, then the @code{OutputFcn} is set to +## @code{odeplot} and the results of the solver are plotted immediately. +## +## If using the @qcode{"Events"} option then three additional outputs may +## be returned. @var{te} holds the time when an Event function returned a +## zero. @var{ye} holds the value of the solution at time @var{te}. @var{ie} +## contains an index indicating which Event function was triggered in the case +## of multiple Event functions. +## +## Example: Solve the scaled @nospell{Robertson} problem +## +## @example +## @group +## frp = @@(@var{t},@var{y}) [-0.04 * @var{y}(1) + 1e4 * @var{y}(2) * @var{y}(3); ... +## 0.04 * @var{y}(1) - 1e4 * @var{y}(2) * @var{y}(3) - 3 * 1e7 * @var{y}(2)^2; ... +## 3 * 1e7 * @var{y}(2)^2]; +## (1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)]; +## +## jac = @@(@var{t},@var{y}) [-0.04, 1e4 * @var{y}(3), 1e4 * @var{y}(2); ... +##0.04, -1e4 * @var{y}(3) - 3 * 1e7 * 2 * @var{y}(2), -1e4 * @var{y}(2); ... +##0 , 3 * 1e7 * 2 * @var{y}(2), 0]; +## +## opt = odeset ("Jacobian", jac, "RelTol", 0.005, "AbsTol", 1e-10); +## +## [@var{t},@var{y}] = ode23tb (frp, [0, 4 * 1e7], [1; 0; 0], opt); +## @end group +## @end example +## @seealso{odeset, odeget, ode23} +## @end deftypefn + +function varargout = ode23tb (fun, trange, init, varargin) + + if (nargin < 3) + print_usage (); + endif + + order = 2; + solver = "ode23tb"; + + if (nargin >= 4) + if (! isstruct (varargin{1})) + ## varargin{1:len} are parameters for fun + odeopts = odeset (); + funarguments = varargin; + elseif (numel (varargin) > 1) + ## varargin{1} is an ODE options structure opt + odeopts = varargin{1}; + funarguments = {varargin{2:numel (varargin)}}; + else # if (isstruct (varargin{1})) + odeopts = varargin{1}; + funarguments = {}; + endif + else # nargin == 3 + odeopts = odeset (); + funarguments = {}; + endif + + if (! isnumeric (trange) || ! isvector (trange)) + error ("Octave:invalid-input-arg", + "ode23tb: TRANGE must be a numeric vector"); + endif + + if (numel (trange) < 2) + error ("Octave:invalid-input-arg", + "ode23tb: TRANGE must contain at least 2 elements"); + elseif (trange(2) == trange(1)) + error ("Octave:invalid-input-arg", + "ode23tb: invalid time span, TRANGE(1) == TRANGE(2)"); + else + direction = sign (trange(2) - trange(1)); + endif + trange = trange(:); + + if (! isnumeric (init) || ! isvector (init)) + error ("Octave:invalid-input-arg", + "ode23tb: INIT must be a numeric vector"); + endif + init = init(:); + + if (ischar (fun)) + try + fun = str2func (fun); + catch + warning (lasterr); + end_try_catch + endif + if (! isa (fun, "function_handle")) + error ("Octave:invalid-input-arg", + "ode23tb: 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. + [defaults, classes, attributes] = odedefaults (numel (init), + trange(1), trange(end)); + + persistent ode23tb_ignore_options = ... + {"BDF", "InitialSlope", "JPattern", "MassSingular", "MaxOrder", "MvPattern", "Vectorized"}; + + defaults = rmfield (defaults, ode23tb_ignore_options); + classes = rmfield (classes, ode23tb_ignore_options); + attributes = rmfield (attributes, ode23tb_ignore_options); + + odeopts = odemergeopts ("ode23tb", odeopts, defaults, classes, attributes); + + odeopts.funarguments = funarguments; + odeopts.direction = direction; + + if (! isempty (odeopts.NonNegative)) + if (isempty (odeopts.Mass)) + odeopts.havenonnegative = true; + else + odeopts.havenonnegative = false; + warning ("Octave:invalid-input-arg", + ['ode23tb: option "NonNegative" is ignored', ... + " when mass matrix is set\n"]); + endif + else + odeopts.havenonnegative = false; + endif + + if (isempty (odeopts.OutputFcn) && nargout == 0) + odeopts.OutputFcn = @odeplot; + odeopts.haveoutputfunction = true; + else + odeopts.haveoutputfunction = ! isempty (odeopts.OutputFcn); + endif + + 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 + + if (! isempty (odeopts.Mass) && isnumeric (odeopts.Mass)) + havemasshandle = false; + mass = odeopts.Mass; # constant mass + elseif (isa (odeopts.Mass, "function_handle")) + havemasshandle = true; # mass defined by a function handle + else # no mass matrix - creating a diag-matrix of ones for mass + havemasshandle = false; # mass = diag (ones (length (init), 1), 0); + endif + + if (isempty (odeopts.Jacobian)) + error ("Octave:invalid-input-arg", + "ode23tb: An analytical Jacobian must be passed"); + elseif (! is_function_handle (odeopts.Jacobian)) + error ("Octave:invalid-input-arg", + "ode23tb: Jacobian must be a function handle"); + endif + + ## Starting the initialization of the core solver ode23tb + + if (havemasshandle) # Handle only the dynamic mass matrix, + if (! strcmp (odeopts.MStateDependence, "none")) + ## 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{:}); + else + mass = @(t) odeopts.Mass (t, odeopts.funarguments{:}); + fun = @(t,x) mass (t, odeopts.funarguments{:}) ... + \ fun (t, x, odeopts.funarguments{:}); + 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_23tb, + order, fun, trange, init, odeopts); + + ## Postprocessing, do whatever when terminating integration algorithm + if (odeopts.haveoutputfunction) # Cleanup plotter + feval (odeopts.OutputFcn, [], [], "done", odeopts.funarguments{:}); + endif + if (! isempty (odeopts.Events)) # Cleanup event function handling + ode_event_handler (odeopts.Events, solution.t(end), + solution.x(end,:).', "done", odeopts.funarguments{:}); + endif + + ## Print additional information if option Stats is set + if (strcmpi (odeopts.Stats, "on")) + nsteps = solution.cntloop; # cntloop from 2..end + nfailed = solution.cntcycles - nsteps; # cntcycl from 1..end + + printf ("Number of successful steps: %d\n", nsteps); + printf ("Number of total failed attempts: %d\n", nfailed); + endif + + if (nargout == 2) + 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 (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 + varargout{1}.xe = solution.event{3}; # Time info when an event occurred + varargout{1}.ye = solution.event{4}; # Results when an event occurred + endif + if (strcmpi (odeopts.Stats, "on")) + varargout{1}.stats = struct (); + varargout{1}.stats.nsteps = nsteps; + varargout{1}.stats.nfailed = nfailed; + endif + elseif (nargout > 2) + varargout = cell (1,5); + varargout{1} = solution.t; + varargout{2} = solution.x; + if (! isempty (odeopts.Events)) + varargout{3} = solution.event{3}; # Time info when an event occurred + varargout{4} = solution.event{4}; # Results when an event occurred + varargout{5} = solution.event{2}; # Index info which event occurred + endif + endif + +endfunction + + +%!demo +%! ## Demonstrate convergence order for ode23tb +%! tol = 1e-5 ./ 10.^[0:8]; +%! for i = 1 : numel (tol) +%! opt = odeset ("RelTol", tol(i), "AbsTol", realmin, "Jacobian", @(t,y) -1); +%! [t, y] = ode23tb (@(t, y) -y, [0, 1], 1, opt); +%! h(i) = 1 / (numel (t) - 1); +%! err(i) = norm (y .* exp (t) - 1, Inf); +%! endfor +%! +%! ## Estimate order visually +%! loglog (h, tol, "-ob", +%! h, err, "-b", +%! h, (h/h(end)) .^ 2 .* tol(end), "k--", +%! h, (h/h(end)) .^ 3 .* tol(end), "k-"); +%! axis tight +%! xlabel ("h"); +%! ylabel ("err(h)"); +%! title ("Convergence plot for ode23tb"); +%! legend ("imposed tolerance", "ode23tb (relative) error", +%! "order 2", "order 3", "location", "northwest"); +%! +%! ## Estimate order numerically +%! p = diff (log (err)) ./ diff (log (h)) + +## We are using the "Van der Pol" implementation 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 ODE +%! ydot = [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%!endfunction +%!function jac = jpol (t, y) # The Jacobian +%! jac = [0, 1; -2 * y(2) * y(1) - 1, (1 - y(1)^2)]; +%!endfunction +%!function ref = fref () # The computed reference sol +%! ref = [0.32331666704577, -1.83297456798624]; +%!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); # 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); # does not seem to be exact +%!endfunction +%!function mas = fmas (t, y, varargin) +%! 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 +%!endfunction +%!function out = fout (t, y, flag, varargin) +%! out = false; +%! if (strcmp (flag, "init")) +%! if (! isequal (size (t), [2, 1])) +%! error ('fout: step "init"'); +%! endif +%! elseif (isempty (flag)) +%! if (! isequal (size (t), [1, 1])) +%! error ('fout: step "calc"'); +%! endif +%! elseif (strcmp (flag, "done")) +%! if (! isempty (t)) +%! warning ('fout: step "done"'); +%! endif +%! else +%! error ("fout: invalid flag <%s>", flag); +%! endif +%!endfunction +%! +%!test # two output arguments +%! opt = odeset ("Jacobian", @jpol); +%! [t, y] = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); +%!test # anonymous function instead of real function +%! opt = odeset ("Jacobian", @jpol); +%! fvdp = @(t,y) [y(2); (1 - y(1)^2) * y(2) - y(1)]; +%! [t, y] = ode23tb (fvdp, [0 2], [2 0], opt); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); +%!test # extra input arguments passed through +%! opt = odeset ("Jacobian", @jpol); +%! [t, y] = ode23tb (@fpol, [0 2], [2 0], opt, 12, 13, "KL"); +%! assert ([t(end), y(end,:)], [2, fref], 1e-2); +%!test # Solve another anonymous function below zero +%! ref = [0, 14.77810590694212]; +%! opt = odeset ("Jacobian", @(t,y) 1); +%! [t, y] = ode23tb (@(t,y) y, [-2 0], 2, opt); +%! assert ([t(end), y(end,:)], ref, 1e-1); +%!test # InitialStep option +%! opt = odeset ("InitialStep", 1e-8, "Jacobian", @jpol); +%! [t, y] = ode23tb (@fpol, [0 0.2], [2 0], opt); +%! assert ([t(2)-t(1)], [1e-8], 1e-9); +%!test # MaxStep option +%! opt = odeset ("MaxStep", 1e-3, "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 0.2], [2 0], opt); +%! assert ([sol.x(5)-sol.x(4)], [1e-3], 1e-4); +%!test # Solve in backward direction starting at t=0 +%! ref = [-1.205364552835178, 0.951542399860817]; +%! opt = odeset ("Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 -2], [2 0], opt); +%! 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]; +%! opt = odeset ("Jacobian", @jpol); +%! sol = ode23tb (@fpol, [2 0 -2], fref, opt); +%! assert ([sol.x(end); sol.y(:,end)], [-2; ref'], 2e-1); +%!test # Solve another anonymous function in backward direction +%! ref = [-1, 0.367879437558975]; +%! opt = odeset ("Jacobian", @(t,y) 1); +%! sol = ode23tb (@(t,y) y, [0 -1], 1, opt); +%! assert ([sol.x(end); sol.y(:,end)], ref', 1e-2); +%!test # Solve another anonymous function below zero +%! ref = [0, 14.77810590694212]; +%! opt = odeset ("Jacobian", @(t,y) 1); +%! sol = ode23tb (@(t,y) y, [-2 0], 2, opt); +%! assert ([sol.x(end); sol.y(:,end)], ref', 1e-1); +%!test # Solve in backward direction starting at t=0 with MaxStep option +%! ref = [-1.205364552835178, 0.951542399860817]; +%! opt = odeset ("MaxStep", 1e-3, "Jacobian", @jpol); +%! sol = ode23tb (@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); +%!test # AbsTol option +%! opt = odeset ("AbsTol", 1e-5, "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-2); +%!test # AbsTol and RelTol option +%! opt = odeset ("AbsTol", 1e-8, "RelTol", 1e-8, "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-2); +%!test # RelTol and NormControl option -- higher accuracy +%! opt = odeset ("RelTol", 1e-8, "NormControl", "on", "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-4); +%!test # Keeps initial values while integrating +%! opt = odeset ("NonNegative", 2, "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! 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, "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%!test # Stats must add further elements in sol +%! opt = odeset ("Stats", "on", "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert (isfield (sol, "stats")); +%! assert (isfield (sol.stats, "nsteps")); +%!test # Events option add further elements in sol +%! opt = odeset ("Events", @feve, "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 10], [2 0], opt); +%! assert (isfield (sol, "ie")); +%! assert (sol.ie(1), 2); +%! assert (isfield (sol, "xe")); +%! assert (isfield (sol, "ye")); +%!test # Events option, now stop integration +%! warning ("off", "integrate_adaptive:unexpected_termination", "local"); +%! opt = odeset ("Events", @fevn, "NormControl", "on", "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 10], [2 0], opt); +%! assert ([sol.ie, sol.xe, sol.ye], +%! [2.0, 2.496110, -0.830550, -2.677589], .5e-1); +%!test # Events option, five output arguments +%! warning ("off", "integrate_adaptive:unexpected_termination", "local"); +%! opt = odeset ("Events", @fevn, "NormControl", "on", "Jacobian", @jpol); +%! [t, y, vxe, ye, vie] = ode23tb (@fpol, [0 10], [2 0], opt); +%! assert ([vie, vxe, ye], [2.0, 2.496110, -0.830550, -2.677589], 1e-1); +%!test # Mass option as function +%! opt = odeset ("Mass", @fmas, "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-2); +%!test # Mass option as matrix +%! opt = odeset ("Mass", eye (2,2), "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-2); +%!test # Mass option as sparse matrix +%! opt = odeset ("Mass", sparse (eye (2,2)), "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-2); +%!test # Mass option as function and sparse matrix +%! opt = odeset ("Mass", @fmsa, "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-2); +%!test # Mass option as function and MStateDependence +%! opt = odeset ("Mass", @fmas, "MStateDependence", "strong", "Jacobian", @jpol); +%! sol = ode23tb (@fpol, [0 2], [2 0], opt); +%! assert ([sol.x(end); sol.y(:,end)], [2; fref'], 1e-2); + +## Note: The following options have no effect on this solver +## therefore it makes no sense to test them here: +## +## "BDF" +## "InitialSlope" +## "JPattern" +## "MassSingular" +## "MaxOrder" +## "MvPattern" +## "Vectorized" + +%!test # Check that imaginary part of solution does not get inverted +%! opt = odeset ("Jacobian", @(x,y) 0); +%! sol = ode23tb (@(x,y) 1, [0 1], 1i, opt); +%! assert (imag (sol.y), ones (size (sol.y))) +%! [x, y] = ode23tb (@(x,y) 1, [0 1], 1i, opt); +%! assert (imag (y), ones (size (y))) + +## Test input validation +%!error ode23tb () +%!error ode23tb (1) +%!error ode23tb (1,2) +%!error ode23tb (@fpol, {[0 25]}, [3 15 1]) +%!error ode23tb (@fpol, [0 25; 25 0], [3 15 1]) +%!error ode23tb (@fpol, [1], [3 15 1]) +%!error ode23tb (@fpol, [1 1], [3 15 1]) +%!error ode23tb (@fpol, [0 25], {[3 15 1]}) +%!error ode23tb (@fpol, [0 25], [3 15 1; 3 15 1]) +%!error ode23tb (1, [0 25], [3 15 1]) diff -r d9ca3f15f739 -r 776c5839997a scripts/ode/private/AbsRel_norm.m --- a/scripts/ode/private/AbsRel_norm.m Mon Jul 31 16:25:19 2017 -0700 +++ b/scripts/ode/private/AbsRel_norm.m Wed Aug 02 17:14:28 2017 +0200 @@ -18,18 +18,19 @@ ## . ## -*- texinfo -*- -## @deftypefn {} {retval =} AbsRel_norm (@var{x}, @var{x_old}, @var{AbsTol}, @var{RelTol}, @var{normcontrol}, @var{y}) +## @deftypefn {} {retval =} AbsRel_norm (@var{x}, @var{x_old}, @var{err_est}, @var{AbsTol}, @var{RelTol}, @var{normcontrol}) ## Undocumented internal function. ## @end deftypefn -function retval = AbsRel_norm (x, x_old, AbsTol, RelTol, normcontrol, y = zeros (size (x))) +function retval = AbsRel_norm (x, x_old, err_est, AbsTol, RelTol, normcontrol) if (normcontrol) sc = max (AbsTol, RelTol * max (sqrt (sumsq (x)), sqrt (sumsq (x_old)))); - retval = sqrt (sumsq ((x - y))) / sc; + retval = sqrt (sumsq (err_est)) / sc; else sc = max (AbsTol, RelTol .* max (abs (x), abs (x_old))); - retval = max (abs (x - y) ./ sc); + retval = max (abs (err_est) ./ sc); endif + endfunction diff -r d9ca3f15f739 -r 776c5839997a scripts/ode/private/integrate_adaptive.m --- a/scripts/ode/private/integrate_adaptive.m Mon Jul 31 16:25:19 2017 -0700 +++ b/scripts/ode/private/integrate_adaptive.m Wed Aug 02 17:14:28 2017 +0200 @@ -5,13 +5,13 @@ ## ## 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. +## 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. +## 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 @@ -101,9 +101,7 @@ endif ## Initialize the EventFcn - have_EventFcn = false; if (! isempty (options.Events)) - have_EventFcn = true; ode_event_handler (options.Events, tspan(1), x, "init", options.funarguments{:}); endif @@ -118,7 +116,6 @@ solution.unhandledtermination = true; ireject = 0; - NormControl = strcmp (options.NormControl, "on"); k_vals = []; iout = istep = 1; @@ -126,19 +123,17 @@ ## 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); + + [t_new, x_new, err, new_k_vals] = ... + stepper (func, t_old, x_old, dt, options, k_vals, t_new); solution.cntcycles += 1; - if (options.havenonnegative) + if (options.havenonnegative) % TO DO: TO BE REWRITTEN x_new(nn, end) = abs (x_new(nn, end)); - x_est(nn, end) = abs (x_est(nn, end)); + % x_est(nn, end) = abs (x_est(nn, end)); endif - err = AbsRel_norm (x_new, x_old, options.AbsTol, options.RelTol, - NormControl, x_est); - ## Accept solution only if err <= 1.0 if (err <= 1) @@ -164,7 +159,7 @@ ## Call Events function only if a valid result has been found. ## Stop integration if eventbreak is true. - if (have_EventFcn) + if (! isempty (options.Events)) break_loop = false; for idenseout = 1:numel (t_caught) id = t_caught(idenseout); @@ -226,7 +221,7 @@ ## Call Events function only if a valid result has been found. ## Stop integration if eventbreak is true. - if (have_EventFcn) + if (! isempty (options.Events)) solution.event = ... ode_event_handler (options.Events, t(istep), x(:, istep), [], options.funarguments{:}); diff -r d9ca3f15f739 -r 776c5839997a scripts/ode/private/runge_kutta_23.m --- a/scripts/ode/private/runge_kutta_23.m Mon Jul 31 16:25:19 2017 -0700 +++ b/scripts/ode/private/runge_kutta_23.m Wed Aug 02 17:14:28 2017 +0200 @@ -22,8 +22,8 @@ ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23 (@var{fun}, @var{t}, @var{x}, @var{dt}, @var{options}) ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23 (@var{fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals}) ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23 (@var{fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals}, @var{t_next}) -## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{x_est}] =} runge_kutta_23 (@dots{}) -## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{x_est}, @var{k_vals_out}] =} runge_kutta_23 (@dots{}) +## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{err}] =} runge_kutta_23 (@dots{}) +## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{err}, @var{k}] =} runge_kutta_23 (@dots{}) ## ## This function can be used to integrate a system of ODEs with a given initial ## condition @var{x} from @var{t} to @var{t+dt}, with the Bogacki-Shampine @@ -40,27 +40,27 @@ ## ## @var{dt} is the timestep, that is the length of the integration interval. ## -## The optional fourth argument @var{options} specifies options for the ODE +## The optional fifth argument @var{options} specifies options for the ODE ## solver. It is a structure generated by @code{odeset}. In particular it ## contains the field @var{funarguments} with the optional arguments to be used ## in the evaluation of @var{fun}. ## -## The optional fifth argument @var{k_vals_in} contains the Runge-Kutta +## The optional sixth argument @var{k_vals} contains the Runge-Kutta ## evaluations of the previous step to use in a FSAL scheme. ## -## The optional sixth argument @var{t_next} (@code{t_next = t + dt}) specifies -## the end of the integration interval. The output @var{x_next} s the higher +## The optional seventh argument @var{t_next} (@code{t_next = t + dt}) specifies +## the end of the integration interval. The output @var{x_next} is the higher ## order computed solution at time @var{t_next} (local extrapolation is ## performed). ## -## Optionally the functions can also return @var{x_est}, a lower order solution -## for the estimation of the error, and @var{k_vals_out}, a matrix containing -## the Runge-Kutta evaluations to use in a FSAL scheme or for dense output. +## Optionally the functions can also return @var{err}, an estimation of the +## error, and @var{k}, a matrix containing the Runge-Kutta evaluations to use +## in a FSAL scheme or for dense output. ## ## @seealso{runge_kutta_45_dorpri} ## @end deftypefn -function [t_next, x_next, x_est, k] = runge_kutta_23 (fun, t, x, dt, +function [t_next, x_next, err, k] = runge_kutta_23 (fun, t, x, dt, options = [], k_vals = [], t_next = t + dt) @@ -102,6 +102,9 @@ k(:,4) = feval (fun, t_next, x_next, args{:}); cc_prime = dt * c_prime; x_est = x + k * cc_prime(:); + err_est = x_next - x_est; + err = AbsRel_norm (x_next, x, err_est, options.AbsTol, options.RelTol, ... + strcmpi (options.NormControl, "on")); endif endfunction diff -r d9ca3f15f739 -r 776c5839997a scripts/ode/private/runge_kutta_23tb.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/ode/private/runge_kutta_23tb.m Wed Aug 02 17:14:28 2017 +0200 @@ -0,0 +1,207 @@ +## Copyright (C) 2017 Fabio Cassini +## +## 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 +## . + +## -*- texinfo -*- +## @deftypefn {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23tb (@var{f}, @var{t}, @var{x}, @var{dt}) +## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23tb (@var{f}, @var{t}, @var{x}, @var{dt}, @var{options}) +## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23tb (@var{f}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals}) +## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_23tb (@var{f}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals}, @var{t_next}) +## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{err}] =} runge_kutta_23tb (@dots{}) +## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{err}, @var{k}] =} runge_kutta_23tb (@dots{}) +## +## This function can be used to integrate a system of stiff ODEs with a given initial +## condition @var{x} from @var{t} to @var{t+dt}, with the TR-BDF2 +## method of order 2. For the definition of this method see +## @url{https://math.la.asu.edu/~gardner/TRBDF.pdf}. +## +## @var{f} is a function handle that defines the ODE: @code{y' = f(tau,y)}. +## The function must accept two inputs where the first is time @var{tau} and +## the second is a column vector of unknowns @var{y}. +## +## @var{t} is the first extreme of integration interval. +## +## @var{x} is the initial condition of the system.. +## +## @var{dt} is the timestep, that is the length of the integration interval. +## +## The optional fifth argument @var{options} specifies options for the ODE +## solver. It is a structure generated by @code{odeset}. In particular it +## contains the field @var{funarguments} with the optional arguments to be used +## in the evaluation of @var{fun}. +## +## The optional sixth argument @var{k_vals} contains the Runge-Kutta +## evaluations of the previous step to use in a FSAL scheme. +## +## The optional seventh argument @var{t_next} (@code{t_next = t + dt}) specifies +## the end of the integration interval. +## +## The output @var{x_next} is the lower order computed solution at time @var{t_next}. +## +## The output @var{err} is an estimation of the error. +## +## The output @var{k} is a matrix containing the Runge-Kutta evaluations to use +## in a FSAL scheme or for dense output. +## +## @seealso{runge_kutta_23} +## @end deftypefn + +function [t_next, x_next, err, k] = runge_kutta_23tb (f, t, x, dt, + options = [], + k_vals = [], + t_next = t + dt) + + persistent gamma = 2 - sqrt (2); + persistent d = gamma / 2; + persistent w = sqrt (2) / 4; + + persistent c = [0, gamma, 1]; + persistent b = [w, w, d]; + persistent bhat = [(1 - w) / 3, (3 * w + 1) / 3, d / 3]; + persistent A = [0, 0, 0; d, d, 0; w, w, d]; + + normcontrol = strcmp(options.NormControl, "on"); + abstol = options.AbsTol; + reltol = options.RelTol; + + maxit = 10; % first maximum number of iterations + newt1 = 0; % number of first newton iterations + newt2 = 0; % number of second newton iterations + + flagc1 = 1; + flagc2 = 1; + + w = length(x); + k = zeros(w,3); + x_next = zeros(w,1); + x_est = zeros(w,1); + + err = realmax; # the step will be rejected unless everything goes ok + + if (! isempty (options)) # extra arguments for function evaluator + args = options.funarguments; + else + args = {}; + endif + + % First step + + if (! isempty (k_vals)) # k values from previous step are passed + k(:,1) = k_vals(:,3); # FSAL property + else + k(:,1) = feval (f ,t, x, args{:}); + endif + + % Second step + + F = @(fi2) fi2 - feval (f ,t + gamma * dt, x + dt * (A(2,1) * k(:,1) + A(2,2) * fi2), args{:}); + + k(:,2) = k(:,1); + + I = speye (w); + jf = feval (options.Jacobian, t + c(2) * dt, x + dt * (A(2,1) * k(:,1) + A(2,2) * k(:,2))); + JF = I - dt * A(2,2) * jf; + [L, U, p] = lu (JF, "vector"); + + delta = - U \ (L \ (feval (F, (k(:,2))(p)))); + + con = (norm (delta, inf)) > (abstol + norm (k(:,2), inf) * reltol) / 10; + + while con && (newt1 < maxit) + newt1 += 1; + k(:,2) = k(:,2) + delta; + delta = - U \ (L \ (feval (F, (k(:,2))(p)))); + c1 = (norm (delta, inf)) > (abstol + norm (k(:,2), inf) * reltol) / 10; + endwhile + + if (con) % 10 iterations are not sufficient + if (norm(delta, inf)) < 10 * (abstol + norm (k(:,2), inf) * reltol) % newton is converging, update jf and go on with other 10 iters + jf = feval (options.Jacobian, t + c(2) * dt, x + dt * (A(2,1) * k(:,1) + A(2,2) * k(:,2))); + JF = I - dt * A(2,2) * jf; + [L, U, p] = lu (JF, "vector"); + + con = (norm (delta, inf)) > (abstol + norm (k(:,2), inf) * reltol) / 10; + while con && (newt1 < 2 * maxit) + newt1 += 1; + k(:,2) = k(:,2) + delta; + delta = - U \ (L \ (feval (F, (k(:,2))(p)))); + con = (norm (delta, inf)) > (abstol + norm (k(:,2), inf) * reltol) / 10; + endwhile + + if (con) % 20 iters are not sufficient, go out + flagc1 = 0; + else + k(:,2) = k(:,2) + delta; + endif + else % newton is not converging, go out + flagc1 = 0; + endif + else + k(:,2) = k(:,2) + delta; + endif + + if flagc1 == 1 % Compute the third step only if second passed + F = @(fi3) fi3 - feval (f, t + dt, x + dt * (A(3,1) * k(:,1) + A(3,2) * k(:,2) + A(3,3) * fi3), args{:}); + + k(:,3) = k(:,2); + delta = - U \ (L \ (feval (F, (k(:,3))(p)))); + + con = norm (delta, inf) > (abstol + norm (k(:,3), inf) * reltol) / 10; + + while con && (newt2 < maxit) + newt2 += 1; + k(:,3) = k(:,3) + delta; + delta = - U \ (L \ (feval (F, (k(:,3))(p)))); + con = norm (delta, inf) > (abstol + norm (k(:,3), inf) * reltol) / 10; + endwhile + + if (con) + if (norm(delta,inf)) < 10 * (abstol + norm (k(:,3), inf) * reltol) + jf = feval (options.Jacobian, t + dt, x + dt * (A(3,1) * k(:,1) + A(3,2) * k(:,2) + A(3,3) * k(:,3))); + JF = I - dt * A(3,3) * jf; + [L, U, p] = lu (JF, "vector"); + con = (norm (delta, inf)) > (abstol + norm (k(:,3), inf) * reltol) / 10; + + while con && (newt2 < 2 * maxit) + newt2 += 1; + k(:,3) = k(:,3) + delta; + delta = - U \ (L \ (feval (F, (k(:,2))(p)))); + con = (norm (delta, inf)) > (abstol + norm (k(:,3), inf) * reltol) / 10; + endwhile + + if (con) + flagc2 = 0; + else + k(:,3) = k(:,3) + delta; + endif + else + flagc2 = 0; + endif + else + k(:,3) = k(:,3) + delta; + endif + + if flagc2 == 1 % everything is ok + x_next = x + dt * (b(1) * k (:,1) + b(2) * k(:,2) + b(3) * k(:,3)); + x_est = x + dt * (bhat(1) * k(:,1) + bhat(2) * k(:,2) + bhat(3) * k(:,3)); + err_est = U \ (L \ ((x_next - x_est)(p))); + err = AbsRel_norm (x_next, x, err_est, abstol, reltol, normcontrol); + endif + + endif + +endfunction diff -r d9ca3f15f739 -r 776c5839997a scripts/ode/private/runge_kutta_45_dorpri.m --- a/scripts/ode/private/runge_kutta_45_dorpri.m Mon Jul 31 16:25:19 2017 -0700 +++ b/scripts/ode/private/runge_kutta_45_dorpri.m Wed Aug 02 17:14:28 2017 +0200 @@ -23,7 +23,7 @@ ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals}) ## @deftypefnx {} {[@var{t_next}, @var{x_next}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals}, @var{t_next}) ## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{x_est}] =} runge_kutta_45_dorpri (@dots{}) -## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{x_est}, @var{k_vals_out}] =} runge_kutta_45_dorpri (@dots{}) +## @deftypefnx {} {[@var{t_next}, @var{x_next}, @var{err}, @var{k}] =} runge_kutta_45_dorpri (@dots{}) ## ## This function can be used to integrate a system of ODEs with a given initial ## condition @var{x} from @var{t} to @var{t+dt} with the Dormand-Prince method. @@ -54,14 +54,13 @@ ## 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. +## Third output parameter is an 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 (fun, t, x, dt, +function [t_next, x_next, err, k] = runge_kutta_45_dorpri (fun, t, x, dt, options = [], k_vals = [], t_next = t + dt) @@ -113,6 +112,9 @@ k(:,7) = feval (fun, t_next, x_next, args{:}); cc_prime = dt * c_prime; x_est = x + k * cc_prime(:); + err_est = x_next - x_est; + err = AbsRel_norm (x_next, x, err_est, options.AbsTol, options.RelTol, ... + strcmpi (options.NormControl, "on")); endif endfunction diff -r d9ca3f15f739 -r 776c5839997a scripts/ode/private/starting_stepsize.m --- a/scripts/ode/private/starting_stepsize.m Mon Jul 31 16:25:19 2017 -0700 +++ b/scripts/ode/private/starting_stepsize.m Wed Aug 02 17:14:28 2017 +0200 @@ -40,14 +40,14 @@ args = {}) ## compute norm of initial conditions - d0 = AbsRel_norm (x0, x0, AbsTol, RelTol, normcontrol); + d0 = AbsRel_norm (x0, 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, y, AbsTol, RelTol, normcontrol); if (d0 < 1e-5 || d1 < 1e-5) h0 = 1e-6; @@ -64,7 +64,7 @@ yh = yh{1}; endif d2 = (1 / h0) * ... - AbsRel_norm (yh - y, yh - y, AbsTol, RelTol, normcontrol); + AbsRel_norm (yh - y, yh - y, yh - y, AbsTol, RelTol, normcontrol); if (max (d1, d2) <= 1e-15) h1 = max (1e-6, h0 * 1e-3);