Attempting to perturb ODE when initial condition is equilibrium point does not work

Question

I have the following system of differential equations:

$$ x' = ax- cy + e1 $$

$$y' = by- dx + e2 $$

for variables $x,y$ and parameters $a,b,c,d,e1,e2$.

I'd like to solve this in python, which is easy enough with odeint.

def f(y, t, params):
    weaponsC1, weaponsC2 = y
    a, b, c, d, e1, e2 = params  
    t_0 = 10
    t_1 = 11
    derivs = [a*x- c*y + e1,       
              b*y - d*x + e2]
    return derivs


#solve
psoln = odeint(f, y0, t, args=(params,), rtol=1e-10)

This returns what you'd expect:

For

$$a = 1\\ b = 1 \\ c = 2 \\ d = 2 \\ e1 = 5 \\ e2 = 5$$

and initial condition: $(2,2)$, the code above returns:

Now I'd like to add a forcing function $g$ such that I can perturb the ODE by some constant amount at some time. After a bit of a struggle with neater methods, I settled on doing it like this:

def f(y, t, params):
    weaponsC1, weaponsC2 = y
    a, b, c, d, e1, e2 = params  
    t_0 = 10
    t_1 = 11
    derivs = [a*weaponsC1- c*weaponsC2 + e1, + (2*(t > t_0)*(t < t_0+0.5)) + (5*(t > t_1)*(t < t_1+0.5)),       
              b*weaponsC2 - d*weaponsC1 + e2 + (2*(t > t_0)*(t < t_0+0.5))+ (5*(t > t_1)*(t < t_1+0.5))]
    return derivs

A graph of just the forcing function is shown, thought at different time value:

def g(t_0):
    return((2*(t > t_0)*(t < t_0+0.5)))

plt.plot(t,(g(3) + g(4)))

when I plot the ODE with this forcing function, it works fine for any initial condition which is not the equilibrium:

However, when I set the equilibrium as the starting condition, try as I may, I can never impact the ODE:

this is the same code that made the above working solution; the only change is the starting condition.

I know ODEint assumes the diffeq is differentiable, and mine is not — could this be the cause of the issue? Or is this the result of some property of the ODE?

Answer 1

It is possible that since you are exactly at equilibrium, your integrator is taking step sizes too large to resolve your forcing function. This is especially bad for your problem because the forcing function is discontinuous with little support, so the chances that the integrator time steps line up with the discontinuities is pretty low.

This is typically handled in integrators either by events or by specifying time points where the solver is forced to resolve the solution to the desired tolerance with small step sizes. odeint does not handle events, but possibly you could set the "tcrit" parameter to be the discontinuity points of your forcing function, which will force the solver to take extra care around those points