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 Ci = y Zi = y # Human are created (exponential) --> y * 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 * 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 = fig.add_subplot(1,1,1) 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
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
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
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!