SDE solver in python: manual determination of integrator step size (dt)

Aim: I am trying to solve a system of SDEs, while using the SDEint package in python 3.x. It is a system of SDEs adapted from and inspired by the Zombie Apocalypse model, published by Munz, 2009. I tried to add an explanation of the terms within the code.

Problem: The problem I would like to solve occurs in the example below (overflow encoutered in double_scalar). Note that the noice term is currently 0, so it is basically a system of ODEs and not yet a system of SDEs.

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

p, f, e, k = 0.7, 0.0002, 0.05, 0.02    # human birth rate, confilt occurence rate, conflict end rate, zombie killing rate

tspan = np.linspace(0, 20., 100)
y0 = np.array([100000., 0., 0., p])

def ff(y, t):
Hi = y[0]
Ci = y[1]
Zi = y[2]

# Human are created (exponential)                       --> y[3] * Hi
# Humans and zombies engage in conflict at rate (f)     --> - f * Hi * Zi
# the conficts end at rate (e) without a victor         --> + e * Ci
f0 = y[3] * Hi - f * Hi * Zi + e * Ci

# Humans and zombies engage in conflict at rate (f)     --> f * Hi * Zi
# the conficts end at rate (e) without a victor         --> - e * Ci
# and zombies kill humans in conflict at rate (k)       --> - k * Ci
f1 = f * Hi * Zi - e * Ci - k * Ci

# Zombies arise from the earth at a fixed rate          --> 10**7
# they arise as victor from a conflict                  --> + k * Ci
# enter conflicts                                       --> - f * Hi * Zi
# and exit conflicts without killing                    --> + e * Ci
f2 = 10**7 + k * Ci - f * Hi * Zi + e * Ci

f3 = 0
print(f0, f1 ,f2, f3)
return np.array([f0, f1, f2, f3])

def GG(y, t):
return np.diag([0, 0, 0, 0])

result = sdeint.itoint(ff, GG, y0, tspan)

fig = plt.figure()
ax.set_yscale('log')
ax.set_ylim(bottom=1, top=10**13)
plt.plot(result)
plt.show()


I have tried to solve the same system of ODEs (without noise) using Scipy's odeint and this seems to work well due to the possibility of adding several arguments to odeint.

result = odeint(ff, y0, tspan, mxstep=10 ** 9, rtol=10 ** (-3),
atol=10 ** (-3), hmax=1)


So I thought maybe SDEint uses too large integration steps, however there is no possibility to manually add additional arguments to sdeint.

Potential Solution(?): So I was wondering, if a solution might be to manually set the set size of the integration. I tried to do so by increasing the amount of steps in np.linspace e.g.

tspan = np.linspace(0, 20., 100000) # 100 --> 100000


However, this does not seem to work properly.

Question: Since I am not really familiar with ODEs/SDEs, I am not sure if I am addressing this problem properly. Is this a common way to address the problem of too large integration steps? And might the solution be to increase the step size or am I completely mistaken?

Any help would be greatly appreciated!

• Could you add your reasoning for the parameter choice? You add to a variable with initial value 0.0001 that is constant in the original version a noise of variance 1e4*dt. With such a perturbation I could imagine that any non-linear model rapidly finds a region of positive feed-back and then divergence. – LutzL Feb 8 at 13:14
• Yes, you are completely right. I tried to come up with a reproducible example, however maybe it was not the best... The main problem is that I translated an ODE system into SDEs. odeint was able to integrate perfectly (small step size, args atol/rtol 10^-3), however the same system (even without noise term) explodes in SDEint. I thought that a manually fixed step size might solve the problem, however I am not completely sure. – user213544 Feb 8 at 14:19
• I tried to come up with an improved reproducible example of the problem I am facing. Hopefully I have succeeded:) In short: I am able to solve an ODE system and I would like to make it into an SDE system, so I used the package SDEint. Without the ability to add arguments to sdeint the system crashed, I think due to a larger integration step size. – user213544 Feb 10 at 11:20
• SDEint without noise should reduce to the Heun or explicit trapezoidal method, a second order method. It is possible that variables become accidentially negative due to the larger step errors and then enter a positive feed-back cycle exploding the system. – LutzL Feb 10 at 11:37
• Is there a way to prevent these accidental negative variables or larger step errors? – user213544 Feb 10 at 13:38