0

I am simulating chemical reactions. That involves simple ODE systems, which I solve with scipy.solve_ivp as shown below (here I convert species y[0] into y[1] at a fixed rate).

import matplotlib as plt
import scipy as sp
import numpy as np

def fun(t, y): 
    return [-0.1 * y[0], +0.1 * y[0]]
#return array with the ODEs of levels of state variable 0, 1...

t_span = [0, 100]
y0 = [5, 0] # initial values
nb_of_steps = 100
t_eval = np.arange(t_span[0], 
                   t_span[1], 
                   (t_span[1]-t_span[0])/nb_of_steps)

ode_solve = sp.integrate.solve_ivp(fun, 
                                   t_span,
                                   y0,
                                   method='RK45',
                                   t_eval=t_eval,
                                   dense_output=False,
                                   events=None,
                                   vectorized=False)

plt.plot(ode_solve.t,ode_solve.y[0, :]);
plt.plot(ode_solve.t,ode_solve.y[1, :]);
plt.show()

I would like to manually set new values for certain species at pre-defined times. For example I would like in this example to add more of the starting molecule at time t=5 to see how the system evolves. Or, even better, define the values of a variable for all values of t (is my input is e.g. a periodic stimulation)?

Mowgli
  • 101
  • You don't. That would be equivalent to adding a delta pulse to the right side of that component, the ODE solvers in contrast rely on the right side to have a certain degree of smoothness. Simply integrate to 5, change the state and restart the integration for the next segment. – Lutz Lehmann Oct 26 '18 at 21:17
  • @LutzL Right, so I can't change the value "brutally", but what about adding a pre-defined variable (continuous, derivable) parameter? Like a periodic oscillation? – Mowgli Oct 26 '18 at 21:43
  • Like for example `def fun(t, y): return [-0.1 * y[0], +0.1 * y[0]+sin(t)]`? – Lutz Lehmann Oct 26 '18 at 21:47

0 Answers0