Change initial condition (add more mass) with time events in scipy `solve_ivp`

Question

I'm trying to add more mass to the system when the time advanced. I'm trying to obtain the solution for a decay equation: \frac{dy}{dt} = -0.32 * y

The initial condition y(t=0) = 100. I want also at t=6 to add 100 to the solution. I tried doing this in events but without any luck. Any suggestions?

def one(t, y, ke): 
    ydot = -ke * y
    return ydot 

def dose(t, y, ke):
    return y[0] + 100*(t==6)



tspan = [0.0, 12.0]
teval = np.array([0, 2, 4, 6,8, 10, 12])
y0 = [100.0]
ke = 0.32
D = 100
sol = solve_ivp(one, tspan, y0, t_eval=teval, args=(ke,), events=(dose), dense_output=True)
sol.y
plt.plot(sol.t, sol.y[0])

Thanks!

The "state change" event action that you intend to use does not exist in the scipy functions, the only implemented actions are "record" and "terminate". You need to stop the integration, change the state vector, and restart the integration. Then join the results, or plot them piecewise. In https://stackoverflow.com/questions/70376807/desolve-cant-understand I tried to define/explain how events work, as actions performed at roots of functions on the time-state space. The event function in scipy is just the function that is tested for roots. — Lutz Lehmann, Apr 13 '22 at 05:27
Thanks for the suggestions. It does make sense. Any example with stopping the integration at specific time and add more mass to the y vector? Thanks! — Mutaz Jaber, Apr 14 '22 at 06:31
I've posted a possible infrastructure for this approach, even if the setup of the problem is different, in https://stackoverflow.com/questions/54767096/odeint-returns-wrong-results — Lutz Lehmann, Apr 14 '22 at 07:42

score 0 · Answer 1 · answered Apr 12 '22 at 20:53

A way to change the output at a given point is to add an impulse to its derivative, here an example with a gaussian pulse (sigma adjusts the width of the pulse). I also used BDF method instead of the default Runge-Kutta, in order to be able to handle narrower pulses, but the limitation of this approach is that a very narrow pulse will not be noticed by the ODE integration

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
def one(t, y, ke): 
    sigma = 0.1
    ydot = -ke * y + 100/(sigma*np.sqrt(2*np.pi))*np.exp(-(t - 6)**2/(2*sigma**2))
    return ydot 



tspan = [0.0, 12.0]
teval = np.array([0, 2, 4, 6,8, 10, 12])
y0 = [100.0]
ke = 0.32
D = 100
sol = solve_ivp(one, tspan, y0, t_eval=np.linspace(0, 12, 1000),
                args=(ke,), dense_output=True, method='BDF')
sol.y
plt.plot(sol.t, sol.y[0])

Thanks for the suggestion. The downfall for this is the inferences! since the mass I want to add represent drug mass and the width of the pulse should be govern by pharmacokinetic behavior of the drug. — Mutaz Jaber, Apr 14 '22 at 06:35
In that case you may want to model this pharmacokinetic behavior as a second differential equation. — Bob, Apr 14 '22 at 14:46

Change initial condition (add more mass) with time events in scipy `solve_ivp`

1 Answers1