1######################################################################## 2## 3## Copyright (C) 2013-2021 The Octave Project Developers 4## 5## See the file COPYRIGHT.md in the top-level directory of this 6## distribution or <https://octave.org/copyright/>. 7## 8## This file is part of Octave. 9## 10## Octave is free software: you can redistribute it and/or modify it 11## under the terms of the GNU General Public License as published by 12## the Free Software Foundation, either version 3 of the License, or 13## (at your option) any later version. 14## 15## Octave is distributed in the hope that it will be useful, but 16## WITHOUT ANY WARRANTY; without even the implied warranty of 17## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 18## GNU General Public License for more details. 19## 20## You should have received a copy of the GNU General Public License 21## along with Octave; see the file COPYING. If not, see 22## <https://www.gnu.org/licenses/>. 23## 24######################################################################## 25 26## -*- texinfo -*- 27## @deftypefn {} {@var{solution} =} integrate_adaptive (@var{@@stepper}, @var{order}, @var{@@func}, @var{tspan}, @var{x0}, @var{options}) 28## 29## This function file can be called by an ODE solver function in order to 30## integrate the set of ODEs on the interval @var{[t0, t1]} with an adaptive 31## timestep. 32## 33## The function returns a structure @var{solution} with two fields: @var{t} 34## and @var{y}. @var{t} is a column vector and contains the time stamps. 35## @var{y} is a matrix in which each column refers to a different unknown 36## of the problem and the row number is the same as the @var{t} row number. 37## Thus, each row of the matrix @var{y} contains the values of all unknowns at 38## the time value contained in the corresponding row in @var{t}. 39## 40## The first input argument must be a function handle or inline function 41## representing the stepper, i.e., the function responsible for step-by-step 42## integration. This function discriminates one method from the others. 43## 44## The second input argument is the order of the stepper. It is needed 45## to compute the adaptive timesteps. 46## 47## The third input argument is a function handle or inline function that 48## defines the ODE: 49## 50## @ifhtml 51## 52## @example 53## @math{y' = f(t,y)} 54## @end example 55## 56## @end ifhtml 57## @ifnothtml 58## @math{y' = f(t,y)}. 59## @end ifnothtml 60## 61## The fourth input argument is the time vector which defines the integration 62## interval, i.e., @var{[tspan(1), tspan(end)]} and all intermediate elements 63## are taken as times at which the solution is required. 64## 65## The fifth argument represents the initial conditions for the ODEs and the 66## last input argument contains some options that may be needed for the 67## stepper. 68## 69## @end deftypefn 70 71function solution = integrate_adaptive (stepper, order, func, tspan, x0, 72 options) 73 74 fixed_times = numel (tspan) > 2; 75 76 t_new = t_old = t = tspan(1); 77 x_new = x_old = x = x0(:); 78 79 ## Get first initial timestep 80 dt = options.InitialStep; 81 if (isempty (dt)) 82 dt = starting_stepsize (order, func, t, x, 83 options.AbsTol, options.RelTol, 84 strcmp (options.NormControl, "on"), 85 options.funarguments); 86 endif 87 88 dir = options.direction; 89 dt = dir * min (abs (dt), options.MaxStep); 90 91 options.comp = 0.0; 92 93 ## Factor multiplying the stepsize guess 94 facmin = 0.8; 95 facmax = 1.5; 96 fac = 0.38^(1/(order+1)); # formula taken from Hairer 97 98 ## Initialize the OutputFcn 99 if (options.haveoutputfunction) 100 if (! isempty (options.OutputSel)) 101 solution.retout = x(options.OutputSel,end); 102 else 103 solution.retout = x; 104 endif 105 feval (options.OutputFcn, tspan, solution.retout, "init", 106 options.funarguments{:}); 107 endif 108 109 ## Initialize the EventFcn 110 have_EventFcn = false; 111 if (! isempty (options.Events)) 112 have_EventFcn = true; 113 ode_event_handler (options.Events, tspan(1), x, 114 "init", options.funarguments{:}); 115 endif 116 117 if (options.havenonnegative) 118 nn = options.NonNegative; 119 endif 120 121 solution.cntloop = 0; 122 solution.cntcycles = 0; 123 solution.cntsave = 2; 124 solution.unhandledtermination = true; 125 ireject = 0; 126 127 NormControl = strcmp (options.NormControl, "on"); 128 k_vals = []; 129 iout = istep = 1; 130 131 while (dir * t_old < dir * tspan(end)) 132 133 ## Compute integration step from t_old to t_new = t_old + dt 134 [t_new, options.comp] = kahan (t_old, options.comp, dt); 135 [t_new, x_new, x_est, new_k_vals] = ... 136 stepper (func, t_old, x_old, dt, options, k_vals, t_new); 137 138 solution.cntcycles += 1; 139 140 if (options.havenonnegative) 141 x_new(nn, end) = abs (x_new(nn, end)); 142 x_est(nn, end) = abs (x_est(nn, end)); 143 endif 144 145 err = AbsRel_norm (x_new, x_old, options.AbsTol, options.RelTol, 146 NormControl, x_est); 147 148 ## Accept solution only if err <= 1.0 149 if (err <= 1) 150 151 solution.cntloop += 1; 152 ireject = 0; # Clear reject counter 153 154 ## if output time steps are fixed 155 if (fixed_times) 156 157 t_caught = find ((dir * tspan(iout:end) > dir * t_old) 158 & (dir * tspan(iout:end) <= dir * t_new)); 159 t_caught = t_caught + iout - 1; 160 161 if (! isempty (t_caught)) 162 t(t_caught) = tspan(t_caught); 163 iout = max (t_caught); 164 x(:, t_caught) = ... 165 runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], 166 t(t_caught), new_k_vals, dt, func, 167 options.funarguments); 168 169 istep += 1; 170 171 ## Call Events function only if a valid result has been found. 172 ## Stop integration if eventbreak is true. 173 if (have_EventFcn) 174 break_loop = false; 175 for idenseout = 1:numel (t_caught) 176 id = t_caught(idenseout); 177 td = t(id); 178 solution.event = ... 179 ode_event_handler (options.Events, t(id), x(:, id), [], 180 options.funarguments{:}); 181 if (! isempty (solution.event{1}) && solution.event{1} == 1) 182 t(id) = solution.event{3}(end); 183 t = t(1:id); 184 x(:, id) = solution.event{4}(end, :).'; 185 x = x(:,1:id); 186 solution.unhandledtermination = false; 187 break_loop = true; 188 break; 189 endif 190 endfor 191 if (break_loop) 192 break; 193 endif 194 endif 195 196 ## Call OutputFcn only if a valid result has been found. 197 ## Stop integration if function returns true. 198 if (options.haveoutputfunction) 199 cnt = options.Refine + 1; 200 approxtime = linspace (t_old, t_new, cnt); 201 approxvals = interp1 ([t_old, t(t_caught), t_new], 202 [x_old, x(:, t_caught), x_new] .', 203 approxtime, "linear") .'; 204 if (isvector (approxvals) && ! isscalar (approxtime)) 205 approxvals = approxvals.'; 206 endif 207 if (! isempty (options.OutputSel)) 208 approxvals = approxvals(options.OutputSel, :); 209 endif 210 stop_solve = false; 211 for ii = 1:numel (approxtime) 212 stop_solve = feval (options.OutputFcn, 213 approxtime(ii), approxvals(:, ii), [], 214 options.funarguments{:}); 215 if (stop_solve) 216 break; # break from inner loop 217 endif 218 endfor 219 if (stop_solve) # Leave main loop 220 solution.unhandledtermination = false; 221 break; 222 endif 223 endif 224 225 endif 226 227 else # not fixed times 228 229 t(++istep) = t_new; 230 x(:, istep) = x_new; 231 iout = istep; 232 233 ## Call Events function only if a valid result has been found. 234 ## Stop integration if eventbreak is true. 235 if (have_EventFcn) 236 solution.event = ... 237 ode_event_handler (options.Events, t(istep), x(:, istep), [], 238 options.funarguments{:}); 239 if (! isempty (solution.event{1}) && solution.event{1} == 1) 240 t(istep) = solution.event{3}(end); 241 x(:, istep) = solution.event{4}(end, :).'; 242 solution.unhandledtermination = false; 243 break; 244 endif 245 endif 246 247 ## Call OutputFcn only if a valid result has been found. 248 ## Stop integration if function returns true. 249 if (options.haveoutputfunction) 250 cnt = options.Refine + 1; 251 approxtime = linspace (t_old, t_new, cnt); 252 approxvals = interp1 ([t_old, t_new], 253 [x_old, x_new] .', 254 approxtime, "linear") .'; 255 if (isvector (approxvals) && ! isscalar (approxtime)) 256 approxvals = approxvals.'; 257 endif 258 if (! isempty (options.OutputSel)) 259 approxvals = approxvals(options.OutputSel, :); 260 endif 261 stop_solve = false; 262 for ii = 1:numel (approxtime) 263 stop_solve = feval (options.OutputFcn, 264 approxtime(ii), approxvals(:, ii), [], 265 options.funarguments{:}); 266 if (stop_solve) 267 break; # break from inner loop 268 endif 269 endfor 270 if (stop_solve) # Leave main loop 271 solution.unhandledtermination = false; 272 break; 273 endif 274 endif 275 276 endif 277 278 ## move to next time-step 279 t_old = t_new; 280 x_old = x_new; 281 k_vals = new_k_vals; 282 283 else # error condition 284 285 ireject += 1; 286 287 ## Stop solving if, in the last 5,000 steps, no successful valid 288 ## value has been found. 289 if (ireject >= 5_000) 290 error (["integrate_adaptive: Solving was not successful. ", ... 291 " The iterative integration loop exited at time", ... 292 " t = %f before the endpoint at tend = %f was reached. ", ... 293 " This happened because the iterative integration loop", ... 294 " did not find a valid solution at this time stamp. ", ... 295 " Try to reduce the value of 'InitialStep' and/or", ... 296 " 'MaxStep' with the command 'odeset'.\n"], 297 t_old, tspan(end)); 298 endif 299 300 endif 301 302 ## Compute next timestep, formula taken from Hairer 303 err += eps; # avoid divisions by zero 304 dt *= min (facmax, max (facmin, fac * (1 / err)^(1 / (order + 1)))); 305 dt = dir * min (abs (dt), options.MaxStep); 306 if (! (abs (dt) > eps (t(end)))) 307 break; 308 endif 309 310 ## Make sure we don't go past tpan(end) 311 dt = dir * min (abs (dt), abs (tspan(end) - t_old)); 312 313 endwhile 314 315 ## Check if integration of the ode has been successful 316 if (dir * t(end) < dir * tspan(end)) 317 if (solution.unhandledtermination == true) 318 warning ("integrate_adaptive:unexpected_termination", 319 [" Solving was not successful. ", ... 320 " The iterative integration loop exited at time", ... 321 " t = %f before the endpoint at tend = %f was reached. ", ... 322 " This may happen if the stepsize becomes too small. ", ... 323 " Try to reduce the value of 'InitialStep'", ... 324 " and/or 'MaxStep' with the command 'odeset'."], 325 t(end), tspan(end)); 326 else 327 warning ("integrate_adaptive:unexpected_termination", 328 ["Solver was stopped by a call of 'break'", ... 329 " in the main iteration loop at time", ... 330 " t = %f before the endpoint at tend = %f was reached. ", ... 331 " This may happen because the @odeplot function", ... 332 " returned 'true' or the @event function returned 'true'."], 333 t(end), tspan(end)); 334 endif 335 endif 336 337 ## Set up return structure 338 solution.t = t(:); 339 solution.x = x.'; 340 341endfunction 342