| 30 | 30 | // Step size control: |
|
| 31 | 31 | SEXP r_step_size_min, SEXP r_step_size_max, |
|
| 32 | 32 | SEXP r_step_size_initial, SEXP r_step_max_n, |
|
| 33 | + | SEXP r_step_size_min_allow, |
|
| 33 | 34 | // Critical times and events: |
|
| 34 | 35 | SEXP r_tcrit, SEXP r_is_event, SEXP r_events, |
|
| 35 | 36 | // Other: |
| 159 | 160 | obj->step_size_max = fmin(fabs(REAL(r_step_size_max)[0]), DBL_MAX); |
|
| 160 | 161 | obj->step_size_initial = REAL(r_step_size_initial)[0]; |
|
| 161 | 162 | obj->step_max_n = INTEGER(r_step_max_n)[0]; |
|
| 163 | + | obj->step_size_min_allow = INTEGER(r_step_size_min_allow)[0]; |
|
| 162 | 164 | ||
| 163 | 165 | obj->stiff_check = INTEGER(r_stiff_check)[0]; |
|
| 164 | 166 |
| 51 | 51 | ##' used will be the largest of the absolute value of this |
|
| 52 | 52 | ##' \code{step_size_min} or \code{.Machine$double.eps}. If the |
|
| 53 | 53 | ##' integration attempts to make a step smaller than this, it will |
|
| 54 | - | ##' throw an error, stopping the integration (note that this differs |
|
| 55 | - | ##' from the treatment of \code{hmin} in \code{deSolve::lsoda}). |
|
| 54 | + | ##' throw an error by default, stopping the integration (note that |
|
| 55 | + | ##' this differs from the treatment of \code{hmin} in |
|
| 56 | + | ##' \code{deSolve::lsoda}). See \code{allow_step_size_min} to change |
|
| 57 | + | ##' this behaviour. |
|
| 56 | 58 | ##' |
|
| 57 | 59 | ##' @param step_size_max The largest step size. By default there is |
|
| 58 | 60 | ##' no maximum step size (Inf) so the solver can take as large a |
| 70 | 72 | ##' the number of steps (or 11 times if using \code{method = |
|
| 71 | 73 | ##' "dopri853"}). |
|
| 72 | 74 | ##' |
|
| 75 | + | ##' @param step_size_min_allow Logical, indicating if when a step size |
|
| 76 | + | ##' is driven down to \code{step_size_min} we should allow it to |
|
| 77 | + | ##' proceed. This is the behaviour in of \code{hmin} in |
|
| 78 | + | ##' \code{deSolve::lsoda}. |
|
| 79 | + | ##' |
|
| 73 | 80 | ##' @param tcrit An optional vector of critical times that the solver |
|
| 74 | 81 | ##' must stop at (rather than interpolating over). This can include |
|
| 75 | 82 | ##' an end time that we can't go past, or points within the |
| 288 | 295 | rtol = 1e-6, atol = 1e-6, |
|
| 289 | 296 | step_size_min = 0, step_size_max = Inf, |
|
| 290 | 297 | step_size_initial = 0, step_max_n = 100000L, |
|
| 298 | + | step_size_min_allow = FALSE, |
|
| 291 | 299 | tcrit = NULL, event_time = NULL, event_function = NULL, |
|
| 292 | 300 | method = "dopri5", |
|
| 293 | 301 | stiff_check = 0, |
| 335 | 343 | assert_scalar(step_size_max) |
|
| 336 | 344 | assert_scalar(step_size_initial) |
|
| 337 | 345 | assert_size(step_max_n) |
|
| 346 | + | assert_scalar_logical(step_size_min_allow) |
|
| 338 | 347 | assert_size(n_history) |
|
| 339 | 348 | ||
| 340 | 349 | assert_scalar_logical(grow_history) |
| 398 | 407 | ## Step control: |
|
| 399 | 408 | step_size_min, step_size_max, |
|
| 400 | 409 | step_size_initial, as.integer(step_max_n), |
|
| 410 | + | step_size_min_allow, |
|
| 401 | 411 | ## Critical and events |
|
| 402 | 412 | as.numeric(tcrit), is_event, event, |
|
| 403 | 413 | ## Other: |
| 78 | 78 | ret->step_size_min = DBL_EPSILON; |
|
| 79 | 79 | ret->step_size_max = DBL_MAX; |
|
| 80 | 80 | ret->step_size_initial = 0.0; |
|
| 81 | + | ret->step_size_min_allow = false; |
|
| 81 | 82 | ret->step_max_n = 100000; // from dopri5.f:212 |
|
| 82 | 83 | ret->step_factor_safe = 0.9; // from dopri5.f:265 |
|
| 83 | 84 |
| 110 | 111 | ret->step_size_max = obj->step_size_max; |
|
| 111 | 112 | ret->step_size_initial = obj->step_size_initial; |
|
| 112 | 113 | ret->step_max_n = obj->step_max_n; |
|
| 114 | + | ret->step_size_min_allow = obj->step_size_min_allow; |
|
| 113 | 115 | ret->step_beta = obj->step_beta; |
|
| 114 | 116 | ret->step_constant = obj->step_constant; |
|
| 115 | 117 |
| 338 | 340 | ||
| 339 | 341 | // TODO: factor into its own thing |
|
| 340 | 342 | while (true) { |
|
| 343 | + | bool force_this_step = false; |
|
| 341 | 344 | if (obj->n_step > obj->step_max_n) { |
|
| 342 | 345 | obj->error = true; |
|
| 343 | 346 | obj->code = ERR_TOO_MANY_STEPS; |
|
| 344 | 347 | break; |
|
| 345 | 348 | } |
|
| 346 | 349 | if (fabs(h) <= obj->step_size_min) { |
|
| 347 | - | obj->error = true; |
|
| 348 | - | obj->code = ERR_STEP_SIZE_TOO_SMALL; |
|
| 349 | - | break; |
|
| 350 | + | if (obj->step_size_min_allow) { |
|
| 351 | + | h = copysign(obj->step_size_min, obj->sign); |
|
| 352 | + | force_this_step = true; |
|
| 353 | + | } else { |
|
| 354 | + | obj->error = true; |
|
| 355 | + | obj->code = ERR_STEP_SIZE_TOO_SMALL; |
|
| 356 | + | break; |
|
| 357 | + | } |
|
| 350 | 358 | } else if (fabs(h) <= fabs(obj->t) * DBL_EPSILON) { |
|
| 351 | 359 | obj->error = true; |
|
| 352 | 360 | obj->code = ERR_STEP_SIZE_VANISHED; |
| 384 | 392 | double h_new = dopri_h_new(obj, fac_old, h, err); |
|
| 385 | 393 | const bool accept = err <= 1; |
|
| 386 | 394 | ||
| 387 | - | if (accept) { |
|
| 395 | + | if (accept || force_this_step) { |
|
| 388 | 396 | // Step is accepted :) |
|
| 389 | 397 | fac_old = fmax(err, 1e-4); |
|
| 390 | 398 | obj->n_accept++; |
| 8 | 8 | #include <Rversion.h> |
|
| 9 | 9 | ||
| 10 | 10 | static const R_CallMethodDef call_methods[] = { |
|
| 11 | - | {"Cdopri", (DL_FUNC) &r_dopri, 26}, |
|
| 11 | + | {"Cdopri", (DL_FUNC) &r_dopri, 27}, |
|
| 12 | 12 | {"Cdopri_continue", (DL_FUNC) &r_dopri_continue, 10}, |
|
| 13 | 13 | {"Cdopri_copy", (DL_FUNC) &r_dopri_copy, 1}, |
|
| 14 | 14 | {"Cylag", (DL_FUNC) &r_ylag, 2}, |
richfitz
