Mercurial > octave
view scripts/ode/private/integrate_adaptive.m @ 30564:796f54d4ddbf stable
update Octave Project Developers copyright for the new year
In files that have the "Octave Project Developers" copyright notice,
update for 2021.
In all .txi and .texi files except gpl.txi and gpl.texi in the
doc/liboctave and doc/interpreter directories, change the copyright
to "Octave Project Developers", the same as used for other source
files. Update copyright notices for 2022 (not done since 2019). For
gpl.txi and gpl.texi, change the copyright notice to be "Free Software
Foundation, Inc." and leave the date at 2007 only because this file
only contains the text of the GPL, not anything created by the Octave
Project Developers.
Add Paul Thomas to contributors.in.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 28 Dec 2021 18:22:40 -0500 |
parents | 0a5b15007766 |
children | e1788b1a315f |
line wrap: on
line source
######################################################################## ## ## Copyright (C) 2013-2022 The Octave Project Developers ## ## See the file COPYRIGHT.md in the top-level directory of this ## distribution or <https://octave.org/copyright/>. ## ## 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 ## <https://www.gnu.org/licenses/>. ## ######################################################################## ## -*- texinfo -*- ## @deftypefn {} {@var{solution} =} integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@func}, @var{tspan}, @var{x0}, @var{options}) ## ## This function file can be called by an ODE solver function in order to ## integrate the set of ODEs on the interval @var{[t0, t1]} with an adaptive ## timestep. ## ## The function returns a structure @var{solution} with two fields: @var{t} ## and @var{y}. @var{t} is a column vector and contains the time stamps. ## @var{y} is a matrix in which each column refers to a different unknown ## of the problem and the row number is the same as the @var{t} row number. ## 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 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 inline function that ## defines the ODE: ## ## @ifhtml ## ## @example ## @math{y' = f(t,y)} ## @end example ## ## @end ifhtml ## @ifnothtml ## @math{y' = f(t,y)}. ## @end ifnothtml ## ## 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. ## ## @end deftypefn 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 = options.InitialStep; if (isempty (dt)) dt = starting_stepsize (order, func, t, x, options.AbsTol, options.RelTol, strcmp (options.NormControl, "on"), options.funarguments); endif dir = options.direction; dt = dir * min (abs (dt), options.MaxStep); options.comp = 0.0; ## Factor multiplying the stepsize guess facmin = 0.8; facmax = 1.5; fac = 0.38^(1/(order+1)); # formula taken from Hairer ## Initialize the OutputFcn if (options.haveoutputfunction) if (! isempty (options.OutputSel)) solution.retout = x(options.OutputSel,end); else solution.retout = x; endif feval (options.OutputFcn, tspan, solution.retout, "init", options.funarguments{:}); 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 if (options.havenonnegative) nn = options.NonNegative; endif solution.cntloop = 0; solution.cntcycles = 0; solution.cntsave = 2; solution.unhandledtermination = true; ireject = 0; NormControl = strcmp (options.NormControl, "on"); k_vals = []; iout = istep = 1; while (dir * t_old < dir * tspan(end)) ## Compute integration step from t_old to t_new = t_old + dt [t_new, options.comp] = kahan (t_old, options.comp, dt); [t_new, x_new, x_est, new_k_vals] = ... stepper (func, t_old, x_old, dt, options, k_vals, t_new); solution.cntcycles += 1; if (options.havenonnegative) x_new(nn, end) = abs (x_new(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) solution.cntloop += 1; ireject = 0; # Clear reject counter ## if output time steps are fixed if (fixed_times) t_caught = find ((dir * tspan(iout:end) > dir * t_old) & (dir * tspan(iout:end) <= dir * t_new)); t_caught = t_caught + iout - 1; if (! isempty (t_caught)) t(t_caught) = tspan(t_caught); iout = max (t_caught); x(:, t_caught) = ... runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], t(t_caught), new_k_vals, dt, func, options.funarguments); istep += 1; ## Call Events function only if a valid result has been found. ## Stop integration if eventbreak is true. if (have_EventFcn) break_loop = false; for idenseout = 1:numel (t_caught) id = t_caught(idenseout); td = t(id); solution.event = ... ode_event_handler (options.Events, t(id), x(:, id), [], options.funarguments{:}); if (! isempty (solution.event{1}) && solution.event{1} == 1) t(id) = solution.event{3}(end); t = t(1:id); x(:, id) = solution.event{4}(end, :).'; x = x(:,1:id); solution.unhandledtermination = false; break_loop = true; break; endif endfor if (break_loop) break; endif endif ## Call OutputFcn only if a valid result has been found. ## Stop integration if function returns true. if (options.haveoutputfunction) cnt = options.Refine + 1; approxtime = linspace (t_old, t_new, cnt); approxvals = interp1 ([t_old, t(t_caught), t_new], [x_old, x(:, t_caught), x_new] .', approxtime, "linear") .'; if (isvector (approxvals) && ! isscalar (approxtime)) approxvals = approxvals.'; endif if (! isempty (options.OutputSel)) approxvals = approxvals(options.OutputSel, :); endif stop_solve = false; for ii = 1:numel (approxtime) stop_solve = feval (options.OutputFcn, approxtime(ii), approxvals(:, ii), [], options.funarguments{:}); if (stop_solve) break; # break from inner loop endif endfor if (stop_solve) # Leave main loop solution.unhandledtermination = false; break; endif endif endif else # not fixed times t(++istep) = t_new; x(:, istep) = x_new; iout = istep; ## Call Events function only if a valid result has been found. ## Stop integration if eventbreak is true. if (have_EventFcn) solution.event = ... ode_event_handler (options.Events, t(istep), x(:, istep), [], options.funarguments{:}); if (! isempty (solution.event{1}) && solution.event{1} == 1) t(istep) = solution.event{3}(end); x(:, istep) = solution.event{4}(end, :).'; solution.unhandledtermination = false; break; endif endif ## Call OutputFcn only if a valid result has been found. ## Stop integration if function returns true. if (options.haveoutputfunction) cnt = options.Refine + 1; approxtime = linspace (t_old, t_new, cnt); approxvals = interp1 ([t_old, t_new], [x_old, x_new] .', approxtime, "linear") .'; if (isvector (approxvals) && ! isscalar (approxtime)) approxvals = approxvals.'; endif if (! isempty (options.OutputSel)) approxvals = approxvals(options.OutputSel, :); endif stop_solve = false; for ii = 1:numel (approxtime) stop_solve = feval (options.OutputFcn, approxtime(ii), approxvals(:, ii), [], options.funarguments{:}); if (stop_solve) break; # break from inner loop endif endfor if (stop_solve) # Leave main loop solution.unhandledtermination = false; break; endif endif endif ## move to next time-step t_old = t_new; x_old = x_new; k_vals = new_k_vals; else # error condition ireject += 1; ## Stop solving if, in the last 5,000 steps, no successful valid ## value has been found. if (ireject >= 5_000) error (["integrate_adaptive: Solving was not successful. ", ... " The iterative integration loop exited at time", ... " 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 endif ## Compute next timestep, formula taken from Hairer err += eps; # avoid divisions by zero dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1)))); dt = dir * min (abs (dt), options.MaxStep); if (! (abs (dt) > eps (t(end)))) break; endif ## Make sure we don't go past tpan(end) dt = dir * min (abs (dt), abs (tspan(end) - t_old)); endwhile ## Check if integration of the ode has been successful if (dir * t(end) < dir * tspan(end)) if (solution.unhandledtermination == true) warning ("integrate_adaptive: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'."], t(end), tspan(end)); else warning ("integrate_adaptive: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'."], t(end), tspan(end)); endif endif ## Set up return structure solution.t = t(:); solution.x = x.'; endfunction