BUG: Fix LSODA interpolation scheme (#14552)

* TST: Enable LSODA tests and add failing test.

* ENH: Select appropriate segment for each method.

* BUG: Fix LSODA interpolant.

* ENH: Document LSODA changes

Co-authored-by: Albert Steppi <albert.steppi@gmail.com>

* ENH: Update segment comment

Co-authored-by: Albert Steppi <albert.steppi@gmail.com>

---------

Co-authored-by: Albert Steppi <albert.steppi@gmail.com>

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,20 @@ 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')
+        # if alt_segment=False, then we prioritize a segment with a lower
+        # index.
+        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 +223,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/ivp.py

@@ -678,9 +678,13 @@ def fun(t, x, fun=fun):
 
     if dense_output:
         if t_eval is None:
-            sol = OdeSolution(ts, interpolants)
+            sol = OdeSolution(
+                ts, interpolants, alt_segment=True if method in [BDF, LSODA] else False
+            )
         else:
-            sol = OdeSolution(ti, interpolants)
+            sol = OdeSolution(
+                ti, interpolants, alt_segment=True if method in [BDF, LSODA] else False
+            )
     else:
         sol = None
 

scipy/integrate/_ivp/lsoda.py

@@ -177,10 +177,33 @@ def _dense_output_impl(self):
         iwork = self._lsoda_solver._integrator.iwork
         rwork = self._lsoda_solver._integrator.rwork
 
-        order = iwork[14]
+        # We want to produce the Nordsieck history array, yh, up to the order
+        # used in the last successful iteration. The step size is unimportant
+        # because it will be scaled out in LsodaDenseOutput. Some additional
+        # work may be required because ODEPACK's LSODA implementation produces
+        # the Nordsieck history in the state needed for the next iteration.
+
+        # iwork[13] contains order from last successful iteration, while
+        # iwork[14] contains order to be attempted next.
+        order = iwork[13]
+
+        # rwork[11] contains the step size to be attempted next, while
+        # rwork[10] contains step size from last successful iteration.
         h = rwork[11]
+
+        # rwork[20:20 + (iwork[14] + 1) * self.n] contains entries of the
+        # Nordsieck array in state needed for next iteration. We want
+        # the entries up to order for the last successful step so use the 
+        # following.
         yh = np.reshape(rwork[20:20 + (order + 1) * self.n],
                         (self.n, order + 1), order='F').copy()
+        if iwork[14] < order:
+            # If the order is set to decrease then the final column of yh
+            # has not been updated within ODEPACK's LSODA
+            # implementation because this column will not be used in the
+            # next iteration. We must rescale this column to make the
+            # associated step size consistent with the other columns.
+            yh[:, -1] *= (h / rwork[10]) ** order
 
         return LsodaDenseOutput(self.t_old, self.t, h, order, yh)
 

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_event_dense_output_LSODA(t, y):
+    return y * (t - 2)
+
+
+def jac_event_dense_output_LSODA(t, y):
+    return t - 2
+
+
+def sol_event_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,48 @@ def early_event(t, y):
         assert res.t_events[0][0] == 7
 
 
+def test_event_dense_output_LSODA():
+    def event_lsoda(t, y):
+        return y[0] - 2.02e-5
+
+    rtol = 1e-3
+    atol = 1e-6
+    y0 = [0.05]
+    t_span = [-2, 2]
+    first_step = 1e-3
+    res = solve_ivp(
+        fun_event_dense_output_LSODA,
+        t_span,
+        y0,
+        method="LSODA",
+        dense_output=True,
+        events=event_lsoda,
+        first_step=first_step,
+        max_step=1,
+        rtol=rtol,
+        atol=atol,
+        jac=jac_event_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.success
+    assert_equal(res.status, 0)
+
+    y_true = sol_event_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_event_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],