Update initial condition in ODE solver each time s

Question

I am wanting to solve a system of ODEs where for the first 30,000 seconds, I want one of my state variables to start from the same initial value. After those 30,000 seconds, I want to change the initial value of that state variable to something different and simulate the system for the rest of time. Here is my code:

def ode_rhs(y, t):
    ydot[0] = -p[7]*y[0]*y[1] + p[8]*y[8] + p[9]*y[8]
    ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]
    ydot[2] = -p[10]*y[2]*y[3] + p[11]*y[9] + p[12]*y[9]
    ydot[3] = -p[13]*y[3]*y[6] + p[14]*y[10] + p[15]*y[10] - p[10]*y[2]*y[3] + p[11]*y[9] + p[9]*y[8] - p[21]*y[3]
    ydot[4] = -p[19]*y[4]*y[5] - p[16]*y[4]*y[5] + p[17]*y[11] - p[23]*y[4] + y[7]*p[20]
    ydot[5] = -p[19]*y[4]*y[5] + p[15]*y[10] - p[16]*y[4]*y[5] + p[17]*y[11] + p[18]*y[11] + p[12]*y[9] - p[22]*y[5]
    ydot[6] = -p[13]*y[3]*y[6] + p[14]*y[10] - p[22]*y[6] - p[25]*y[6] - p[23]*y[6]
    ydot[7] = 0
    ydot[8] = p[7]*y[0]*y[1] - p[8]*y[8] - p[9]*y[8]
    ydot[9] = p[10]*y[2]*y[3] - p[11]*y[9] - p[12]*y[9] - p[21]*y[9]
    ydot[10] = p[13]*y[3]*y[6] - p[14]*y[10] - p[15]*y[10] - p[22]*y[10] - p[21]*y[10] - p[23]*y[10]
    ydot[11] = p[19]*y[4]*y[5] + p[16]*y[4]*y[5] - p[17]*y[11] - p[18]*y[11] - p[22]*y[11] - p[23]*y[11]
    ydot[12] = p[22]*y[10] + p[22]*y[11] + p[22]*y[5] + p[22]*y[6] + p[21]*y[10] + p[21]*y[3] + p[21]*y[9] + p[24]*y[13] + p[25]*y[6] + p[23]*y[10] + p[23]*y[11] + p[23]*y[4] + p[23]*y[6]
    ydot[13] = p[15]*y[10] + p[18]*y[11] - p[24]*y[13]
    return ydot

pysb.bng.generate_equations(model)
alias_model_components()
p = np.array([k.value for k in model.parameters])
ydot = np.zeros(len(model.odes))
y0 = np.zeros(len(model.odes))
y0[0:7] = p[0:7]
t = np.linspace(0.0,1000000.0,100000)
r = odeint(ode_rhs,y0,t)

So, in other words, I want to set y0[1] to the same value (100) each time odeint is called for the first 30,000 seconds. I'm effectively trying to let the system equilibrate for an amount of time before inputing a signal into the system. I thought about doing something like if t < 30000: y0[1] = 100 as the first line of my ode_rhs() function, but I'm not quite sure that works.

Answer 1

It sounds like you want y1(t) to remain constant (with the value 100) for the equilibration phase. You can do this by ensuring that dy1(t)/dt = 0 during this phase. There are (at least) two ways you can accomplish that. The first is to modify the calculation of ydot[1] in ode_rhs as follows:

if t < 30000:
    ydot[1] = 0.0
else:
    ydot[1] = -p[7]*y[0]*y[1] + p[8]*y[8]

and use the intitial condition 100 for y[1].

Note that this introduces a discontinuity in the right-hand side of the system, but the adaptive solver used by odeint (the Fortran code LSODA) is usually robust enough to handle it.

Here's a self-contained example. I've made p and t1 arguments to ode_rhs. t1 is the duration of the equilibration phase.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def ode_rhs(y, t, p, t1):
    ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
    if t < t1:
        ydot[1] = 0.0
    else:
        ydot[1] = -p[0]*y[0]*y[1] + p[1]*y[2]
    ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
    return ydot


ydot = np.zeros(3)
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t = np.linspace(0, 200, 2001)
t1 = 20.0

sol = odeint(ode_rhs, y0, t, args=(p, t1))


plt.figure(1)
plt.clf()

plt.subplot(3, 1, 1)
plt.plot(t, sol[:, 0])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[0]')


plt.subplot(3, 1, 2)
plt.plot(t, sol[:, 1])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)

plt.subplot(3, 1, 3)
plt.plot(t, sol[:, 2])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')

plt.show()

A slight variation of the above method is to modify the system by adding a parameter that is either 0 or 1. When the parameter is 0, the equlibration system is solved, and when the parameter is 1, the full system is solved. The code for ydot[1] (in my smaller example) is then

ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])

where full is the parameter.

To handle the equilibration phase, the system is solved once on 0 <= t < t1 with full=0. Then the final value of the equilibration solution is used as the initial condition to the second solution, run with full=1. The advantage of this method is that you are not forcing the solver to deal with the discontinuity.

Here's how it looks in code.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt


def ode_rhs(y, t, p, full):
    ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
    ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
    ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
    return ydot


ydot = np.zeros(3)
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t1 = 20.0  # Equilibration time
tf = 200.0  # Final time

# Solve the equilibration phase.
teq = np.linspace(0, t1, 100)
full = 0
soleq = odeint(ode_rhs, y0, teq, args=(p, full))

# Solve the full system, using the final point of the
# equilibration phase as the initial condition.
y0 = soleq[-1]
# Note: the system is autonomous, so we could just as well start
# at t0=0.  But starting at t1 makes the plots (below) align without
# any additional shifting of the time arrays.
t = np.linspace(t1, tf, 2000)
full = 1
sol = odeint(ode_rhs, y0, t, args=(p, full))

plt.figure(2)
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(teq, soleq[:, 0], t, sol[:, 0])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[0]')

plt.subplot(3, 1, 2)
plt.plot(teq, soleq[:, 1], t, sol[:, 1])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)

plt.subplot(3, 1, 3)
plt.plot(teq, soleq[:, 2], t, sol[:, 2])
plt.axvline(t1, color='r')
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')

plt.show()

And here's the plot that it generates (the plot from the first example is almost exactly the same):

Answer 2

This answer did not work me, as I needed to alter my initial conditions periodically. Therefore, I would like to propose alternative solution, which is to alternate conditions for differential equation inside the function itself:

We look at the value of t and adjust it:

if int(t) > 20:
    full = 1
else:
    full = 0

Here it is inside a function:

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

def ode_rhs(y, t, p, full):

    if int(t) > 20:
        full = 1
    else:
        full = 0

    ydot[0] = -p[0]*y[0]*y[1] + p[1]*y[2] + p[2]*y[2]
    ydot[1] = full * (-p[0]*y[0]*y[1] + p[1]*y[2])
    ydot[2] = p[0]*y[0]*y[1] - p[1]*y[2] - p[2]*y[2]
    return ydot

ydot = np.zeros(3)

# intial conditions
p = np.array([0.01, 0.25, 0.1])
y0 = [20.0, 100.0, 0.0]
t = np.linspace(0, 200, 100)
full = 0

# solve equation
solution = odeint(ode_rhs, y0, t, args=(p, full))

plt.figure()
plt.clf()
plt.subplot(3, 1, 1)
plt.plot(t, solution[:, 0])
plt.axvline(20, color='r')  # vertical line
plt.grid(True)
plt.ylabel('y[0]')

plt.subplot(3, 1, 2)
plt.plot(t, solution[:, 1])
plt.axvline(20, color='r')  # vertical line
plt.grid(True)
plt.ylabel('y[1]')
plt.ylim(0, 110)

plt.subplot(3, 1, 3)
plt.plot(t, solution[:, 2])
plt.axvline(20, color='r')  # x=20 vertical line
plt.grid(True)
plt.ylabel('y[2]')
plt.xlabel('t')

plt.show()

That allows to call the function to solve the equation exactly once.

• You don't have to mess with initial conditions from the previous step
• Easier to plot
• Code is more cleaner and easier to manage

And more importantly, you can now adjust parameters periodically inside the equation. For instance, say you have t = [0:200] and you want to change value of full every 20 steps, you can do so like this:

if int(t/20) % 2 == 0:
    full = 1
else:
    full = 0