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!

Answer 1

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])