Mercurial > octave
changeset 20621:b92f8e148936
maint: Continued clean-up of functions in ode/private dir.
* AbsRel_Norm.m: Use retval as return variable. Use sumsq() rather than explicit
squaring of vector and sum(). Combine multiple lines where possible.
* integrate_adaptive.m: Rewrite docstring. Only call starting_stepsize() if
InitialStep option is empty.
* integrate_const.m: Rewrite docstring. Remove useless commented out code.
Combine multiple lines where possible.
* integrate_n_steps.m: Rewrite docstring. Remove useless commented out code.
Combine multiple lines where possible.
* kahan.m: Remove excessive 4-space indentation, use 2-space indentation.
* ode_rk_interpolate.m: Use parentheses around condition for switch stmt.
Combine multiple lines where possible.
* ode_struct_value_check.m: Remove comma from Copyright statement that Octave
doesn't use.
* odepkg_event_handle.m: Remove comma from Copyright statement that Octave
doesn't use.
* odepkg_structure_check.m: Remove comma from Copyright statement that Octave
doesn't use.
* runge_kutta_45_dorpri.m: Remove comma from Copyright statement that Octave
doesn't use. Improve docstring. Match variable names in documentation to
those in code.
* starting_stepsize.m: Rewrite docstring. Use spaces between function name
and opening parenthesis.
| author | Rik <rik@octave.org> |
|---|---|
| date | Wed, 14 Oct 2015 10:35:53 -0700 |
| parents | b5d2ca6a690c |
| children | caa5dabc7913 |
| files | scripts/ode/private/AbsRel_Norm.m scripts/ode/private/integrate_adaptive.m scripts/ode/private/integrate_const.m scripts/ode/private/integrate_n_steps.m scripts/ode/private/kahan.m scripts/ode/private/ode_rk_interpolate.m scripts/ode/private/ode_struct_value_check.m scripts/ode/private/odepkg_event_handle.m scripts/ode/private/odepkg_structure_check.m scripts/ode/private/runge_kutta_45_dorpri.m scripts/ode/private/starting_stepsize.m |
| diffstat | 11 files changed, 234 insertions(+), 246 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ode/private/AbsRel_Norm.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/AbsRel_Norm.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2014, Jacopo Corno <jacopo.corno@gmail.com> -## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2014-15 Jacopo Corno <jacopo.corno@gmail.com> +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> ## ## This file is part of Octave. ## @@ -17,34 +17,28 @@ ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. -function res = AbsRel_Norm (x, x_old, AbsTol, RelTol, normcontrol, y) +function retval = AbsRel_Norm (x, x_old, AbsTol, RelTol, normcontrol, y) n = length (x); if (nargin == 5) y = zeros (size (x)); - elseif (nargin != 6) - error ("OdePkg:InvalidArgument", - "invalid number of input arguments"); endif - if (length (x_old) != n - || length (y) != n) - error ("OdePkg:InvalidArgument", - "invalid dimensions of input arguments"); + if (length (x_old) != n || length (y) != n) + error ("OdePkg:InvalidArgument", "invalid dimensions of input arguments"); endif if ((length (AbsTol) != 1 && length (AbsTol) != n) || (length (RelTol) != 1 && length (RelTol) != n)) - error ("OdePkg:InvalidArgument", - "invalid dimensions of input arguments"); + error ("OdePkg:InvalidArgument", "invalid dimensions of input arguments"); endif sc = AbsTol + max (abs (x), abs (x_old)) .* RelTol; if (normcontrol) - res = max (abs (x - y) ./ sc); + retval = max (abs (x - y) ./ sc); else - res = sqrt ((1 / n) * sum (((x - y) ./ sc).^2)); + retval = sqrt ((1 / n) * sumsq ((x - y) ./ sc)); endif endfunction
--- a/scripts/ode/private/integrate_adaptive.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/integrate_adaptive.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2015, Carlo de Falco -## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2015 Carlo de Falco +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> ## ## This file is part of Octave. ## @@ -18,28 +18,29 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{t}, @var{y}] =} integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@fun}, @var{tspan}, @var{x0}, @var{options}) +## @deftypefn {Function File} {@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 +## integrate the set of ODEs on the interval @var{[t0, t1]} with an ## adaptive timestep. ## -## This function must be called with two output arguments: @var{t} and @var{y}. -## Variable @var{t} is a column vector and contains the time stamps, instead +## The function returns a structure @var{solution} with two fieldss: @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 rows number is the same of @var{t} rows number so -## that each row of @var{y} contains the values of all unknowns at the time -## value contained in the corresponding row in @var{t}. +## of the problem and the row number is the same as the @var{t} row number. +## Thus, each row of the matrix @var{y} contains the values of all unknowns at +## the time value contained in the corresponding row in @var{t}. ## -## The first input argument must be a function_handle or an inline function -## representing the stepper, that is the function responsible for step-by-step +## The first input argument must be a function handle or inline function +## representing the stepper, i.e., the function responsible for step-by-step ## integration. This function discriminates one method from the others. ## ## The second input argument is the order of the stepper. It is needed ## to compute the adaptive timesteps. ## -## The third input argument is a function_handle or an inline function that -## defines the set of ODE: +## The third input argument is a function handle or inline function that +## defines the ODE: +## ## @ifhtml ## @example ## @math{y' = f(t,y)} @@ -49,9 +50,9 @@ ## @math{y' = f(t,y)}. ## @end ifnothtml ## -## The fourth input argument is the time vector which defines integration -## interval, that is @var{[tspan(1),tspan(end)]} and all the intermediate -## elements are taken as times at which the solution is required. +## The fourth input argument is the time vector which defines the integration +## interval, i.e, @var{[tspan(1), tspan(end)]} and all intermediate elements +## are taken as times at which the solution is required. ## ## The fifth argument represents the initial conditions for the ODEs and the ## last input argument contains some options that may be needed for the stepper. @@ -60,17 +61,20 @@ ## ## @seealso{integrate_const, integrate_n_steps} -function solution = integrate_adaptive (stepper, order, func, tspan, x0, options) +function solution = integrate_adaptive (stepper, order, func, tspan, x0, + options) fixed_times = numel (tspan) > 2; t_new = t_old = t = tspan(1); x_new = x_old = x = x0(:); - ## get first initial timestep - dt = starting_stepsize (order, func, t, x, options.AbsTol, - options.RelTol, options.vnormcontrol); - dt = odeget (options, "InitialStep", dt, "fast_not_empty"); + ## Get first initial timestep + dt = odeget (options, "InitialStep", [], "fast"); + if (isempty (dt)) + dt = starting_stepsize (order, func, t, x, options.AbsTol, options.RelTol, + options.vnormcontrol); + endif dir = odeget (options, "vdirection", [], "fast"); dt = dir * min (abs (dt), options.MaxStep); @@ -116,7 +120,7 @@ ## 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.vcntcycles++; @@ -128,7 +132,7 @@ err = AbsRel_Norm (x_new, x_old, options.AbsTol, options.RelTol, options.vnormcontrol, x_est); - ## Accepted solution only if err <= 1.0 + ## Accept solution only if err <= 1.0 if (err <= 1) solution.vcntloop++; @@ -145,24 +149,23 @@ t(t_caught) = tspan(t_caught); iout = max (t_caught); x(:, t_caught) = ... - ode_rk_interpolate (order, [t_old t_new], [x_old x_new], - tspan(t_caught), new_k_vals, dt, - options.vfunarguments{:}); + ode_rk_interpolate (order, [t_old t_new], [x_old x_new], + tspan(t_caught), new_k_vals, dt, + options.vfunarguments{:}); istep++; + ## Call Events function only if a valid result has been found. + ## Stop integration if veventbreak is true. if (options.vhaveeventfunction) - ## Call event on each dense output timestep. - ## Stop integration if veventbreak is true break_loop = false; for idenseout = 1:numel (t_caught) id = t_caught(idenseout); td = t(id); solution.vevent = ... - odepkg_event_handle (options.Events, t(id), x(:, id), [], - options.vfunarguments{:}); - if (! isempty (solution.vevent{1}) - && solution.vevent{1} == 1) + odepkg_event_handle (options.Events, t(id), x(:, id), [], + options.vfunarguments{:}); + if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) t(id) = solution.vevent{3}(end); t = t(1:id); x(:, id) = solution.vevent{4}(end, :).'; @@ -177,8 +180,8 @@ endif endif - ## Call plot. Stop integration if plot function - ## returns true. + ## Call OutputFcn only if a valid result has been found. + ## Stop integration if function returns false. if (options.vhaveoutputfunction) vcnt = options.Refine + 1; vapproxtime = linspace (t_old, t_new, vcnt); @@ -207,23 +210,22 @@ x(:, istep) = x_new; iout = istep; - ## Call event handler on new timestep. - ## Stop integration if veventbreak is true + ## Call Events function only if a valid result has been found. + ## Stop integration if veventbreak is true. if (options.vhaveeventfunction) solution.vevent = ... - odepkg_event_handle (options.Events, t(istep), x(:, istep), [], - options.vfunarguments{:}); - if (! isempty (solution.vevent{1}) - && solution.vevent{1} == 1) - t(istep) = solution.vevent{3}(end); - x(:, istep) = solution.vevent{4}(end, :).'; - solution.vunhandledtermination = false; - break; - endif + odepkg_event_handle (options.Events, t(istep), x(:, istep), [], + options.vfunarguments{:}); + if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) + t(istep) = solution.vevent{3}(end); + x(:, istep) = solution.vevent{4}(end, :).'; + solution.vunhandledtermination = false; + break; + endif endif - ## Call plot. Stop integration if plot function - ## returns true. + ## Call OutputFcn only if a valid result has been found. + ## Stop integration if function returns false. if (options.vhaveoutputfunction) vcnt = options.Refine + 1; vapproxtime = linspace (t_old, t_new, vcnt); @@ -250,22 +252,22 @@ x_old = x_new; k_vals = new_k_vals; - solution.vcntloop = solution.vcntloop + 1; + solution.vcntloop += 1; vcntiter = 0; else ireject++; - ## Stop solving because in the last 5000 steps no successful valid + ## Stop solving because, in the last 5,000 steps, no successful valid ## value has been found - if (ireject >= 5000) - error (["integrate_adaptive: Solving has not been successful. ", - " The iterative integration loop exited at time", - " t = %f before endpoint at tend = %f was reached. ", - " This happened because the iterative integration loop", - " does not find a valid solution at this time", - " stamp. Try to reduce the value of 'InitialStep' and/or", + if (ireject >= 5_000) + error (["integrate_adaptive: Solving was not successful. ", ... + " The iterative integration loop exited at time", ... + " t = %f before the endpoint at tend = %f was reached. ", ... + " This happened because the iterative integration loop", ... + " did not find a valid solution at this time stamp. ", ... + " Try to reduce the value of 'InitialStep' and/or", ... " 'MaxStep' with the command 'odeset'.\n"], t_old, tspan(end)); endif @@ -273,7 +275,7 @@ endif ## Compute next timestep, formula taken from Hairer - err += eps; # avoid divisions by zero + err += eps; # avoid divisions by zero dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1)))); dt = dir * min (abs (dt), options.MaxStep); @@ -282,26 +284,24 @@ endwhile - ## Check if integration of the ode has been successful if (dir * t(end) < dir * tspan(end)) if (solution.vunhandledtermination == true) error ("integrate_adaptive:unexpected_termination", - [" Solving has not been successful. The iterative", ... - " integration loop exited at time t = %f ", ... - " before endpoint at tend = %f was reached. This may", ... - " happen if the stepsize grows too small. ", ... + [" Solving was not successful. ", ... + " The iterative integration loop exited at time", ... + " t = %f before the endpoint at tend = %f was reached. ", ... + " This may happen if the stepsize becomes too small. ", ... " Try to reduce the value of 'InitialStep'", ... " and/or 'MaxStep' with the command 'odeset'.\n"], t(end), tspan(end)); else warning ("integrate_adaptive:unexpected_termination", - ["Solver has been stopped by a call of 'break' ", ... - " in the main iteration loop at time t = %f before", ... - " endpoint at tend = %f was reached. This may", ... - " happen because the @odeplot function", ... - " returned 'true' or the @event function returned", ... - " 'true'.\n"], + ["Solver was stopped by a call of 'break'", ... + " in the main iteration loop at time ", ... + " t = %f before the endpoint at tend = %f was reached. ", ... + " This may happen because the @odeplot function ", ... + " returned 'true' or the @event function returned 'true'.\n"], t(end), tspan(end)); endif endif
--- a/scripts/ode/private/integrate_const.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/integrate_const.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,4 +1,4 @@ -## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> ## ## This file is part of Octave. ## @@ -17,28 +17,28 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{t}, @var{y}] =} integrate_const (@var{@@stepper}, @var{@@fun}, @var{tspan}, @var{x0}, @var{dt}, @var{options}) +## @deftypefn {Function File} {@var{solution} =} integrate_const (@var{@@stepper}, @var{@@func}, @var{tspan}, @var{x0}, @var{dt}, @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 a constant ## timestep @var{dt}. ## -## This function must be called with two output arguments: @var{t} and @var{y}. -## Variable @var{t} is a column vector and contains the time stamps, instead -## @var{y} is a matrix in which each column refers to a different unknown of -## the problem and the rows number is the same of @var{t} rows number so that -## each row of @var{y} contains the values of all unknowns at the time value -## contained in the corresponding row in @var{t}. +## The function returns a structure @var{solution} with two fieldss: @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. +## Thus, each row of the matrix @var{y} contains the values of all unknowns at +## the time value contained in the corresponding row in @var{t}. ## -## The first input argument must be a function_handle or an inline function -## representing the stepper, that is the function responsible for step-by-step +## The first input argument must be a function handle or inline function +## representing the stepper, i.e., the function responsible for step-by-step ## integration. This function discriminates one method from the others. ## ## The second input argument is the order of the stepper. It is needed to ## compute the adaptive timesteps. ## -## The third input argument is a function_handle or an inline function that -## defines the set of ODE: +## The third input argument is a function handle or inline function that +## defines the ODE: ## ## @ifhtml ## @example @@ -49,9 +49,9 @@ ## @math{y' = f(t,y)}. ## @end ifnothtml ## -## The third input argument is the time vector which defines integration -## interval, that is @var{[tspan(1),tspan(end)]} and all the intermediate -## elements are taken as times at which the solution is required. +## The fourth input argument is the time vector which defines the integration +## interval, i.e, @var{[tspan(1), tspan(end)]} and all intermediate elements +## are taken as times at which the solution is required. ## ## The fourth argument contains the initial conditions for the ODEs. ## @@ -101,8 +101,6 @@ solution.vcntloop = 2; solution.vcntcycles = 1; - #vu = vinit; - #vk = vu.' * zeros(1,6); vcntiter = 0; solution.vunhandledtermination = true; solution.vcntsave = 2; @@ -168,8 +166,7 @@ endwhile ## if new time requested is not out of this interval - if ((counter <= k) - && vdirection * z(end) > vdirection * tspan(counter)) + if (counter <= k && vdirection * z(end) > vdirection * tspan(counter)) ## update the counter i++; else @@ -180,19 +177,17 @@ endwhile endif - x = [x,u(:,2:end)]; t = [t;z(2:end)]; - solution.vcntsave = solution.vcntsave + 1; - solution.vcntloop = solution.vcntloop + 1; + solution.vcntsave += 1; + solution.vcntloop += 1; vcntiter = 0; - ## Call plot only if a valid result has been found, therefore this - ## code fragment has moved here. Stop integration if plot function - ## returns false + ## Call OutputFcn only if a valid result has been found. + ## Stop integration if function returns false. if (options.vhaveoutputfunction) for vcnt = 0:options.Refine # Approximation between told and t - if (options.vhaverefine) # Do interpolation + if (options.vhaverefine) # Do interpolation vapproxtime = (vcnt + 1) / (options.Refine + 2); vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ... + (vapproxtime) * y(:,end); @@ -216,13 +211,12 @@ endif endif - ## Call event only if a valid result has been found, therefore this - ## code fragment has moved here. Stop integration if veventbreak is true + ## Call Events function only if a valid result has been found. + ## Stop integration if veventbreak is true. if (options.vhaveeventfunction) solution.vevent = odepkg_event_handle (options.Events, t(end), x(:,end), [], options.vfunarguments{:}); - if (! isempty (solution.vevent{1}) - && solution.vevent{1} == 1) + if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) t(solution.vcntloop-1,:) = solution.vevent{3}(end,:); x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)'; solution.vunhandledtermination = false; @@ -231,17 +225,18 @@ endif ## Update counters that count the number of iteration cycles - solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics - vcntiter = vcntiter + 1; # Needed to find iteration problems + solution.vcntcycles += 1; # Needed for cost statistics + vcntiter += 1; # Needed to find iteration problems - ## Stop solving because the last 1000 steps no successful valid + ## Stop solving because, in the last 5,000 steps, no successful valid ## value has been found - if (vcntiter >= 5000) - error (["Solving has not been successful. The iterative", - " integration loop exited at time t = %f before endpoint at", - " tend = %f was reached. This happened because the iterative", - " integration loop does not find a valid solution at this time", - " stamp. Try to reduce the value of 'InitialStep' and/or", + if (vcntiter >= 5_000) + error (["integrate_const: Solving was not successful. ", ... + " The iterative integration loop exited at time", ... + " t = %f before the endpoint at tend = %f was reached. ", ... + " This happened because the iterative integration loop", ... + " did not find a valid solution at this time stamp. ", ... + " Try to reduce the value of 'InitialStep' and/or", ... " 'MaxStep' with the command 'odeset'.\n"], s(end), tspan(end)); endif @@ -250,23 +245,26 @@ if (counter > k) j = length (z); endif + endwhile ## Check if integration of the ode has been successful if (vdirection * z(end) < vdirection * tspan(end)) if (solution.vunhandledtermination == true) - error ("OdePkg:InvalidArgument", - ["Solving has not been successful. The iterative integration" - " loop exited at time t = %f before endpoint at tend = %f was", - " reached. This may happen if the stepsize grows smaller than", - " defined in vminstepsize. Try to reduce the value of", - " 'InitialStep' and/or 'MaxStep' with the command 'odeset'.\n"], + error ("integrate_const:unexpected_termination", + [" Solving was not successful. ", ... + " The iterative integration loop exited at time", ... + " t = %f before the endpoint at tend = %f was reached. ", ... + " This may happen if the stepsize becomes too small. ", ... + " Try to reduce the value of 'InitialStep'", ... + " and/or 'MaxStep' with the command 'odeset'.\n"], z(end), tspan(end)); else - warning ("OdePkg:InvalidArgument", - ["Solver has been stopped by a call of 'break' in the main", - " iteration loop at time t = %f before endpoint at tend = %f", - " was reached. This may happen because the @odeplot function", + warning ("integrate_const:unexpected_termination", + ["Solver was stopped by a call of 'break'", ... + " in the main iteration loop at time ", ... + " t = %f before the endpoint at tend = %f was reached. ", ... + " This may happen because the @odeplot function ", ... " returned 'true' or the @event function returned 'true'.\n"], z(end), tspan(end)); endif
--- a/scripts/ode/private/integrate_n_steps.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/integrate_n_steps.m Wed Oct 14 10:35:53 2015 -0700 @@ -17,28 +17,28 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {[@var{t}, @var{y}] =} integrate_n_steps (@var{@@stepper}, @var{@@fun}, @var{t0}, @var{x0}, @var{dt}, @var{n}, @var{options}) +## @deftypefn {Function File} {@var{solution} =} integrate_n_steps (@var{@@stepper}, @var{@@func}, @var{t0}, @var{x0}, @var{dt}, @var{n}, @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,t0 + n*dt]} with a ## constant timestep dt and on a fixed number of steps. ## -## This function must be called with two output arguments: @var{t} and @var{y}. -## Variable @var{t} is a column vector and contains the time stamps, instead -## @var{y} is a matrix in which each column refers to a different unknown of -## the problem and the rows number is the same of @var{t} rows number so that -## each row of @var{y} contains the values of all unknowns at the time -## value contained in the corresponding row in @var{t}. +## The function returns a structure @var{solution} with two fieldss: @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. +## Thus, each row of the matrix @var{y} contains the values of all unknowns at +## the time value contained in the corresponding row in @var{t}. ## -## The first input argument must be a function_handle or an inline function -## representing the stepper, that is the function responsible for step-by-step +## The first input argument must be a function handle or inline function +## representing the stepper, i.e., the function responsible for step-by-step ## integration. This function discriminates one method from the others. ## ## The second input argument is the order of the stepper. It is needed to ## compute the adaptive timesteps. ## -## The third input argument is a function_handle or an inline function that -## defines the set of ODE: +## The third input argument is a function handle or inline function that +## defines the ODE: ## ## @ifhtml ## @example @@ -56,8 +56,8 @@ ## The fifth input argument represents the fixed timestep while the sixth ## contains the number of integration steps. ## -## The last argument is a struct with the options that may be -## needed by the stepper. +## The last argument is a struct with the options that may be needed by the +## stepper. ## @end deftypefn ## ## @seealso{integrate_adaptive, integrate_const} @@ -98,8 +98,6 @@ solution.vcntloop = 2; solution.vcntcycles = 1; - #vu = vinit; - #vk = vu.' * zeros(1,6); vcntiter = 0; solution.vunhandledtermination = true; solution.vcntsave = 2; @@ -131,15 +129,15 @@ x = [x,u(:,2:end)]; t = [t;z(2:end)]; - solution.vcntsave = solution.vcntsave + 1; - solution.vcntloop = solution.vcntloop + 1; + solution.vcntsave += 1; + solution.vcntloop += 1; vcntiter = 0; - ## Call plot only if a valid result has been found, therefore this code - ## fragment has moved here. Stop integration if plot function returns false + ## Call OutputFcn only if a valid result has been found. + ## Stop integration if function returns false. if (options.vhaveoutputfunction) for vcnt = 0:options.Refine # Approximation between told and t - if (options.vhaverefine) # Do interpolation + if (options.vhaverefine) # Do interpolation vapproxtime = (vcnt + 1) / (options.Refine + 2); vapproxvals = (1 - vapproxtime) * vSaveVUForRefine ... + (vapproxtime) * y(:,end); @@ -163,14 +161,12 @@ endif endif - ## Call event only if a valid result has been found, therefore this - ## code fragment has moved here. Stop integration if veventbreak is - ## true + ## Call Events function only if a valid result has been found. + ## Stop integration if veventbreak is true. if (options.vhaveeventfunction) solution.vevent = odepkg_event_handle (options.Events, t(end), x(:,end), [], options.vfunarguments{:}); - if (! isempty (solution.vevent{1}) - && solution.vevent{1} == 1) + if (! isempty (solution.vevent{1}) && solution.vevent{1} == 1) t(solution.vcntloop-1,:) = solution.vevent{3}(end,:); x(:,solution.vcntloop-1) = solution.vevent{4}(end,:)'; solution.vunhandledtermination = false; @@ -179,10 +175,10 @@ endif ## Update counters that count the number of iteration cycles - solution.vcntcycles = solution.vcntcycles + 1; # Needed for cost statistics - vcntiter = vcntiter + 1; # Needed to find iteration problems + solution.vcntcycles += 1; # Needed for cost statistics + vcntiter += 1; # Needed to find iteration problems - ## Stop solving because the last 1000 steps no successful valid + ## Stop solving because, in the last 5,000 steps, no successful valid ## value has been found if (vcntiter >= 5000) error (["Solving has not been successful. The iterative", @@ -198,21 +194,21 @@ ## Check if integration of the ode has been successful #if (vdirection * z(end) < vdirection * tspan(end)) # if (solution.vunhandledtermination == true) - # error ("OdePkg:InvalidArgument", - # ["Solving has not been successful. The iterative", - # " integration loop exited at time t = %f", - # " before endpoint at tend = %f was reached. This may", - # " happen if the stepsize grows smaller than defined in", - # " vminstepsize. Try to reduce the value of 'InitialStep'", - # " and/or 'MaxStep' with the command 'odeset'.\n"], + # error ("integrate_n_steps:unexpected_termination", + # [" Solving was not successful. ", ... + # " The iterative integration loop exited at time", ... + # " t = %f before the endpoint at tend = %f was reached. ", ... + # " This may happen if the stepsize becomes too small. ", ... + # " Try to reduce the value of 'InitialStep'", ... + # " and/or 'MaxStep' with the command 'odeset'.\n"], # z(end), tspan(end)); # else - # warning ("OdePkg:InvalidArgument", - # ["Solver has been stopped by a call of 'break' in the main", - # " iteration loop at time t = %f before endpoint at tend = %f", - # " was reached. This may happen because the @odeplot function", - # " returned 'true' or the @event function returned", - # " 'true'.\n"], + # warning ("integrate_n_steps:unexpected_termination", + # ["Solver was stopped by a call of 'break'", ... + # " in the main iteration loop at time ", ... + # " t = %f before the endpoint at tend = %f was reached. ", ... + # " This may happen because the @odeplot function ", ... + # " returned 'true' or the @event function returned 'true'.\n"], # z(end), tspan(end)); # endif #endif
--- a/scripts/ode/private/kahan.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/kahan.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,4 +1,4 @@ -## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> ## ## This file is part of Octave. ## @@ -17,37 +17,36 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{sum} =} kahan (@var{sum}, @var{comp}, @var{temp}) -## @deftypefnx {Function File} {[@var{sum}, @var{comp}] =} kahan (@var{sum}, @var{comp}, @var{temp}) +## @deftypefn {Function File} {@var{sum} =} kahan (@var{sum}, @var{comp}, @var{term}) +## @deftypefnx {Function File} {[@var{sum}, @var{comp}] =} kahan (@var{sum}, @var{comp}, @var{term}) ## -## This function is the implementation of the Kahan summation algorithm, -## also known as compensated summation. +## This function implements the Kahan summation algorithm, also known as +## compensated summation. ## -## It significantly reduces the numerical error in the total obtained by adding -## a sequence of finite precision floating point numbers, compared to the -## obvious approach. For more details +## The algorithm significantly reduces the numerical error in the total +## 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 in @command{integrate_adaptive} and in +## This function is called by @command{integrate_adaptive} and ## @command{integrate_const} 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 kahan function. +## next calls to @code{kahan} function. ## ## The second input argument contains the compensation term and it is returned -## as second output argument so that it can be reused in the next calls of the -## same computation. +## as the second output argument so that it can be reused in future calls of +## the same summation. ## -## The third input argument is the variable that contains the term to be added -## to @var{Sum}. +## The third input argument @var{term} is the variable to be added to @var{sum}. ## @end deftypefn -function [Sum, comp] = kahan (Sum, comp, term) +function [sum, comp] = kahan (sum, comp, term) - y = term - comp; - t = Sum + y; - comp = (t - Sum) - y; - Sum = t; + y = term - comp; + t = sum + y; + comp = (t - sum) - y; + sum = t; endfunction
--- a/scripts/ode/private/ode_rk_interpolate.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/ode_rk_interpolate.m Wed Oct 14 10:35:53 2015 -0700 @@ -19,8 +19,9 @@ function u_interp = ode_rk_interpolate (order, z, u, t, k_vals, dt, args) - switch order + switch (order) + ## FIXME: Can interpolations for orders 1-4 simply be deleted? 2015-10-14. #{ case 1 u_interp = linear_interpolation (z, u, t); @@ -40,13 +41,11 @@ 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); + 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 @@ -58,21 +57,21 @@ ## [k_vals(:,1) f_half ... ## k_vals(:,end)], ## tspan(counter)); + otherwise - warning ("High order interpolation not yet implemented: ", - "using cubic interpolation instead"); - der(:,1) = feval (func, z(1) , u(:,1), args); - der(:,2) = feval (func, z(2) , u(:,2), args); + warning (["High order interpolation not yet implemented: ", ... + "using cubic interpolation instead"]); + der(:,1) = feval (func, z(1), u(:,1), args); + der(:,2) = feval (func, z(2), u(:,2), args); u_interp = hermite_cubic_interpolation (z, u, der, t); + endswitch 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. +## 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 = ... @@ -80,9 +79,8 @@ (-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. + ## 4th order approximation of y in t+dt/2 as proposed by Shampine in + ## Lawrence, Shampine, "Some Practical Runge-Kutta Formulas", 1986. dt = t(2) - t(1); u_half = x(:,1) + (1/2) * dt * (der(:,1:7) * coefs_u_half); @@ -97,10 +95,11 @@ ## H4 = s.^2 - 3*s.^3 + 2*s.^4; x_out = zeros (rows (x), length (t_out)); - x_out = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4) .* x(:,1) + ... - ( s - 4*s.^2 + 5*s.^3 - 2*s.^4) .* (dt * der(:,1)) + ... - ( 16*s.^2 - 32*s.^3 + 16*s.^4) .* u_half + ... - ( - 5*s.^2 + 14*s.^3 - 8*s.^4) .* x(:,2) + ... - ( s.^2 - 3*s.^3 + 2*s.^4) .* (dt * der(:,end)); + x_out = (1 - 11*s.^2 + 18*s.^3 - 8*s.^4) .* x(:,1) + ... + ( s - 4*s.^2 + 5*s.^3 - 2*s.^4) .* (dt * der(:,1)) + ... + ( 16*s.^2 - 32*s.^3 + 16*s.^4) .* u_half + ... + ( - 5*s.^2 + 14*s.^3 - 8*s.^4) .* x(:,2) + ... + ( s.^2 - 3*s.^3 + 2*s.^4) .* (dt * der(:,end)); endfunction +
--- a/scripts/ode/private/ode_struct_value_check.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/ode_struct_value_check.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> -## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net> ## ## This file is part of Octave. ## @@ -379,7 +379,7 @@ %! # the command odeset. %! %! ode_struct_value_check (odeset); -%! + %!demo %! # Create the OdePkg options structure A with odeset and check it %! # with odepkg_structure_check. This actually is unnecessary
--- a/scripts/ode/private/odepkg_event_handle.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/odepkg_event_handle.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,4 +1,4 @@ -## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net> ## ## This file is part of Octave. ## @@ -153,5 +153,6 @@ vretcell = cell (1,4); endif + endfunction
--- a/scripts/ode/private/odepkg_structure_check.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/odepkg_structure_check.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> -## Copyright (C) 2006-2012, Thomas Treichl <treichl@users.sourceforge.net> +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2006-2012 Thomas Treichl <treichl@users.sourceforge.net> ## ## This file is part of Octave. ##
--- a/scripts/ode/private/runge_kutta_45_dorpri.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/runge_kutta_45_dorpri.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,5 +1,5 @@ -## Copyright (C) 2015, Carlo de Falco -## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2015 Carlo de Falco +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> ## ## This file is part of Octave. ## @@ -23,7 +23,7 @@ ## @deftypefnx {Function File} {[@var{t_next}, @var{x_next}, @var{x_est}, @var{k_vals_out}] =} runge_kutta_45_dorpri (@var{@@fun}, @var{t}, @var{x}, @var{dt}, @var{options}, @var{k_vals_in}) ## ## 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. +## condition @var{x} from @var{t} to @var{t+dt} with the Dormand-Prince method. ## For the definition of this method see ## @url{http://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method}. ## @@ -36,7 +36,7 @@ ## error. ## ## Fourth output parameter is matrix containing the Runge-Kutta evaluations -## to use in a FSAL scheme or for dense output. +## to use in an FSAL scheme or for dense output. ## ## First input argument is the function describing the system of ODEs to be ## integrated. @@ -52,14 +52,15 @@ ## adapt the computation to what is needed. ## ## Sixth input parameter is optional and describes the Runge-Kutta evaluations -## of the previous step to use in a FSAL scheme. +## of the previous step to use in an FSAL scheme. ## @end deftypefn ## ## @seealso{odepkg} -function [t_out, x_out, x_est, k] = ... - runge_kutta_45_dorpri (f, t, x, dt, opts = [], k_vals = [], - t_out = t + dt) +function [t_next, x_next, x_est, k] = runge_kutta_45_dorpri (f, t, x, dt, + options = [], + k_vals = [], + t_next = t + dt) persistent a = [0 0 0 0 0 0; 1/5 0 0 0 0 0; @@ -70,43 +71,43 @@ 35/384 0 500/1113 125/192 -2187/6784 11/84]; 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)]; - ## x_est according to Shampine 1986: + persistent c_prime = [(5179/57600) 0 (7571/16695) (393/640), ... + (-92097/339200) (187/2100) (1/40)]; + ## According to Shampine 1986: ## persistent c_prime = [(1951/21600) 0 (22642/50085) (451/720), ... ## (-12231/42400) (649/6300) (1/60)]; - persistent c_prime = [(5179/57600) 0 (7571/16695) (393/640), ... - (-92097/339200) (187/2100) (1/40)]; s = t + dt * b; cc = dt * c; aa = dt * a; k = zeros (rows (x), 7); - if (! isempty (opts)) # extra arguments for function evaluator - args = opts.vfunarguments; + if (! isempty (options)) # extra arguments for function evaluator + args = options.vfunarguments; else args = {}; endif - if (! isempty (k_vals)) # k values from previous step are passed + if (! isempty (k_vals)) # k values from previous step are passed k(:,1) = k_vals(:,end); # FSAL property else k(:,1) = feval (f, t, x, args{:}); endif - k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).', args{:}); + k(:,2) = feval (f, s(2), x + k(:,1) * aa(2, 1).' , args{:}); k(:,3) = feval (f, s(3), x + k(:,1:2) * aa(3, 1:2).', args{:}); k(:,4) = feval (f, s(4), x + k(:,1:3) * aa(4, 1:3).', args{:}); k(:,5) = feval (f, s(5), x + k(:,1:4) * aa(5, 1:4).', args{:}); k(:,6) = feval (f, s(6), x + k(:,1:5) * aa(6, 1:5).', args{:}); ## compute new time and new values for the unknowns - ## t_out = t + dt; - x_out = x + k(:,1:6) * cc(:); # 5th order approximation + ## t_next = t + dt; + x_next = x + k(:,1:6) * cc(:); # 5th order approximation ## if the estimation of the error is required if (nargout >= 3) ## new solution to be compared with the previous one - k(:,7) = feval (f, t_out, x_out, args{:}); + k(:,7) = feval (f, t_next, x_next, args{:}); cc_prime = dt * c_prime; x_est = x + k * cc_prime(:); endif
--- a/scripts/ode/private/starting_stepsize.m Wed Oct 14 09:25:04 2015 -0700 +++ b/scripts/ode/private/starting_stepsize.m Wed Oct 14 10:35:53 2015 -0700 @@ -1,4 +1,4 @@ -## Copyright (C) 2013, Roberto Porcu' <roberto.porcu@polimi.it> +## Copyright (C) 2013 Roberto Porcu' <roberto.porcu@polimi.it> ## ## This file is part of Octave. ## @@ -17,20 +17,20 @@ ## <http://www.gnu.org/licenses/>. ## -*- texinfo -*- -## @deftypefn {Function File} {@var{h} =} starting_stepsize (@var{order}, @var{@@fun}, @var{t0}, @var{x0}) +## @deftypefn {Function File} {@var{h} =} starting_stepsize (@var{order}, @var{@@func}, @var{t0}, @var{x0}, @var{AbsTol}, @var{RelTol}, @var{normcontrol}) ## -## This function file can be used to determine a good initial step for an ODE -## solver of order @var{order}. The algorithm is that one described in [1]. +## Determine a good initial timestep for an ODE solver of order @var{order} +## using the algorithm described in reference [1]. ## -## Second input argument, which is @var{@@fun}, is the function describing -## the differential equations, @var{t0} is the initial time and @var{x0} is -## the initial condition. +## 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. ## -## This function returns a good guess for the initial timestep @var{h}. -## ## References: ## [1] E. Hairer, S.P. Norsett and G. Wanner, -## "Solving Ordinary Differential Equations I: Nonstiff Problems", Springer. +## @cite{Solving Ordinary Differential Equations I: Nonstiff Problems}, +## Springer. ## @end deftypefn ## ## @seealso{odepkg} @@ -59,7 +59,7 @@ AbsRel_Norm (func (t0+h0, x1) - y, func (t0+h0, x1) - y, AbsTol, RelTol, normcontrol); - if (max(d1, d2) <= 1e-15) + if (max (d1, d2) <= 1e-15) h1 = max (1e-6, h0*1e-3); else h1 = (1e-2 / max (d1, d2)) ^(1 / (order+1));