Unexpected behavior with Scipy.integrate.odeint on discontinuous function

user_15400465 0
Viewed 118

I am trying to solve the following ordinary differential equation: f'(t) = -2 if 300 <= t <= 303 else 0

import scipy.integrate as integr
import matplotlib.pyplot as plt
Y0=25
def f(Y,t):
    a= -2 if 300 <=t <= 303 else 0
    return a
    
    
T=np.linspace(0,500,5000)
sol=integr.odeint(f,Y0,T)

plt.plot(T,sol)
plt.show()

However the result is only a flat line :
[1]: https://i.stack.imgur.com/KAK4F.png

Whereas it works fine if the interval is bigger : 150 <= t <= 350 instead of 300 <= t <= 303

Any idea why?
Thanks in advance

scipy numpy python
1 Answers

odeint (which is a wrapper to LSODA) or any ODE solvers can't really deal with discontinuities like this in one fellow swoop. ODE solvers normally assume a smooth solution. The solution is flat because odeint is likely taking such large timesteps that it is stepping past the t = 300 to 303 region.

In these situations it is best to do a single integration for each smooth part. So integrate from t = 0 to 300, then stop. Use solution at 300 as initial conditions to integrate from 300 to 303, then stop. Then use solution at t = 303 to integrate from 303 to 500.

There are some packages which can automate this... I know Assimulo has this feature. But I don't know how to use it.

Here is a clunky solution with odeint

import scipy.integrate as integr
import matplotlib.pyplot as plt

def f1(y,t):
    return 0.0

def f2(y,t):
    return -2.0
    
y0 = np.array([25])    
t1=np.linspace(0,300,1000)
sol1=integr.odeint(f1,y0,t1)

t2=np.linspace(300,303,1000)
sol2=integr.odeint(f2,sol1[-1],t2)

t3=np.linspace(303,500,1000)
sol3=integr.odeint(f1,sol2[-1],t3)

sol = np.concatenate((sol1, sol2, sol3))
t = np.concatenate((t1, t2, t3))

plt.plot(t,sol)
plt.show()
user_118153070
Related