Trouble with backwards time integration in Python

Question

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.

A few comments/questions:

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.

Answer 1

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,