How to solve in python a set of ODEs where one of them contains a dirac delta function (or stepwise increment)?

Question

I'm trying to visualise an oscillator system in python. As long as there is no build-up in one of the variables (with the use of a dirac-delta function), the system works. However, when I try to include it, no change appears to the system. Underneath you can find the process I went through together with some questions I have.

I have written some code in python that is able to correctly visualise the regular oscillator behaviour. However, when including a step-wise increment in one of the variables, the system doesn't react to it. Below, you are able to find the version with the increment included. The code runs but gives the output as if the increment wasn't there. Under the code, you are able to find my problems with getting a correct solution.

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

a4 = 5.68
aAPD = 9.09
deltaZ = 0.41
gamma = 1.33
negV = 1
n=500

def f(V):
    if V<=-60:
        return V/20 + 4
    elif V>=20:
        return V/20 - 1
    else:
        return -V/80 + 1/4

def alpha(V):
    if V<0:
        return a4
    else:
        return -aAPD

def g(q):
    return q/(q+1)

def step(V):
    global negV
    if V>0 and negV == 1:
        negV = 0
        return deltaZ
    elif V<=0 and negV == 0:
        negV = 1
        return 0
    else:
        return 0

# function that returns dz/dt
def model(z,t):
    V = z[0]
    y = z[1]
    q = z[2]
    dVdt = 25000*(y - f(V))
    dydt = alpha(V) - a4*g(q)
    dqdt = -gamma*g(q) + step(V)
    dzdt = [dVdt,dydt,dqdt]
    return dzdt

# initial condition
z0 = [0,0,0]

# time points
t = np.linspace(0,0.5,n)

# solve ODE
z = odeint(model,z0,t)

# plot results
plt.plot(t,z[:,0],label='Voltage')
plt.plot(t,z[:,1],label='y')
plt.plot(t,z[:,2],label='Z')
plt.ylabel('response')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

The problem lies in the function step(V). The oscillator crosses at a certain moment the value V=0. At this timepoint, the value of q should be incremented by deltaZ.

There are several points in the odesolver of python I dn't understand or don't know how to get access to the correct values. Therefore here a short list with problems/questions:

• To find the point where V crosses 0, I make use of a global value negV. However, it would be nicer if I would be able to include a statement along the lines of if V[i]>0 and V[i-1]<0: return deltaZ. This information should be available somewhere since you plot the values of V in the end, but I don't see how I can access previous values of V.
• the increment of the q value should happen instantaneously. Therefore one solution might be to first solve dqdt, and then increment q by the value deltaZ. I don't see how this can be accomplished however. Is there a way to directly access the values of V,y and q? Later on I would like to alter the value of V manually as well, so this would be useful.
• another solution might be to keep the deltaZ value in dqdt but then I should know the timestep with which the differential equations get solved, That way I am able to divide deltaZ=0.41 by this timestep such that the increment of q is deltaZ=0.41. When looking at the odeint information I saw however that the timestep is variable, so perhaps this is not that good of a solution...

Preferably I would get all three points solved, but if you are able to provide the answer to one of them, I would also already be helped a lot.

Answer 1

I highly recommend using the solve_ivp function in scipy using an event to detect the V=0 crossing. See here.

If you use a terminal event, then the solution will stop precisely at the point of V=0. You can then update the solution outside of the ODE integrator with your instantaneous change and use this as the initial condition for a second call to solve_ivp to finish the integration. I assume that the crossing of V=0 only happens once, but you could build in logic to handle multiple crossings as well with events.