Scipy odeint breaks on discontinuous dynamics function e.g. xdot = -sgn(x)

Question

I am trying to use Python and scipy.integrate.odeint to simulate the following dynamical system:

But this integration breaks numerically in Python resulting in the following and similar images (usually even worse than this):

Generated using the following in iPython/Jupyter notebook:

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

f = lambda x,t: -np.sign(x)
x0 = 3
ts = np.linspace(0,10,1000)
xs = odeint(f,x0,ts)
plt.plot(ts,xs)
plt.show()

Any advice how to better simulate such a system with discontinuous dynamics?

Edit #1:

Example result when ran with smaller timestep, ts = np.linspace(0,10,1000000), in response to @Hun's answer. This is also an incorrect result according to my expectations.

Just to add that I get different images upon repeating the execution of the code… — Daniel Naftalovich, Apr 08 '16 at 05:16
I am using python 3.5.1, numpy 1.10.4, scipy 0.17.0. It just worked fine when I ran your code exactly as shown. — Hun, Apr 08 '16 at 16:47

score 0 · Answer 1 · answered Apr 07 '16 at 22:44

0

I ran your code exactly as shown and this is what I got.

answered Apr 07 '16 at 22:44

Hun

3,707
2
15
15

Can you run it again several times please and report if you always get this image? – Daniel Naftalovich Apr 08 '16 at 05:23
(This is what I expect it to look like) – Daniel Naftalovich Apr 08 '16 at 05:57
ts = np.linspace(0,10,1000) # make this time step smaller. For example, make it 1e4. – Hun Apr 08 '16 at 14:59
Does not solve the problem. I ran it on `https://try.jupyter.org/` with `Python 2` with the line `ts = np.linspace(0,10,1000000)` and got a related issue (see image added to my original post, because I can't include it within a comment). – Daniel Naftalovich Apr 08 '16 at 16:34

score 0 · Answer 2 · answered Apr 08 '16 at 01:16

0

I also ran it exactly as shown and got:

answered Apr 08 '16 at 01:16

Ari Cooper-Davis

3,374
3
26
43

This is another type of result that I get sometime… – Daniel Naftalovich Apr 08 '16 at 05:17
Do you have any more information for that dynamical system, is it called anything that I may be able to look up? – Ari Cooper-Davis Apr 09 '16 at 09:36

score 0 · Answer 3 · edited May 23 '17 at 12:24

One approach that can work better is by implementing a custom sgn() function that has epsilon-tolerance around zero, instead of expecting precise equality with zero.

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

sgn = lambda x,eps: (x<-eps)*(-1) + (x>=-eps and x<eps)*0 + (x>=eps)*1
f = lambda x,t: -sgn(x,1e-7)
x0 = 3
ts = np.linspace(0,10,1000)
xs = odeint(f,x0,ts)
plt.plot(ts,xs)
plt.show()

Note that an eps parameter of 1e-7 worked fine for me, but smaller than that seemed to crash. It seems from numpy.finfo and demonstrated here that the machine epsilon of numpy.float32 is around 1e-7, while float is around 1e-16.

Scipy odeint breaks on discontinuous dynamics function e.g. xdot = -sgn(x)

3 Answers3