TST: Update segment selection and add a test for LSODA interpolant.

scipy/integrate/_ivp/common.py

@@ -144,13 +144,21 @@ class OdeSolution:
     interpolants : list of DenseOutput with n_segments elements
         Local interpolants. An i-th interpolant is assumed to be defined
         between ``ts[i]`` and ``ts[i + 1]``.
+    alt_segment : boolean
+        Requests the alternative interpolant segment selection scheme. At each
+        solver integration point, two interpolant segments are available. The
+        default (False) and alternative (True) behaviours select the segment
+        for which the requested time corresponded to ``t`` and ``t_old``,
+        respectively. This functionality is only relevant for testing the
+        interpolants' accuracy: different integrators use different
+        construction strategies.
 
     Attributes
     ----------
     t_min, t_max : float
         Time range of the interpolation.
     """
-    def __init__(self, ts, interpolants):
+    def __init__(self, ts, interpolants, alt_segment=False):
         ts = np.asarray(ts)
         d = np.diff(ts)
         # The first case covers integration on zero segment.
@@ -169,20 +177,19 @@ def __init__(self, ts, interpolants):
             self.t_min = ts[0]
             self.t_max = ts[-1]
             self.ascending = True
+            self.side = "right" if alt_segment else "left"
             self.ts_sorted = ts
         else:
             self.t_min = ts[-1]
             self.t_max = ts[0]
             self.ascending = False
+            self.side = "left" if alt_segment else "right"
             self.ts_sorted = ts[::-1]
 
     def _call_single(self, t):
         # Here we preserve a certain symmetry that when t is in self.ts,
         # then we prioritize a segment with a lower index.
-        if self.ascending:
-            ind = np.searchsorted(self.ts_sorted, t, side='left')
-        else:
-            ind = np.searchsorted(self.ts_sorted, t, side='right')
+        ind = np.searchsorted(self.ts_sorted, t, side=self.side)
 
         segment = min(max(ind - 1, 0), self.n_segments - 1)
         if not self.ascending:
@@ -215,10 +222,7 @@ def __call__(self, t):
         t_sorted = t[order]
 
         # See comment in self._call_single.
-        if self.ascending:
-            segments = np.searchsorted(self.ts_sorted, t_sorted, side='left')
-        else:
-            segments = np.searchsorted(self.ts_sorted, t_sorted, side='right')
+        segments = np.searchsorted(self.ts_sorted, t_sorted, side=self.side)
         segments -= 1
         segments[segments < 0] = 0
         segments[segments > self.n_segments - 1] = self.n_segments - 1

scipy/integrate/_ivp/tests/test_ivp.py

@@ -1,5 +1,5 @@
 from itertools import product
-from numpy.testing import (assert_, assert_allclose,
+from numpy.testing import (assert_, assert_allclose, assert_array_less,
                            assert_equal, assert_no_warnings, suppress_warnings)
 import pytest
 from pytest import raises as assert_raises
@@ -135,6 +135,18 @@ def sol_complex(t):
     return y.reshape((1, -1))
 
 
+def fun_dense_output_LSODA(t, y):
+    return y * (t - 2)
+
+
+def jac_dense_output_LSODA(t, y):
+    return t - 2
+
+
+def sol_dense_output_LSODA(t):
+    return np.exp(t ** 2 / 2 - 2 * t + np.log(0.05) - 6)
+
+
 def compute_error(y, y_true, rtol, atol):
     e = (y - y_true) / (atol + rtol * np.abs(y_true))
     return np.linalg.norm(e, axis=0) / np.sqrt(e.shape[0])
@@ -201,11 +213,7 @@ def test_integration():
         e = compute_error(yc, yc_true, rtol, atol)
         assert_(np.all(e < 5))
 
-        # LSODA for some reasons doesn't pass the polynomial through the
-        # previous points exactly after the order change. It might be some
-        # bug in LSOSA implementation or maybe we missing something.
-        if method != 'LSODA':
-            assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)
+        assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)
 
 
 def test_integration_complex():
@@ -538,9 +546,7 @@ def test_max_step():
             e = compute_error(yc, yc_true, rtol, atol)
             assert_(np.all(e < 5))
 
-            # See comment in test_integration.
-            if method is not LSODA:
-                assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)
+            assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)
 
             assert_raises(ValueError, method, fun_rational, t_span[0], y0,
                           t_span[1], max_step=-1)
@@ -584,9 +590,7 @@ def test_first_step():
             e = compute_error(yc, yc_true, rtol, atol)
             assert_(np.all(e < 5))
 
-            # See comment in test_integration.
-            if method is not LSODA:
-                assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)
+            assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)
 
             assert_raises(ValueError, method, fun_rational, t_span[0], y0,
                           t_span[1], first_step=-1)
@@ -711,6 +715,37 @@ def early_event(t, y):
         assert res.t_events[0][0] == 7
 
 
+def test_dense_output_LSODA():
+    rtol = 1e-3
+    atol = 1e-6
+    y0 = [0.05]
+    t_span = [-2, 2]
+    first_step = 1e-3
+    res = solve_ivp(fun_dense_output_LSODA, t_span, y0, method="LSODA",
+                    dense_output=True, first_step=first_step, max_step=1,
+                    rtol=rtol, atol=atol, jac=jac_dense_output_LSODA)
+
+    assert_equal(res.t[0], t_span[0])
+    assert_equal(res.t[-1], t_span[-1])
+    assert_allclose(first_step, np.abs(res.t[1] - t_span[0]))
+    assert res.t_events is None
+    assert res.success
+    assert_equal(res.status, 0)
+
+    y_true = sol_dense_output_LSODA(res.t)
+    e = compute_error(res.y, y_true, rtol, atol)
+    assert_array_less(e, 5)
+
+    tc = np.linspace(*t_span)
+    yc_true = sol_dense_output_LSODA(tc)
+    yc = res.sol(tc)
+
+    e = compute_error(yc, yc_true, rtol, atol)
+    assert_array_less(e, 5)
+
+    assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15)
+
+
 def test_no_integration():
     for method in ['RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA']:
         sol = solve_ivp(lambda t, y: -y, [4, 4], [2, 3],