hgrecco / numbakit-ode

Leveraging numba to speed up ODE integration
Other
68 stars 3 forks source link

Calling solver.run twice produce the error: Cannot interpolate at t=0.0 as it is smaller than the current smallest value in history (9.0) #29

Closed ivvaan closed 5 months ago

ivvaan commented 12 months ago
import nbkode
def func(t, y):
    return -0.1 * y
    t0 = 0.
    y0 = 1.
solver = nbkode.ForwardEuler(func, t0, y0)
ts, ys = solver.run([0., 5., 10.])
ts, ys = solver.run([0., 5., 10.])

ValueError                                Traceback (most recent call last)
Cell In[126], line 7
      5 solver = nbkode.ForwardEuler(func, t0, y0)
      6 ts, ys = solver.run([0., 5., 10.])
----> 7 ts, ys = solver.run([0., 5., 10.])

File ~\AppData\Roaming\Python\Python311\site-packages\nbkode\core.py:333, in Solver.run(self, t)
    316 def run(self, t: Union[Real, np.ndarray]) -> Tuple[np.ndarray, np.ndarray]:
    317     """Integrates the ODE interpolating at each of the timepoints `t`.
    318 
    319     Parameters
   (...)
    331         One of the timepoints provided is outside the valid range.
    332     """
--> 333     return self.run_events(t, None)[:2]

File ~\AppData\Roaming\Python\Python311\site-packages\nbkode\core.py:395, in Solver.run_events(self, t, events)
    392     t = t[ndx]
    394 if t[0] < self.cache.ts[0]:
--> 395     raise ValueError(
    396         f"Cannot interpolate at t={t[0]} as it is smaller "
    397         f"than the current smallest value in history ({self.cache.ts[0]})"
    398     )
    400 self._check_time(np.max(t))
    402 to_interpolate = t <= self.t

ValueError: 

> Cannot interpolate at t=0.0 as it is smaller than the current smallest value in history (9.0)
maurosilber commented 12 months ago

Each solver instance has an internal state (the last time t and state y), and successive calls to .run start from the last time. So, you can do this:

solver = nbkode.ForwardEuler(func, t0, y0)
ts1, ys1 = solver.run([0., 5., 10.])
ts2, ys2 = solver.run([15., 20., 25.])

but you have to create a new instance to "restart" it:

solver = nbkode.ForwardEuler(func, t0, y0)
ts, ys = solver.run([0., 5., 10.])
solver = nbkode.ForwardEuler(func, t0, y0)
ts, ys = solver.run([0., 5., 10.])