diff --git a/docs/changelog.rst b/docs/changelog.rst index 34dac7d..a1685cc 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -13,8 +13,6 @@ Features ^^^^^^^^ * Added an example with the fractional Lorenz system (:ghpr:`13`). -* Support setting a constant for - :meth:`~pycaputo.fode.FractionalDifferentialEquationMethod.predict_time_step`. * Add a guess for the number of corrector iterations for :class:`~pycaputo.fode.CaputoPECEMethod` from [Garrappa2010]_. * Implement :class:`~pycaputo.quadrature.RiemannLiouvilleSimpsonMethod`, a diff --git a/docs/fode.rst b/docs/fode.rst index 7ad9004..88b60c4 100644 --- a/docs/fode.rst +++ b/docs/fode.rst @@ -17,13 +17,27 @@ History Handling .. autoclass:: FixedState .. autoclass:: FixedSizeHistory -Interface ---------- +Time Interval Handling +---------------------- + +.. exception:: StepEstimateError + +.. autoclass:: TimeSpan +.. autoclass:: FixedTimeSpan +.. autoclass:: GradedTimeSpan +.. autoclass:: FixedLipschitzTimeSpan +.. autoclass:: LipschitzTimeSpan + +Time Stepping Events +-------------------- .. autoclass:: Event .. autoclass:: StepFailed .. autoclass:: StepCompleted +Time Stepping Interface +----------------------- + .. autoclass:: FractionalDifferentialEquationMethod .. autoclass:: ProductIntegrationMethod @@ -31,9 +45,6 @@ Interface .. autofunction:: advance .. autofunction:: make_initial_condition -.. autofunction:: make_predict_time_step_fixed -.. autofunction:: make_predict_time_step_graded - Caputo Derivative FODEs ======================= diff --git a/docs/references.rst b/docs/references.rst index 33eb39c..706e467 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -31,6 +31,10 @@ References International Journal of Computer Mathematics, Vol. 87, pp. 2281--2290, 2010, `DOI `__. +.. [Baleanu2012] D. Baleanu, K. Diethelm, E. Scalas, J. J. Trujillo, + *Fractional Calculus - Models and Numerical Methods*, + World Scientific, 2012. + .. [Garrappa2015] R. Garrappa, *Numerical Evaluation of Two and Three Parameter Mittag-Leffler Functions*, SIAM Journal on Numerical Analysis, Vol. 53, pp. 1350--1369, 2015, diff --git a/examples/brusselator-predictor-corrector.py b/examples/brusselator-predictor-corrector.py index 197637d..0cb65e4 100644 --- a/examples/brusselator-predictor-corrector.py +++ b/examples/brusselator-predictor-corrector.py @@ -24,7 +24,7 @@ def brusselator(t: float, y: Array, *, a: float, mu: float) -> Array: # {{{ solve -from pycaputo.fode import CaputoPECEMethod +from pycaputo.fode import CaputoPECEMethod, FixedTimeSpan alpha = 0.8 a = 1.0 @@ -33,9 +33,8 @@ def brusselator(t: float, y: Array, *, a: float, mu: float) -> Array: stepper = CaputoPECEMethod( derivative_order=(alpha,), - predict_time_step=1.0e-2, + tspan=FixedTimeSpan.from_data(1.0e-2, tstart=0.0, tfinal=50.0), source=partial(brusselator, a=a, mu=mu), - tspan=(0, 50), y0=(y0,), corrector_iterations=1, ) diff --git a/examples/fractional-lorenz.py b/examples/fractional-lorenz.py index 3b1ea99..28cb713 100644 --- a/examples/fractional-lorenz.py +++ b/examples/fractional-lorenz.py @@ -36,7 +36,7 @@ def lorenz_jac(t: float, y: Array, *, sigma: float, rho: float, beta: float) -> # {{{ solve -from pycaputo.fode import CaputoWeightedEulerMethod +from pycaputo.fode import CaputoWeightedEulerMethod, FixedTimeSpan # NOTE: order example taken from https://doi.org/10.1016/j.chaos.2009.03.016 alpha = (0.985, 0.99, 0.99) @@ -49,10 +49,9 @@ def lorenz_jac(t: float, y: Array, *, sigma: float, rho: float, beta: float) -> stepper = CaputoWeightedEulerMethod( derivative_order=alpha, - predict_time_step=1.0e-2, + tspan=FixedTimeSpan.from_data(1.0e-2, tstart=0.0, tfinal=75.0), source=partial(lorenz, sigma=sigma, rho=rho, beta=beta), source_jac=partial(lorenz_jac, sigma=sigma, rho=rho, beta=beta), - tspan=(0, 75), y0=(y0,), theta=1.0, ) diff --git a/pycaputo/fode/__init__.py b/pycaputo/fode/__init__.py index c7492c7..acc90c4 100644 --- a/pycaputo/fode/__init__.py +++ b/pycaputo/fode/__init__.py @@ -3,14 +3,18 @@ from pycaputo.fode.base import ( Event, + FixedLipschitzTimeSpan, + FixedTimeSpan, FractionalDifferentialEquationMethod, + GradedTimeSpan, + LipschitzTimeSpan, StepCompleted, + StepEstimateError, StepFailed, + TimeSpan, advance, evolve, make_initial_condition, - make_predict_time_step_fixed, - make_predict_time_step_graded, ) from pycaputo.fode.caputo import ( CaputoForwardEulerMethod, @@ -47,6 +51,12 @@ "FixedState", "FractionalDifferentialEquationMethod", "History", + "StepEstimateError", + "TimeSpan", + "FixedTimeSpan", + "FixedLipschitzTimeSpan", + "GradedTimeSpan", + "LipschitzTimeSpan", "ProductIntegrationMethod", "ProductIntegrationState", "State", @@ -56,6 +66,4 @@ "advance", "evolve", "make_initial_condition", - "make_predict_time_step_fixed", - "make_predict_time_step_graded", ) diff --git a/pycaputo/fode/base.py b/pycaputo/fode/base.py index 7b74b1f..864d889 100644 --- a/pycaputo/fode/base.py +++ b/pycaputo/fode/base.py @@ -5,13 +5,15 @@ from abc import ABC, abstractmethod from dataclasses import dataclass -from functools import cached_property, singledispatch +from functools import cached_property, partial, singledispatch from typing import Iterator +import numpy as np + from pycaputo.derivatives import FractionalOperator from pycaputo.fode.history import History from pycaputo.logging import get_logger -from pycaputo.utils import Array, CallbackFunction, ScalarStateFunction, StateFunction +from pycaputo.utils import Array, CallbackFunction, StateFunction logger = get_logger(__name__) @@ -19,119 +21,244 @@ # {{{ time step -def make_predict_time_step_fixed(dt: float) -> ScalarStateFunction: - """ - :returns: a callable returning a fixed time step *dt*. - """ +class StepEstimateError(RuntimeError): + """An exception raised when a time step estimate has failed.""" + + +@dataclass(frozen=True) +class TimeSpan(ABC): + """A description of a discrete time interval.""" + + #: Start of the time interval. + tstart: float + #: End of the time interval (leave *None* for infinite time stepping). + tfinal: float | None + #: Number of time steps (leave *None* for infinite time stepping). + nsteps: int | None - if dt < 0: - raise ValueError(f"Time step should be positive: {dt}") + if __debug__: + + def __post_init__(self) -> None: + if self.tfinal is not None and self.tstart > self.tfinal: + raise ValueError("Invalid time interval: 'tstart' > 'tfinal'") + + if self.nsteps is not None and self.nsteps <= 0: + raise ValueError( + f"Number of iterations must be positive: {self.nsteps}" + ) + + def finished(self, n: int, t: float) -> bool: + """Check if the evolution should finish at iteration *n* and time *t*.""" + if self.tfinal is not None and t >= self.tfinal: + return True + + if self.nsteps is not None and n >= self.nsteps: + return True + + return False + + @abstractmethod + def get_next_time_step_raw(self, n: int, t: float, y: Array) -> float: + """Get a raw estimate of the time step at the given state values. + + This function is meant to be implemented by subclasses, but + :meth:`~pycaputo.fode.TimeSpan.get_next_time_step` should be used to + ensure a consistent time stepping. + """ + + def get_next_time_step(self, n: int, t: float, y: Array) -> float: + r"""Get an estimate of the time step at the given state values. + + :arg n: current iteration. + :arg t: starting time, i.e. the time step is estimated for the interval + :math:`[t, t + \Delta t]`. + :arg y: state value at the time *t*. + :returns: an estimate for the time step :math:`\Delta t`. + :raises StepEstimateError: When the next time step cannot be estimated + or is not finite. + """ + try: + dt = self.get_next_time_step_raw(n, t, y) + except Exception as exc: + raise StepEstimateError("Failed to get next time step") from exc + + if not np.isfinite(dt): + raise StepEstimateError(f"Time step is not finite: {dt}") + + # add eps to ensure that t >= tfinal is true + if self.tfinal is not None: + eps = float(5.0 * np.finfo(y.dtype).eps) + dt = min(dt, self.tfinal - t) + eps + + # TODO: Would be nice to have some smoothing over time steps so that + # they don't vary too much here, but that may be unwanted? - def predict_time_step(t: float, y: Array) -> float: return dt - return predict_time_step +@dataclass(frozen=True) +class FixedTimeSpan(TimeSpan): + """A :class:`TimeSpan` with a fixed time step.""" + + #: The fixed time step that should be used. + dt: float + + def get_next_time_step_raw(self, n: int, t: float, y: Array) -> float: + """See :meth:`pycaputo.fode.TimeSpan.get_next_time_step_raw`.""" + return self.dt + + @classmethod + def from_data( + cls, + dt: float, + tstart: float = 0.0, + tfinal: float | None = None, + nsteps: int | None = None, + ) -> FixedTimeSpan: + """Create a consistent time span with a fixed time step. + + This ensures that the following relation holds: + ``tfinal = tstart + nsteps * dt`` for all given values. This can be + achieved by small modifications to either *nsteps* or *dt*. + + :arg dt: desired time step (chosen time step may be slightly smaller). + :arg tstart: start of the time span. + :arg tfinal: end of the time span, which is not required. + :arg nsteps: number of time steps in the span, which is not required. + """ + + if tfinal is not None: + nsteps = int((tfinal - tstart) / dt) + 1 + dt = (tfinal - tstart) / nsteps + elif nsteps is not None: + tfinal = tstart + nsteps * dt + else: + # nsteps and tfinal are None, so nothing we can do here + pass + + return cls(tstart=tstart, tfinal=tfinal, nsteps=nsteps, dt=dt) -def make_predict_time_step_graded( - tspan: tuple[float, float], maxit: int, r: int = 2 -) -> ScalarStateFunction: - r"""Construct a time step that is smaller around the initial time. + +@dataclass(frozen=True) +class GradedTimeSpan(TimeSpan): + r"""A :class:`TimeSpan` with a variable graded time step. This graded grid of time steps is described in [Garrappa2015b]_. It essentially takes the form .. math:: - \Delta t_n = \frac{T - t_0}{N^r} ((n + 1)^r - n^r), + \Delta t_n = \frac{t_f - t_s}{N^r} ((n + 1)^r - n^r), - where the time interval is :math:`[t_0, T]` and :math:`N` time steps are + where the time interval is :math:`[t_s, t_f]` and :math:`N` time steps are taken. This graded grid can give full second-order convergence for certain methods such as the Predictor-Corrector method (e.g. implemented by :class:`~pycaputo.fode.CaputoPECEMethod`). - - :arg tspan: time interval :math:`[t_0, T]`. - :arg maxit: maximum number of iterations to take in the interval. - :arg r: a grading exponent that controls the clustering of points at :math:`t_0`. - - :returns: a callable returning a graded time step. """ - if maxit <= 0: - raise ValueError(f"Negative number of iterations not allowed: {maxit}") - if r < 1: - raise ValueError(f"Gradient exponent must be >= 1: {r}") + #: A grading exponent that controls the clustering of points at :math:`t_s`. + r: int - if tspan[0] > tspan[1]: - raise ValueError(f"Invalid time interval: {tspan}") + if __debug__: - h = (tspan[1] - tspan[0]) / maxit**r + def __post_init__(self) -> None: + if self.tfinal is None: + raise ValueError("'tfinal' must be given for the graded estimate.") - def predict_time_step(t: float, y: Array) -> float: - # FIXME: this works, but seems a bit overkill just to get the iteration - n = round(((t - tspan[0]) / h) ** (1 / r)) + if self.nsteps is None: + raise ValueError("'nsteps' must be given for the graded estimate") - return float(h * ((n + 1) ** r - n**r)) + super().__post_init__() + if self.r < 1: + raise ValueError(f"Exponent must be >= 1: {self.r}") - return predict_time_step + def get_next_time_step_raw(self, n: int, t: float, y: Array) -> float: + """See :meth:`pycaputo.fode.TimeSpan.get_next_time_step_raw`.""" + assert self.tfinal is not None + assert self.nsteps is not None + h = (self.tfinal - self.tstart) / self.nsteps**self.r + return float(h * ((n + 1) ** self.r - n**self.r)) -def make_predict_time_step_lipschitz( - f: StateFunction, - alpha: float, - yspan: tuple[float, float], - *, - theta: float = 0.9, -) -> ScalarStateFunction: - r"""Estimate a time step based on the Lipschitz constant. - This method estimates the time step using +@dataclass(frozen=True) +class FixedLipschitzTimeSpan(TimeSpan): + r"""A :class:`TimeSpan` that uses the Lipschitz constant to estimate the time step. + + This method estimates the time step using (Theorem 2.5 from [Baleanu2012]_) .. math:: - \Delta t = \theta \frac{1}{(\Gamma(2 - \alpha) L)^{\frac{1}{\alpha}}}, + \Delta t = \frac{1}{(\Gamma(2 - \alpha) L)^{\frac{1}{\alpha}}}, + + where :math:`L` is the Lipschitz constant estimated by :attr:`lipschitz_constant`. + Note that this is essentially a fixed time step estimate. + """ + + #: An estimate of the Lipschitz contants of the right-hand side :math:`f(t, y)` + #: for all times :math:`t \in [t_s, t_f]`. + lipschitz_constant: float + #: Fractional order of the derivative. + alpha: float + + def get_next_time_step_raw(self, n: int, t: float, y: Array) -> float: + """See :meth:`pycaputo.fode.TimeSpan.get_next_time_step_raw`.""" + if y.shape != (1,): + raise ValueError(f"Only scalar functions are supported: {y.shape}") + + from math import gamma - where :math:`L` is the Lipschitz constant estimated by - :func:`estimate_lipschitz_constant` and :math:`\theta` is a coefficient that - should be smaller than 1. + L = self.lipschitz_constant + return float((gamma(2 - self.alpha) * L) ** (-1.0 / self.alpha)) + + +@dataclass(frozen=True) +class LipschitzTimeSpan(TimeSpan): + r"""A :class:`TimeSpan` that uses the Lipschitz constant to estimate the time step. + + This uses the same logic as :class:`FixedLipschitzTimeSpan`, but computes the + Lipschitz constant using the estimate from + :func:`~pycaputo.lipschitz.estimate_lipschitz_constant` at each time :math:`t`. .. warning:: This method requires many evaluations of the right-hand side function *f* and will not be efficient in practice. Ideally, the Lipschitz constant is approximated by some theoretical result. - - :arg f: the right-hand side function used in the differential equation. - :arg alpha: fractional order of the derivative. - :arg yspan: an expected domain of the state variable. The Lipschitz constant - is estimated by sampling in this domain. - :arg theta: a coefficient used to control the time step. - - :returns: a callable returning an estimate of the time step based on the - Lipschitz constant. """ - from functools import partial - from math import gamma - if not 0.0 < alpha < 1.0: - raise ValueError(f"Derivative order must be in (0, 1): {alpha}") + #: The right-hand side function used in the differential equation. + f: StateFunction + #: Fractional order of the derivative. + alpha: float + #: An expected domain of the state variable. The Lipschitz constant is + #: estimated by sampling in this domain. + yspan: tuple[float, float] - if theta <= 0.0: - raise ValueError(f"Theta constant cannot be negative: {theta}") + if __debug__: + + def __post_init__(self) -> None: + super().__post_init__() - if yspan[0] > yspan[1]: - yspan = (yspan[1], yspan[0]) + if not 0.0 < self.alpha < 1.0: + raise ValueError(f"Derivative order must be in (0, 1): {self.alpha}") - from pycaputo.lipschitz import estimate_lipschitz_constant + if self.yspan[0] >= self.yspan[1]: + raise ValueError("Invalid yspan interval: ystart > yfinal") - def predict_time_step(t: float, y: Array) -> float: + def get_next_time_step_raw(self, n: int, t: float, y: Array) -> float: + """See :meth:`pycaputo.fode.TimeSpan.get_next_time_step_raw`.""" if y.shape != (1,): raise ValueError(f"Only scalar functions are supported: {y.shape}") - L = estimate_lipschitz_constant(partial(f, t), yspan[0], yspan[1]) - return float((gamma(2 - alpha) * L) ** (-1.0 / alpha)) + from math import gamma - return predict_time_step + from pycaputo.lipschitz import estimate_lipschitz_constant + + L = estimate_lipschitz_constant( + partial(self.f, t), self.yspan[0], self.yspan[1] + ) + return float((gamma(2 - self.alpha) * L) ** (-1.0 / self.alpha)) # }}} @@ -202,15 +329,11 @@ class FractionalDifferentialEquationMethod(ABC): #: The fractional derivative order used for the derivative. derivative_order: tuple[float, ...] - #: A callable used to predict the time step. - predict_time_step: float | ScalarStateFunction + #: An instance describing the discrete time span being simulated. + tspan: TimeSpan #: Right-hand side source term. source: StateFunction - #: The initial and final times of the simulation. The final time can be - #: *None* if evolving the equation to an unknown final time, e.g. by - #: setting *maxit* in :func:`evolve`. - tspan: tuple[float, float | None] #: Values used to reconstruct the required initial conditions. y0: tuple[Array, ...] @@ -237,11 +360,6 @@ def __post_init__(self) -> None: def largest_derivative_order(self) -> float: return max(self.derivative_order) - @property - def is_constant_time_step(self) -> bool: - """A flag for whether the method uses a constant time step.""" - return not callable(self.predict_time_step) - @property def name(self) -> str: """An identifier for the method.""" @@ -264,7 +382,6 @@ def evolve( *, callback: CallbackFunction | None = None, history: History | None = None, - maxit: int | None = None, verbose: bool = True, ) -> Iterator[Event]: """Evolve the fractional-order ordinary differential equation in time. @@ -274,7 +391,6 @@ def evolve( :arg y0: initial conditions for the FODE. :arg history: a :class:`History` instance that handles checkpointing the necessary state history for the method *m*. - :arg maxit: the maximum number of iterations. :arg verbose: print additional iteration details. :returns: an :class:`Event` (usually a :class:`StepCompleted`) containing @@ -292,7 +408,7 @@ def advance( ) -> Array: """Advance the solution *y* by to the time *t*. - This function takes ``(history[t_0, ... t_n], t_{n + 1}, y_n)`` and is + This function takes ``(history[t_s, ... t_n], t_{n + 1}, y_n)`` and is expected to update the history with values at :math:`t_{n + 1}`. The time steps can be recalculated from the history if necessary. diff --git a/pycaputo/fode/caputo.py b/pycaputo/fode/caputo.py index 4ee97b4..b7e80c4 100644 --- a/pycaputo/fode/caputo.py +++ b/pycaputo/fode/caputo.py @@ -76,7 +76,7 @@ def _advance_caputo_forward_euler( alpha = m.derivative_order dy = np.zeros_like(y) - dy = _update_caputo_initial_condition(dy, t - m.tspan[0], m.y0) + dy = _update_caputo_initial_condition(dy, t - m.tspan.tstart, m.y0) dy = _update_caputo_forward_euler(dy, history, alpha, len(history)) history.append(t, m.source(t, dy)) @@ -269,7 +269,7 @@ def _advance_caputo_weighted_euler( # add explicit terms fnext = np.zeros_like(y) - fnext = _update_caputo_initial_condition(fnext, t - m.tspan[0], m.y0) + fnext = _update_caputo_initial_condition(fnext, t - m.tspan.tstart, m.y0) fnext, omega = _update_caputo_weighted_euler(fnext, history, alpha, m.theta, n) # solve implicit equation @@ -447,7 +447,7 @@ def _advance_caputo_predictor_corrector( # add initial conditions y0 = np.zeros_like(y) - y0 = _update_caputo_initial_condition(y0, t - m.tspan[0], m.y0) + y0 = _update_caputo_initial_condition(y0, t - m.tspan.tstart, m.y0) # predictor step (forward Euler) yp = np.copy(y0) @@ -512,7 +512,7 @@ def _advance_caputo_modified_pece( # compute common terms dy = np.zeros_like(y) - dy = _update_caputo_initial_condition(dy, t - m.tspan[0], m.y0) + dy = _update_caputo_initial_condition(dy, t - m.tspan.tstart, m.y0) if n == 1: yp = _update_caputo_forward_euler(dy, history, m.derivative_order, len(history)) diff --git a/pycaputo/fode/history.py b/pycaputo/fode/history.py index c2006ed..a37015d 100644 --- a/pycaputo/fode/history.py +++ b/pycaputo/fode/history.py @@ -5,6 +5,9 @@ from abc import ABC, abstractmethod from dataclasses import dataclass, field +from typing import Any + +import numpy as np from pycaputo.logging import get_logger from pycaputo.utils import Array @@ -12,6 +15,19 @@ logger = get_logger(__name__) +def make_history(state: Array, nsteps: int | None = None) -> History: + """Make a history checkpointing instance for the given *state* type. + + :arg state: an array of the same shape and :class:`~numpy.dtype` as the + results generated by the time stepper. + :arg nsteps: number of time steps that will be taken, if known in advance. + """ + if nsteps is None: + return VariableProductIntegrationHistory() + else: + return FixedSizeHistory.empty(nsteps, state.shape, dtype=state.dtype) + + # {{{ interface @@ -131,7 +147,7 @@ class FixedSizeHistory(History): history: Array = field(repr=False) ts: Array = field(repr=False) - filled: int = 0 + filled: int = field(default=0, init=False) def __bool__(self) -> bool: return self.filled == 0 @@ -158,5 +174,23 @@ def append(self, t: float, y: Array) -> None: self.history[k] = y self.ts[k] = t + @classmethod + def empty( + cls, + n: int, + shape: tuple[int, ...], + *, + dtype: np.dtype[Any] | None = None, + ) -> FixedSizeHistory: + """Construct a :class:`FixedSizeHistory` for given sizes. + + :arg n: number of time steps that will be stored. + :arg d: dimension of each state array. + """ + return cls( + history=np.empty((n, *shape), dtype=dtype), + ts=np.empty(n, dtype=dtype), + ) + # }}} diff --git a/pycaputo/fode/product_integration.py b/pycaputo/fode/product_integration.py index 6800fd4..d378929 100644 --- a/pycaputo/fode/product_integration.py +++ b/pycaputo/fode/product_integration.py @@ -13,6 +13,7 @@ from pycaputo.fode.base import ( Event, FractionalDifferentialEquationMethod, + StepEstimateError, evolve, make_initial_condition, ) @@ -40,23 +41,12 @@ def _evolve_pi( *, callback: CallbackFunction | None = None, history: History | None = None, - maxit: int | None = None, raise_on_fail: bool = False, ) -> Iterator[Event]: - from pycaputo.fode.base import ( - StepCompleted, - StepFailed, - advance, - make_predict_time_step_fixed, - ) - - if callable(m.predict_time_step): - predict_time_step = m.predict_time_step - else: - predict_time_step = make_predict_time_step_fixed(m.predict_time_step) + from pycaputo.fode.base import StepCompleted, StepFailed, advance n = 0 - t, tfinal = m.tspan + t = m.tspan.tstart y = make_initial_condition(m) if history is None: @@ -73,36 +63,21 @@ def _evolve_pi( if callback is not None and callback(t, y): break - if tfinal is not None and t >= tfinal: + if m.tspan.finished(n, t): break - if maxit is not None and n >= maxit: - break - - # next iteration - n += 1 - # next time step try: - dt = predict_time_step(t, y) - except Exception as exc: + dt = m.tspan.get_next_time_step(n, t, y) + except StepEstimateError as exc: logger.error("Failed to predict time step.", exc_info=exc) if raise_on_fail: raise exc yield StepFailed(t=t, iteration=n, reason=str(exc)) - if not np.isfinite(dt): - logger.error("Predicted time step is not finite: %g", dt) - if raise_on_fail: - raise ValueError(f"Predicted time step is not finite: {dt}") - - yield StepFailed(t=t, iteration=n, reason="time step is not finite") - - if tfinal is not None: - # NOTE: adding eps to ensure that t >= tfinal is true - dt = min(dt, tfinal - t) + float(5 * np.finfo(y.dtype).eps) - + # next iteration + n += 1 t += dt # advance diff --git a/tests/test_fode_caputo.py b/tests/test_fode_caputo.py index 4d84e00..fc76e1b 100644 --- a/tests/test_fode_caputo.py +++ b/tests/test_fode_caputo.py @@ -65,31 +65,31 @@ def garrappa2009_source_jac(t: float, y: Array, *, alpha: float) -> Array: # }}} -# {{{ test_predict_time_step_graded +# {{{ test_graded_time_span -def test_predict_time_step_graded() -> None: - from pycaputo.fode import make_predict_time_step_graded +def test_graded_time_span() -> None: + from pycaputo.fode import GradedTimeSpan - maxit = 100 - tspan = (-1.5, 3.0) + nsteps = 100 + tstart, tfinal = (-1.5, 3.0) r = 3 - predict_time_step = make_predict_time_step_graded(tspan, maxit, r) + tspan = GradedTimeSpan(tstart, tfinal, nsteps=nsteps, r=r) - n = np.arange(maxit) - t_ref = tspan[0] + (n / maxit) ** r * (tspan[1] - tspan[0]) + n = np.arange(nsteps) + t_ref = tstart + (n / nsteps) ** r * (tfinal - tstart) t = np.empty_like(t_ref) dummy = np.empty(3) - t[0] = tspan[0] - for i in range(1, maxit): - dt = predict_time_step(t[i - 1], dummy) + t[0] = tstart + for i in range(1, nsteps): + dt = tspan.get_next_time_step(i - 1, t[i - 1], dummy) t[i] = t[i - 1] + dt error = la.norm(t - t_ref) / la.norm(t_ref) logger.info("error: %.12e", error) - assert error < 1.0e-15 + assert error < 1.0e-13 # }}} @@ -135,9 +135,8 @@ def wrapper( return cls( derivative_order=alpha, - predict_time_step=dt, + tspan=fode.FixedTimeSpan.from_data(dt, tstart=tspan[0], tfinal=tspan[1]), source=partial(source, alpha=alpha), - tspan=tspan, y0=(y0,), **kwargs, )