Repository: scipy/scipy
Pull request: BUG: Fix LSODA interpolation scheme
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->
Closes gh-14480

#### What does this implement/fix?
<!--Please explain your changes.-->
As explained in #14480, `LsodaDenseOutput` fails to reproduce `y(t_old)`. It looks to me that this only occurs when LSODA plans on changing the method order (i.e. `NQU != NQCUR`, equivalent to `self._lsoda_solver._integrator.iwork[13] != self._lsoda_solver._integrator.iwork[14]`). To be clear, this is not exactly a bug in Scipy - Scipy perfectly implements Taylor's theorem and applies the scaling explained in DLSODA. Instead, it is likely a mistake either in the FORTRAN implementation or in its associated docstring.

Taking #14480 again as an example, the implemented formula in Scipy is equivalent to:
```
(1.71524262e-06
 + -1.24130739e-06 * (0.06342019780923039 - 0.06402588986309915) / 0.00019480871245159827
 + 1.88400737e-07 * ((0.06342019780923039 - 0.06402588986309915) / 0.00019480871245159827) ** 2)
```
which yields `7.395919849794794e-06`.

The proposed formula reads:
```
(1.71524262e-06
 + -1.24130739e-06 * (0.06342019780923039 - 0.06402588986309915) / 0.00019480871245159827
 + 1.88400737e-07 * ((0.06342019780923039 - 0.06402588986309915) / 0.00019480871245159827) ** 2
 + -4.87455807e-07 * ((0.06342019780923039 - 0.06402588986309915) / 0.0006056920538687616) ** 3)
```
and yields `7.883375656794793e-06`, which is the correct result.

In the formulas above, `[1.71524262e-06, -1.24130739e-06, 1.88400737e-07, -4.87455807e-07]` is the Nordsieck history array `YH`, obtained in that case through `self._lsoda_solver._integrator.rwork[20:24]`; `0.06342019780923039` is `t_old`; `0.06402588986309915` is `t`; `0.00019480871245159827` is `self._lsoda_solver._integrator.rwork[11]`, `HCUR`; and `0.0006056920538687616` is `self._lsoda_solver._integrator.rwork[10]`, `HU`.

Therefore, it looks like one term is missing when `NQCUR < NQU`. According to the docstring, the number of terms included in the formula should be controlled by `NQCUR`. However, I have found that it is really `NQU` that should be used. In the case where `NQU == NQCUR`, nothing changes. When `NQCUR < NQU`, as highlighted above, one term is added, and, interestingly, it must use `HU` instead of `HCUR`. Finally, when `NQCUR > NQU`, the proposed implementation discards one term compared to the current implementation. I have tested locally all the possible scenarios with both method order and step size changing; the proposed implementation yields the right result every time.

I do not have an explanation for the why behind all of this. Perhaps someone motivated can have a careful look through DSTODA - I feel that the reason lies there somewhere.

#### Additional information
<!--Any additional information you think is important.-->
ODEPACK is available [here](https://computing.llnl.gov/projects/odepack/software).

Problem statement:
LSODA implementation of dense output yields incorrect result
<!-- 
Thank you for taking the time to file a bug report. 
Please fill in the fields below, deleting the sections that 
don't apply to your issue. You can view the final output
by clicking the preview button above.
Note: This is a comment, and won't appear in the output.
-->

I am using solve_ivp with LSODA to integrate an ODE in the presence of an event; the event function and the solve_ivp call look like
```
def interrupt_integration(t, y, *tuple_of_args):
        return y[0] - initial_value / 1e5

interrupt_integration.terminal = True
interrupt_integration.direction = -1

sol = solve_ivp(func, integration_range, [initial_value], method="LSODA", t_eval=t_values,
                events=interrupt_integration, args=tuple_of_args,
                first_step=1e-9, max_step=np.inf, rtol=5e-4, atol=5e-7,
                jac=jac, lband=0, uband=0, min_step=1e-11)
```
I am experiencing an issue following the occurrence of the event. Specifically, the value of y calculated by LSODA differs from the one returned by the call to the interpolant generated by LsodaDenseOutput. I obtain the following
```
0.002851754215052383 [1.21218232e-05] 1.2e-05 1.218231563465635e-07
0.0029674726381464737 [7.54816189e-06] 1.2e-05 -4.451838114260822e-06
0.002851754215052383 [1.19125649e-05] 1.2e-05 -8.743505742940739e-08
0.0029674726381464737 [7.54816189e-06] 1.2e-05 -4.451838114260822e-06
```
The values printed correspond to t, y, initial_value / 1e5 and y[0] - initial_value / 1e5. On the first line, one can see that the returned value from the event is positive, hence the event is inactive. On the second line, however, the returned value is negative and the event becomes active. As a result, solve_event_equation within ivp.py is called and the third and fourth lines are produced. The fourth line is identical to the second line, which is the expected behaviour but, in the third line, the value of y differs from that in the first line, meaning that the value calculated by the interpolant from LsodaDenseOutput is incorrect.

I have had a look at _dense_output_impl within lsoda.py, but this is beyond my knowledge. Does anyone have an idea of why the interpolant yields an erroneous value at the lower end of the interpolation range?

#### Scipy/Numpy/Python version information:
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```
-->
Python 3.9.6, Numpy 1.21.1, Scipy 1.7.0

Hints from issue discussion:
Hi,
do you have a minimal code sample that we can directly copy/paste to reproduce such an issue ?

Thanks for looking into this. I have sort of a reproducer, but considering some values within `func` are calculated through complex functions, I figured it would be best to just dump them into dictionaries. As a result, I have this:
```
import numpy as np
from scipy.integrate import solve_ivp


def func(t, y):
    return y / (1 - t) - y / (foo[t] - bar[t] * t)


def jac(t, y):
    return 1 / (1 - t) - 1 / (foo[t] - bar[t] * t)


def stop(t, y):
    print(t, y, y[0] - y0 / 1e5)
    return y[0] - y0 / 1e5


stop.terminal = True
stop.direction = -1

foo = {0.0: 0.029476812948995148, 1e-09: 0.02947681040077826,
       2e-09: 0.02947680785256066, 1.0002e-05: 0.02945128433533103,
       2.0002000000000003e-05: 0.029425678101783845,
       3.0002000000000002e-05: 0.029399989088615703,
       0.000130002: 0.029230452355455988,
       0.00023000199999999998: 0.029191774167678885,
       0.000330002: 0.029152915220043835,
       0.00043000199999999997: 0.029113875310867024,
       0.001430002: 0.02876766834368568, 0.002430002:
       0.02847185665295307, 0.003430002: 0.028202751722422884,
       0.004430002000000001: 0.027952662825450645,
       0.008519605788205807: 0.022905984940179824,
       0.007281063006042693: 0.027300979492089986,
       0.010132124012085386: 0.02025345412862362,
       0.008660010923042845: 0.021214768184557208,
       0.007556852589442723: 0.02724175078517937,
       0.007832642172842753: 0.02718479976564211,
       0.008384221339642813: 0.02527129401858587,
       0.008935800506442873: 0.019528292281649904,
       0.009487379673242933: 0.02037676066704014,
       0.00934000544054084: 0.01989416388865358,
       0.009744210374638808: 0.02032775867521952,
       0.010148415308736775: 0.020250325840727494,
       0.010866168034960004: 0.020113393229186118,
       0.011583920761183234: 0.019981633614217835,
       0.012301673487406463: 0.01984881183755286,
       0.014538351128642633: 0.019447584430548222,
       0.016775028769878804: 0.019058597729832413,
       0.01901170641111497: 0.018678455557851265,
       0.02423145283713465: 0.017814548563380436,
       0.02945119926315433: 0.016969990116954593,
       0.02759309012139257: 0.0172698529171475,
       0.03095472740565049: 0.016729649690162277,
       0.03431636468990841: 0.016192717307204542,
       0.03370112503574072: 0.01629065815929285,
       0.03644752266583095: 0.0158532553587958,
       0.03919392029592118: 0.015417187561523906,
       0.03864347589924936: 0.015504231290303634,
       0.04083942913266777: 0.015157383813710517,
       0.04303538236608618: 0.014809684936033558,
       0.042660331506766305: 0.014868650294859016,
       0.04448123388086484: 0.014581292474113251,
       0.04630213625496338: 0.014293458984378008,
       0.04951943560969117: 0.013785937182894536,
       0.048852300487073494: 0.0138912794370001,
       0.05140246471918361: 0.013488961724187341,
       0.05395262895129373: 0.013085220816138231,
       0.05650279318340384: 0.012679840416492304,
       0.05459017000932126: 0.012984291473056365,
       0.05522771106734879: 0.01288307272299171,
       0.055865252125376325: 0.01278156463465979,
       0.05650279318340386: 0.0126798404164923,
       0.057777875299458915: 0.012478571735661786,
       0.059052957415513974: 0.012276204069044047,
       0.06032803953156903: 0.012072821582693105,
       0.06160312164762409: 0.011870364583737523,
       0.0606468100605828: 0.012022308152748701,
       0.060965580589596564: 0.011971727504618139,
       0.06128435111861033: 0.011921079645771786,
       0.0616031216476241: 0.011870364583737521,
       0.062208813701492865: 0.011773816371676687,
       0.06281450575536163: 0.011677025613681216,
       0.06342019780923039: 0.011579992363122005,
       0.06402588986309915: 0.01148271667413238}
bar = {0.0: 0.0, 1e-09: 3.4086535954356014, 2e-09: 3.4086542978849566,
       1.0002e-05: 3.4159427995254577,
       2.0002000000000003e-05: 3.4232346887012572,
       3.0002000000000002e-05: 3.4305299687566397,
       0.000130002: 0.5756975084092468,
       0.00023000199999999998: 0.5771698633237826,
       0.000330002: 0.5786426287815926,
       0.00043000199999999997: 0.580115804987956,
       0.001430002: 0.5024993883592048,
       0.002430002: 0.4589098021277861,
       0.003430002: 0.430273623161857,
       0.004430002000000001: 0.4013637154514233,
       0.008519605788205807: 0.29037630343847937,
       0.007281063006042693: 0.3592207609391992,
       0.010132124012085386: 0.24566161928918345,
       0.008660010923042845: 0.2685977026683152,
       0.007556852589442723: 0.3570232800585041,
       0.007832642172842753: 0.35153277854189896,
       0.008384221339642813: 0.3207871254518165,
       0.008935800506442873: 0.2468874833078626,
       0.009487379673242933: 0.25151150476846235,
       0.00934000544054084: 0.24708375042521946,
       0.009744210374638808: 0.24843356077119336,
       0.010148415308736775: 0.24568640999109434,
       0.010866168034960004: 0.2452380786270683,
       0.011583920761183234: 0.2388955900978546,
       0.012301673487406463: 0.2366961569546367,
       0.014538351128642633: 0.22747068350665495,
       0.016775028769878804: 0.2209836035536901,
       0.01901170641111497: 0.2150464347549513,
       0.02423145283713465: 0.2052438158186673,
       0.02945119926315433: 0.19818362781401422,
       0.02759309012139257: 0.20039935988879123,
       0.03095472740565049: 0.1965611123586582,
       0.03431636468990841: 0.19323644827923409,
       0.03370112503574072: 0.19374659715079648,
       0.03644752266583095: 0.19099864724820126,
       0.03919392029592118: 0.18874966613397598,
       0.03864347589924936: 0.18944350874581367,
       0.04083942913266777: 0.1873057119526678,
       0.04303538236608618: 0.1854254416360881,
       0.042660331506766305: 0.1859252389311718,
       0.04448123388086484: 0.18458656503930682,
       0.04630213625496338: 0.1831337661765724,
       0.04951943560969117: 0.18076074209302737,
       0.048852300487073494: 0.1812626936758732,
       0.05140246471918361: 0.1794895940007729,
       0.05395262895129373: 0.17846279401370146,
       0.05650279318340384: 0.17759866407692976,
       0.05459017000932126: 0.17807580722623323,
       0.05522771106734879: 0.17789876937610216,
       0.055865252125376325: 0.17775444177480337,
       0.05650279318340386: 0.17759866407692948,
       0.057777875299458915: 0.17668665156727248,
       0.059052957415513974: 0.1759666871808635,
       0.06032803953156903: 0.17555528345585283,
       0.06160312164762409: 0.17483524999613184,
       0.0606468100605828: 0.17537267580783927,
       0.060965580589596564: 0.17519182162596608,
       0.06128435111861033: 0.17501268983427154,
       0.0616031216476241: 0.17483524999613087,
       0.062208813701492865: 0.17450265408093713,
       0.06281450575536163: 0.17420312944304941,
       0.06342019780923039: 0.1741284993925215,
       0.06402588986309915: 0.17405523810115953}

y0 = 0.772
sol = solve_ivp(func, [0.0, 0.07], [y0], method="LSODA",
                events=stop, first_step=1e-9, max_step=np.inf, rtol=5e-4,
                atol=5e-7, jac=jac, lband=0, uband=0, min_step=1e-11)
```
This fails for me:
```
0.06342019780923039 [7.88337564e-06] 1.6337564465948222e-07
0.06402588986309915 [1.71524262e-06] -6.004757381989139e-06
0.06342019780923039 [7.39591984e-06] -3.240801624861349e-07
0.06402588986309915 [1.71524262e-06] -6.004757381989139e-06
Traceback (most recent call last):
  File "/home/thomas/mwe.py", line 152, in <module>
    sol = solve_ivp(func, [0.0, 0.07], [y0], method="LSODA",
  File "/home/thomas/.local/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py", line 601, in solve_ivp
    root_indices, roots, terminate = handle_events(
  File "/home/thomas/.local/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py", line 108, in handle_events
    roots = [solve_event_equation(events[event_index], sol, t_old, t)
  File "/home/thomas/.local/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py", line 108, in <listcomp>
    roots = [solve_event_equation(events[event_index], sol, t_old, t)
  File "/home/thomas/.local/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py", line 77, in solve_event_equation
    return brentq(lambda t: event(t, sol(t)), t_old, t,
  File "/home/thomas/.local/lib/python3.9/site-packages/scipy/optimize/zeros.py", line 776, in brentq
    r = _zeros._brentq(f, a, b, xtol, rtol, maxiter, args, full_output, disp)
ValueError: f(a) and f(b) must have different signs
```
where you can clearly see that y values associated to t = 0.06342019780923039 are not consistent.
Reading through the [FORTRAN code](https://computing.llnl.gov/projects/odepack/software), Scipy's implementation is correct. Using the values denoted T, TCUR and HCUR to calculate S, NQCUR + 1 to obtain L, and the YH values, DKY as coded in DINTDY yields the exact same value as given by Scipy. The mathematical formula is correct, but I am wondering if there could be an error somewhere in the scaling of YH; it is mentioned multiple times in DSTODA that YH has to be rescaled.

```
>>> S = (0.06342019780923039 - 0.06402588986309915) / 0.00019480871245159827
>>> L = 3
>>> YH = [1.71524262e-06, -1.24130739e-06, 1.88400737e-07]
>>> DKY = YH[-1]
>>> for JB in range(1, L):
...     DKY = [1.71524262e-06, -1.24130739e-06,  1.88400737e-07][2 - JB] + S * DKY
... 
>>> DKY
7.395919849794794e-06
```
Thanks Scipy for making really awesome stuff.
I also ran into this problem. LSODA Fortran intdy() routine uses yh data (rwork[20 to end] sort of) to interpolate solution (states) only for within the last time stoda subroutine was called, so the last solver step. Scipy has written an interpolant object which allows the user to provide a time value which can be earlier than the time range of the last stoda step (and no warning is provided to user), while intdy() will error out for the same use case. I even had an exchange with Alan Hindmarsh about the operational scope of LSODA solution interpolation and he confirmed it's only the last step. Since Scipy essentially only allows the user to interact with LSODA in its itask=1 mode (integrate to some future time t_out), Scipy should then keep track of the time span of the last step, then restrict the user interpolation requests to that time span. or, let the user interact with LSODA with itask=2 (single stoda step mode).
@albert-mathews Thank you for contributing to the discussion. I am not sure I understand why you assume the SciPy interpolant can reach outside of the timestep. Is there a different timestep between `LSODA` and `STODA`? If not, all that SciPy does is concatenate each step-specific interpolant as the integration progresses, which then allows for interpolation across the full integration range. The issue described here corresponds to an erroneous choice in terms of how many terms should be included in Taylor's theorem and an erroneous scaling of one of the required terms. Out of curiosity, what do you think about the proposed changes in the PR associated with this issue?
@Patol75 
> If not, all that SciPy does is concatenate each step-specific interpolant as the integration progresses, which then allows for interpolation across the full integration range.

thanks for the clarification, i was not aware of that. It was not clear to me from looking at the code. it looks like Scipy is just fetching the latest rwork data from LSODA memory:
![image](https://user-images.githubusercontent.com/80334291/134340758-f5773d51-3b47-4063-98da-faea9ad51a84.png)

I read the PR, and I cannot comment as I'm familiar with the C implementation of LSODA, and not much the Fortran implementation, though it's highly likely the C implementation would share bugs if any. What I do know is that when I remove the time range check on intdy(), i get bogus interpolation back. When write my outer loops so as to ensure I'm only ever making calls to intdy() such that i can guarantee that time is within the last stoda step, i never get inaccurate interpolation results.

The way I understand it, when LSODA is requested to advance to some arbitrary future time tf starting from t, it may take any number of single steps, each step is performed by subroutine STODA. Typically the last step taken passes tf, and LSODA uses intdy() to find and return the solution at tf. For each STODA subroutine step (one step of length h), there is yh data written to the working memory, and this is used by intdy() to interpolate the solution from the last step. If t=0.0, tf=1.0, and the last step was 0.05s in length, and ended at 1.02, then the yh data in memory is only valid for interpolation of the solution between 0.97s and 1.02s. If the Fortran user calls intdy() with t=0.8, intdy() returns an error. If the Scipy user calls sol.sol with t=0.8s, Scipy returns a value.
My understanding is that the yh data is specific to the method (functional iteration or chord method), the order of the method, and the time range of the step (tn - h to tn, where tn is time after the completed step, and h is the step size). So if Scipy is saving each version of yh data after each step and concatenating (thus growing yh history with each stoda step), I would think it would also have to save the order, and time range associated with each given yh data, and for an arbitrary time provided by the user, it would first have to identify the set of yh and order data associated with the time range that contains the time provided by the user. If this is what Scipy is doing, then it was not evident to me from looking at Scipy code that it was doing this, so sorry for missing that. 

Still, using sol.sol with LSODA, at least for me, results slightly inaccurate results in some cases (likely those for which the time request was not within the last stoda step), where as running the exact same code by switching solve_ivp to RK45, the problem goes away. If this is not the right place for reporting this, sorry for muddying the waters, please point me to the right place, and if it's helpful, ill provide the code to recreate. 

thanks again for the response, and of course for creating such amazing python libraries. I have benefitted so much from Scipy, I am truly grateful. cheers. 
That is valuable feedback, thanks for that.

Regarding the implementation in SciPy, I can give some details about `solve_ivp`. SciPy uses a [`while` loop](https://github.com/scipy/scipy/blob/97ea4e506c7a4c6fdd144c09e00522134dd64c94/scipy/integrate/_ivp/ivp.py#L575-L638) to make progress through the integration. The loop depends on the status of the integration, initially set to `None` and which can take the value `0`, `-1`, or `1` at any loop step if integration successfully reaches the requested final `t` value, fails or encounters a terminating event, respectively. At each loop step, `solver.step()` is called, and thereby, integration marches one step further. The rest of the loop mainly corresponds to dealing with keyword arguments provided by the user, such as `t_eval`, `dense_output` and `events`. What we are interested in here is `solver.step()` and `solver.dense_output()`, with the latter called when `dense_output=True`. `solver.step()`, if using `LSODA`, will call [`_step_impl`](https://github.com/scipy/scipy/blob/97ea4e506c7a4c6fdd144c09e00522134dd64c94/scipy/integrate/_ivp/lsoda.py#L140-L161) in `lsoda.py`, and the [FORTRAN code](https://github.com/scipy/scipy/tree/master/scipy/integrate/odepack) will run for one integration step. After the step completes, `solve_ivp` will call `solver.dense_output()` if requested, which, in turn, calls [`_dense_output_impl`](https://github.com/scipy/scipy/blob/97ea4e506c7a4c6fdd144c09e00522134dd64c94/scipy/integrate/_ivp/lsoda.py#L163-L172) in `lsoda.py`. As a result, `rwork` is effectively accessed after every integration step. The result of `_dense_output_impl` is then stored in `interpolants` in `solve_ivp` at each loop step, and the final full interpolant is built through `OdeSolution` after execution left the `while` loop. I hope this clarifies what is being done.

Considering that you too experience inaccurate results through `sol.sol` with `LSODA`, could you compile the [branch](https://github.com/Patol75/scipy/tree/Fix_LSODA_Dense_Output) associated with the linked PR and report if results improve? That would definitely be helpful.
@Patol75 
thanks  for the detailed explanation. this makes sense now. 

We have stopped development on our python implementation and have switched to C. I cannot promise I'll will have time to go back and do this test, but if I do, I'll report my findings here.

Commit message: BUG: Fix LSODA interpolant.

Patch to review:
diff --git a/scipy/integrate/_ivp/common.py b/scipy/integrate/_ivp/common.py
index 7f6d2331100a..271e73ea9905 100644
--- a/scipy/integrate/_ivp/common.py
+++ b/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
diff --git a/scipy/integrate/_ivp/ivp.py b/scipy/integrate/_ivp/ivp.py
index 3d515ac66bb0..08a9011484eb 100644
--- a/scipy/integrate/_ivp/ivp.py
+++ b/scipy/integrate/_ivp/ivp.py
@@ -667,9 +667,13 @@ def solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False,
 
     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
 
diff --git a/scipy/integrate/_ivp/lsoda.py b/scipy/integrate/_ivp/lsoda.py
index 3a31efb90238..616d874b8028 100644
--- a/scipy/integrate/_ivp/lsoda.py
+++ b/scipy/integrate/_ivp/lsoda.py
@@ -168,10 +168,12 @@ def _dense_output_impl(self):
         iwork = self._lsoda_solver._integrator.iwork
         rwork = self._lsoda_solver._integrator.rwork
 
-        order = iwork[14]
+        order = iwork[13]
         h = rwork[11]
         yh = np.reshape(rwork[20:20 + (order + 1) * self.n],
                         (self.n, order + 1), order='F').copy()
+        if iwork[14] < order:
+            yh[:, -1] *= (h / rwork[10]) ** order
 
         return LsodaDenseOutput(self.t_old, self.t, h, order, yh)
 
diff --git a/scipy/integrate/_ivp/tests/test_ivp.py b/scipy/integrate/_ivp/tests/test_ivp.py
index 2fb294823af5..cdf864285663 100644
--- a/scipy/integrate/_ivp/tests/test_ivp.py
+++ b/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],


Code review comments:
diff_hunk: @@ -168,10 +168,12 @@ def _dense_output_impl(self):
         iwork = self._lsoda_solver._integrator.iwork
         rwork = self._lsoda_solver._integrator.rwork
 
-        order = iwork[14]
+        order = iwork[13]
         h = rwork[11]
         yh = np.reshape(rwork[20:20 + (order + 1) * self.n],
                         (self.n, order + 1), order='F').copy()
+        if iwork[14] < order:
path: scipy/integrate/_ivp/lsoda.py
text: @user1:
Looks like it can be achieved simpler by assigning `h = rwork[10]` in this case. Am I right?

@author:
Without having a detailed look, I do not think so. What is happening here is that only the last term needs rescaling, so changing `h` would not properly account for that. This last-term problem is likely a consequence of a mistake somewhere in the associated Fortran code.

diff_hunk: @@ -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.
path: scipy/integrate/_ivp/common.py
text: @user1:
This comment should be updated for the new behavior since now it's only true when `alt_segment=False`. 
```suggestion
        # Here we preserve a certain symmetry that when t is in self.ts,
        # if alt_segment=False, then we prioritize a segment with a lower
        # index.
```

Feel free to edit to your liking.

diff_hunk: @@ -168,10 +168,12 @@ def _dense_output_impl(self):
         iwork = self._lsoda_solver._integrator.iwork
         rwork = self._lsoda_solver._integrator.rwork
 
-        order = iwork[14]
+        order = iwork[13]
         h = rwork[11]
         yh = np.reshape(rwork[20:20 + (order + 1) * self.n],
                         (self.n, order + 1), order='F').copy()
+        if iwork[14] < order:
+            yh[:, -1] *= (h / rwork[10]) ** order
path: scipy/integrate/_ivp/lsoda.py
text: @user1:
```suggestion
        # 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
```

@user1:
Good to have some detailed comments explaining what's happening. I didn't spend too much time on this, so feel free to edit to your liking.

Edit: I had some mistakes in how I originally wrote these comments. Updated now.