# Trouble with backwards time integration in Python

I am struggling with a rather basic numerical integration task: Using Python's scipy.integrate.solve_ivp module to integrate an ODE sytem backwards in time. As a test, I am using the following ODE system (an SEIR model for an infectious disease): \begin{align*} \frac{dS}{dt} &= - \beta SI \\[4pt] \frac{dE}{dt} &= \beta SI - \sigma E \\[4pt] \frac{dI}{dt} &= \sigma E - MI \\[4pt] \frac{dR}{dt} &= MI \end{align*}

To verify that the backwards time integration works, I first integrated it forwards in time, using the initial condition $$y(0) = (S(0), E(0), I(0), R(0)) = (58500, 800, 200, 500)$$ and the parameter values $$(\beta, \sigma, M) = (.1493, 0.1917, 0.2016).$$ The solution for $$I(t)$$ looks as follows:

I used the forward-time solution value at $$t = 70$$ as the "initial" condition for the backwards-time integration. I found that $$y(70) = (56349.39, 54.42, 61.62, 3534.57)$$. If my understanding of integration is correct, integrating backwards from $$t = 70$$ to $$t = 0$$, with $$y(70)$$ as above, should give us $$y(0) = (58500, 800, 200, 500)$$. Unfortunately, my attempt gave some very different values (see the graph and time series values at the end of the post). Here is my code in Python:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import pandas as pd

def SEIR_bw(t,y,params):

S,E,I,R = y
beta,sigma,M = params

dS = -beta*S*I
dE = beta*S*I - sigma*E
dI = sigma*E - M*I
dR = M*I

slope = -1*np.asarray([dS, dE, dI, dR])
return slope

#The values at t = 70
N0 = 60000
S70 = 56349.385363
E70 = 54.421832
I70 = 61.622725
R70 = 3534.570079

#The "initial" condition (scaled by N0)
IC = np.asarray([S70, E70, I70, R70])/N0

#The parameter values
beta = .1493
sigma = 0.1917
M = 0.2016

params = np.asarray([beta, sigma, M])
t_vals = np.arange(70,-1,-1)

#Perform the backwards in time integration
out = solve_ivp(fun = SEIR_bw, t_span = [70,0], y0 = IC, args = (params,),
t_eval = t_vals, method = 'RK45')

#Put soln in data frame and scale values back to population size
y = N0*pd.DataFrame(out.y).T
y.columns = ['S','E','I','R']
y.insert(0,'t',t_vals)

print(y)

plt.figure()
plt.plot(y['t'],y['I'])
plt.xlabel(r'$$t$$')
plt.title(r'$$I(t)$$')


Clearly I am doing something wrong because this is the plot it produced:

What am I doing wrong here? Either my understanding of integration is wrong, or my understanding of scipy.integrate.solve_ivp is wrong--I'm not sure which...

Note that in the code I scaled the values of $$S,E,I,R$$ by the total population size $$N(0) = 60,000$$ before performing the integration.

1. Note that I passed the function $$-dy/dt$$ to the integrator. I figure that taking a backwards time step means the slope is the opposite sign of what it would be for a forwards time step. Is this correct?
2. Is t_span = [70,0] the correct way to pass this argument for backwards integration?
3. Some values from the solution are shown below. Note that the y(70) value is correct, but the rest doesn't seem to make sense...
     t             S          E          I            R
0   70  56349.385363  54.421832  61.622725  3534.570079
1   69  56340.883713  52.660602  59.661170  3546.794514
2   68  56332.654245  50.959057  57.757360  3558.629337
3   67  56324.689528  49.312579  55.912825  3570.085066
4   66  56316.990018  47.699918  54.149199  3581.160863
..  ..           ...        ...        ...          ...
66   4  56116.234714   6.188500   7.091672  3870.485113
67   3  56115.244935   6.029139   6.811800  3871.914125
68   2  56114.285489   5.878736   6.536408  3873.299367
69   1  56113.365513   5.715364   6.291479  3874.627642
70   0  56112.494673   5.516053   6.104277  3875.884996


If anyone can shed some light on why my results aren't making sense (and how to fix the code accordingly), I would greatly appreciate it.

• Is your parameter beta correct conceptually? The size of it us usually correct for density simulations, for population number simulations you would need to divide the corresponding term by N, either in the beta initialization or in the ODE function. Dec 18 '20 at 13:28

You were trying to be too helpful. While it is true that for $$u(t)=y(T-t)$$ you get the equation $$\dot u(t)=-\dot y(T-t)=-F(u(t))$$ for an autonomous system, this means that to get the correct solution you now need to integrate the modified bw system forward in time, with increasing time points.

The second variant that you used is simply to call the solver with a descending time array. Then the step size is negative, which already performs the reversion of the derivatives vector. So no modification of the ODE function is necessary.

In total, by applying both methods, what you effectively did was continuing the forward solution for the time span $$[70,140]$$. The plot is for the reversed time $$[70,0]$$, thus you get an ascending graph instead of a descending one.

With unchanged beta you get a non-sensical fall of S to zero immediately at the start of the integration. Change that to beta = .1493/N0 at the initialization. Then the integration calls

T=70.0
t_vals = np.arange(0,T+0.1,1)
atol, rtol = 1e-15, 1e-12
#Perform the forwards-in-time integration
out1 = solve_ivp(fun = SEIR, t_span = [0,T], y0 = IC, args = (params,),
t_eval = t_vals, method = 'RK45', atol=atol, rtol=rtol)
print(out1.message)
print(list(out1.y[:,-1]))
#Perform the backwards-in-time integration
out2 = solve_ivp(fun = SEIR, t_span = [T,0], y0 = out1.y[:,-1], args = (params,),
t_eval = t_vals[::-1], method = 'RK45', atol=atol, rtol=rtol)
print(out2.message)
print(list(out2.y[:,-1]))


give the output

The solver successfully reached the end of the integration interval.
[56349.39797738794, 54.377880676416005, 61.67579473861903, 3534.548347196981]
The solver successfully reached the end of the integration interval.
[58498.29977697898, 803.591935779938, 195.75357247214114, 502.3547147688503]


and the plot for the I component, blue forward, red backward,

• Thank you, Lutz! Your explanation makes a lot of sense! There is one thing that is still bothering me: I copied your code and reproduced your numbers (and using beta/N0 as you did, not scaling the initial conditions). But then I tried replacing y0 = out1.y[:,-1] with [56349.39,54.42,61.62,3534.57] in the second call to solve_ivp() (just truncating a few decimals), and I got the following output: "Required step size is less than spacing between numbers. [86717973.99119039, -90012490.01218799, 3945560.6281798254, -591044.6071822459]..." Huh? Why does it produce nonsense now? Dec 18 '20 at 21:29
• Because in the forward direction the non-linear terms give a negative feed-back, forcing the system towards the equilibrium point, the same term in the backwards direction gives positive feed-back, making the solution path unstable. You get similar errors when you take the exact last point but increase the error tolerances even by 2 orders. These large values and the possible divergence then make the system stiff, forcing the step size down. This problem might be reduced a little by using an implicit solver method like Radau or LSODA. Dec 18 '20 at 21:34
• Thanks for the clarification. (Though this is unfortunate news for the ideas I had in mind with this model...) I'll have to read up on stiff ODEs, so I can get a better understanding of the phenomenon you described, as I admittedly don't understand the concept that well. Dec 18 '20 at 21:40
• No one understands it really well. Essentially, there is a region of step sizes or error tolerances where the combination of ODE system, method and numerical number type works like the leading error term prescribes. This whole situation is stiff if that working range is too small, only one or two doublings large. 64bit floating point is practically fixed, and this problem is much more pronounced with explicit methods, thus changing to implicit often helps somewhat. Dec 18 '20 at 22:24