# Stochastic differential equation system (SDE) : overflow encountered in double-scalars

I'm trying to integrate the following SDE system from Dekker et al. [1]

$$\begin{cases} \frac{dx}{dt}=a_1x^3+a_2x+\phi+\zeta_x\\ \frac{dy}{dt}=b_1z+b_2(\kappa(x)-(y^2+z^2))y+\zeta_y\\ \frac{dz}{dt}=c_1y+b_2(\kappa(x)-(y^2+z^2))z+\zeta_z \end{cases}$$

To do that I'm using the Euler-Maruyama method which is basically a forward-Euler scheme plus a Gaussian noise term where the mean and the standard deviation parameters are picked up from the Table 2 in Dekker et al. article.

We have

• Noise mean : 0
• Noise variance : 0.1
• Integration time : 500
• Time step : 0.5

The parameters $$a_1, a_2,...$$ as well as the coupling term $$\kappa(x)$$, the integration domain, the step time and the initial conditions are all picked from Table 2.

However, running that simulations I get an RuntimeWarning: overflow encountered in double-scalars. Seems like step time would be too big but that's what is recommended in the article.

I'm wondering why I get this error and how I could fix it?

As you can see, after a few iterations, the decimal of x, y and z become exponentially big. It doesn't come from the stochastic term as the print (i.e the array [x, y, z]) in the below capture is the print of what comes out from the forward-euler alone.

## Code

# 1D EDO for Fold bifurcation paramaters
a1 = -1
a2 = 1

# 2D EDO for Hopf bifurcation parameters
b1 = b2 = 1
c1 = -1
c2 = 1

# coupling parameters
gamma1 = -0.1
gamma2 = 0.12

# initial conditions
x0 = -0.5
y0 = 1.0
z0 = -1.0

r0 = 5.0

# time range
t_init = 0
t_fin = 500
time_step = 0.01

# gaussian noise parameters
mean = 0.0
variance = 0.1

# parameters for stochastic plot
a1_stoch = -1
a2_stoch = 1
b1_stoch = 0.1
b2_stoch = 1
c1_stoch = -0.5
c2_stoch = 1
gamma1_stoch = -0.2
gamma2_stoch = 0.3
time_step_stoch = 0.5

import numpy as np
from consts import (a1_stoch, a2_stoch, b1_stoch, b2_stoch, c1_stoch, c2_stoch, gamma1_stoch, gamma2_stoch)

# linear coupling parameter
# proposed by Dekker et al. article
def gamma(x):
return gamma1_stoch + (gamma2_stoch * x)

def fold_hopf_stoch(v, phi):
print(v, phi)
return np.array([
a1_stoch * (v[0] ** 3) + a2_stoch * v[0] + phi,
b1_stoch * v[2] + b2_stoch * (gamma(v[0]) - (v[1]**2 + v[2]**2))*v[1],
c1_stoch * v[1] + c2_stoch * (gamma(v[0]) - (v[1]**2 + v[2]**2))*v[2]
])


import numpy as np
from consts import (mean, variance)

"""
forward euler method with a stochastical term
x_{i+1} = x_i * dt + zeta + sqrt(dt)
(forward-euler + gaussian noise * sqrt(dt))
"""
def forward_euler_maruyama(edo, v, dt, *args):
zeta = np.random.normal(loc=mean, scale=np.sqrt(variance))
print(zeta, edo(v, *args) * dt)
return edo(v, *args) * dt + zeta * np.sqrt(dt)


import numpy as np
import matplotlib.pyplot as plt

from tqdm import tqdm

from .edo import (fold_hopf)
from .edo_stoch import (fold_hopf_stoch)
from .rk4 import rk4
from .euler import forward_euler_maruyama

from consts import (x0, y0, z0, t_init, t_fin, time_step, time_step_stoch)

plt.style.use("seaborn-whitegrid")
np.set_printoptions(precision=3, suppress=True)

class time_series():
def __init__(self):
# initial conditions
self.initial_conditions = np.array([[x0, y0, z0]])
# forcing parameter
self.phi = -2
# time
self.t0 = t_init
self.tN = t_fin
# stochastic number of iterations
self.niter = 10
# legend
self.legends = ["$$x$$ (leading)", "$$y$$ (following)", "$$z$$ (following)"]

def solve(self, solver, edo, dt, nt):
# [x0, y0, z0]
#v = self.initial_conditions[0]

# vector mesh -- will receive the result
v_mesh = np.ones((nt, 3)) #
# set inital conditions
v_mesh[0] = self.initial_conditions[0]
print(edo, dt)

for t in tqdm(range(0, nt - 1)):
v_mesh[t + 1] = v_mesh[t] + solver(edo, v_mesh[t], dt, self.phi)

# increase forcing parameter
self.phi += 0.001

return v_mesh

def basic(self):
dt = time_step
nt = int((self.tN - self.t0) / dt)
time_mesh_basic = np.linspace(start=self.t0, stop=self.tN, num=nt)

# [ [x0, y0, z0], [x1, y1, z1], ..., [xN, yN, zN] ]
results = self.solve(rk4, fold_hopf, dt, nt)
return time_mesh_basic, results

def stochastic(self):
dt = time_step_stoch
nt = int((self.tN - self.t0) / dt)
time_mesh_stoch = np.linspace(start=self.t0, stop=self.tN, num=nt)

stochastic_results = np.ones((self.niter, nt, 3))
# compute a lot of simulations
for i in tqdm(range(0, self.niter)):
results = self.solve(forward_euler_maruyama, fold_hopf_stoch, dt, nt)
stochastic_results[i] = results
return
return time_mesh_stoch, stochastic_results

def plot(self):
#time_mesh_basic, basic_results = self.basic()
time_mesh_stoch, stochastic_results = self.stochastic()
return

# 3 lines, 1 column
fig, ((ax1), (ax2)) = plt.subplots(2, 1, figsize=(15, 7))
fig.suptitle("Série temporelle")

ax1.plot(time_mesh_basic, basic_results[:,0])
ax1.plot(time_mesh_basic, basic_results[:,1])
ax1.plot(time_mesh_basic, basic_results[:,2])

ax1.set_xlabel("$$t$$")
ax1.set_ylabel("$$x$$, $$y$$, $$z$$")
ax1.set_xlim(0,500)
ax1.set_ylim(-1.5,1.5)
ax1.legend(self.legends, loc="upper right")
ax1.set_title("Basique")

print(stochastic_results)

ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,0])
ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,1])
ax2.stackplot(time_mesh_stoch, stochastic_results[:,:,2])

ax2.set_xlabel("$$t$$")
ax2.set_ylabel("$$x$$, $$y$$, $$z$$")
ax2.set_xlim(0,500)
ax2.set_ylim(-1.5,1.5)
ax2.legend(self.legends, loc="center left", bbox_to_anchor=(1,0.5))
ax2.set_title("Stochastique")

plt.tight_layout()
plt.show()



### References:

1. Dekker, M. M.; von Der Heydt, A. S.; Dijkstra, H. A. Cascading transitions in the climate system. Earth Syst. Dyn. 2018, 9 (4), 1243–1260. DOI: 10.5194/esd-9-1243-2018.
• (It is slightly irritating that the kappa in the formulas becomes gamma in the code.) By Lipschitz or stability reasoning you need $(3x^2-1)h<2$ for the Euler method not to explode. For $h=0.5$ this is the case for about $|x|<1.3$. Whenever the integration goes over that limit, the sequence for $x$ will explode, rapidly. Mar 15, 2021 at 20:33
• (Ok, the kappa is a gamma in all the other formula sets.) I do not think that in that paper the perturbed equations were considered in the sense of Ito SDE. As I read it, the integration interval is divided into segments of length 0.5, and in each segment the random term is constant. Then apply RK4 or or whatever with an adapted step size to integrate over each segment. Mar 15, 2021 at 21:22
• The paper doesn't seem to give many details about their discretization scheme and there's no supplementary material. It's possible they used an implicit method for the deterministic part of the ODE. Probably the best thing to do is email the corresponding author and ask for the scripts they used; the software and data availability section at the end says it's available upon request and in any case they'll probably be thrilled that someone's interested enough in their work to try it out. Mar 15, 2021 at 22:12