# SciPy odeint fails in unpredictable ways on deterministic system of ODEs

I've been trying to solve the following (relatively simple) system of Lotka-Volterra ODEs in Python using SciPy's odeint:

$$\dot{z_1} = z_1 \left(- \sigma z_1 + \sigma z_2 + \rho z_3 - z_4 - z_5\right)$$ $$\dot{z_3} = z_3 \left(\sigma z_1 - \sigma z_2 - \rho z_3 + z_4 + z_5\right)$$ $$\dot{z_4} = z_4 \left(\sigma z_1 - \sigma z_2 - \rho z_3 + z_4 + z_5 + z_6 - \beta z_7\right)$$ $$\dot{z_6} = z_6 \left(\sigma z_1 - \sigma z_2 + \rho z_3 - z_4 - z_5 - z_6 + \beta z_7\right)$$ $$\dot{z_2}=\dot{z_5}=\dot{z_7}=0$$

Where $$z_i$$ for $$i \in [1,7]$$ are the variables I'm solving for and $$\rho,\sigma,$$ and $$\beta$$ are positive constants.

However, whenever I run the code something seems to blow up and the results make no sense (seemingly independent of initial conditions). To make matters even more odd, the way that the integrator seems to fail changes with each run. As an example of how non-sensical these failures are, below this I'm including a list of values being output by odeint as the first 100 timesteps of the "solution" to $$z_2$$ (note that $$\dot{z_2} = 0$$ and so it should never move from its initial value of $$1$$):

[1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000 0.00000000e+000 nan 8.16878018e-302 4.63067615e-273 1.57791392e+161 1.88927554e-299 4.17015339e+098 2.78103105e+180 0.00000000e+000 1.31355181e-308 5.12621759e+199 2.64429328e+180 7.69843740e+218 2.64522460e+185 1.42853987e+248 6.01346953e-154 1.22187976e+224 8.78421503e+247 9.18149096e+170 8.49025266e+175 1.94862288e-153 1.39806877e-152 7.17236781e+252 2.28505331e+170 2.63114242e+207 3.68015657e+180 1.03564511e-308 2.21688273e-313 6.46991602e-067 9.09136099e+276 1.45858877e-303 3.85693955e-310 7.93726941e-301 5.53200078e+199 6.01347002e-154 6.01347002e-154 9.42162457e+221 3.06567899e-313 6.35538057e+149 1.91084646e+214 1.59186898e-308 2.32999935e+031 1.95281786e-258 2.42519188e+214 4.62444261e-273 1.45544717e-291 5.84028529e-302 1.89988452e-308 1.28625507e+248 1.45913393e-152 1.71870550e+161 5.96083817e+175 5.55603592e+180 2.19980295e-152 1.71870549e+161]

And again, a list of the values being output by odeint as the first 100 timesteps of the "solution" to $$z_2$$ after re-running the program with all else the same:

[ 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.00000000e+000 1.15846863e+171 6.17113003e-315 4.03901783e-270 3.79160818e-302 4.41534453e+199 2.42766858e-154 6.01334434e-154 1.12947293e-094 6.96407958e+252 8.76427658e+252 2.59338567e+161 1.69599801e+161 5.36835912e+252 1.69288051e+190 2.61191759e+209 2.43812975e-152 1.06263900e+248 9.49690810e-154 2.44085720e-154 2.31633990e-152 1.75533829e+156 5.16596794e-109 8.80088337e+199 1.48889777e+195 1.21697995e-152 8.80088337e+199 2.08063217e-115 4.47593816e-091 3.15474910e+180 1.03417828e-028 1.46921998e+195 3.81388253e+180 1.23478424e-259 2.31633990e-152 1.19461721e+190 6.01334671e-154 1.21698002e-152 6.01334512e-154 3.99727850e+252 4.56969612e+233 1.94862175e-153 1.27873527e-152 1.27990068e-152 6.01334435e-154 2.65141193e+180 1.54069978e+228 6.01334512e-154 1.10901207e+200 1.81590214e-152 2.63114487e+207 5.31258342e-193 4.40505856e-193 1.96192454e+289 -1.79811598e-284 -3.32653116e-111 -3.38460706e+125 4.04547221e-311 4.16867019e-290]

Note that on both runs odeint seems to run without a hitch for a while, but ultimately changes $$z_2$$ and then oscillates wildly. The other variables all behave in a similarly erratic way and what's weird is that none of the values are particularly large/small on the time steps before failure.

I'm really not sure what is going on (nor how to fix it). The root of my problem seems to be different from other things I've seen online (e.g. I'm not dividing by zero or using some misbehaved/non-smooth function) and the solutions proposed to those problems were usually not applicable (or didn't work). Any advise, intuitions, or solutions is appreciated. I'll include my actual code down below this so that y'all can have a look. Also, for what it is worth, I have also tried implementing this system in Mathematica and got errors there too.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.interpolate import LinearNDInterpolator as interpolate
from mpl_toolkits.mplot3d import Axes3D

#Defining some constants and a matrix for my Lotka-Volterra system
rho = 28.0
sigma = 10.0
beta = 8.0 / 3.0

lv_mat = np.array([[-sigma,sigma,rho,-1.0,-1.0,0,0],[0,0,0,0,0,0,0],[sigma,-sigma,-rho,1.0,1.0,0,0],[sigma,-sigma,-rho,1.0,1.0,1.0,-beta],[0,0,0,0,0,0,0],[sigma,-sigma,rho,-1.0,-1.0,-1.0,beta],[0,0,0,0,0,0,0]])

def lv_lorenz(state, t):
z0, z1, z2, z3, z4, z5, z6 = state  # Unpack the state vector

fitness = lv_mat.dot(state) #Compute term for derivatives

#Compute all derivatives
z0_dot = z0 * fitness[0]
z1_dot = z1 * fitness[1]
z2_dot = z2 * fitness[2]
z3_dot = z3 * fitness[3]
z4_dot = z4 * fitness[4]
z5_dot = z5 * fitness[5]
z6_dot = z6 * fitness[6]

return z0_dot, z1_dot, z2_dot, z3_dot, z4_dot, z5_dot, z6_dot

if __name__ == '__main__':
init_state = np.array([1.,1.,1.,1.,1.,1.,1.]) #Initial state

t = np.arange(0.0, 1.0, 0.01) #Range of time to integrate over
states = odeint(lv_lorenz, init_state, t) #Call odeint. I also tried setting mxstep to multiple different values

#Plot the "solution" to this system
plt.plot(t, states[:, 0], label='$$z_1$$ values')
plt.plot(t, states[:, 1], label='$$z_2$$ values')
plt.plot(t, states[:, 2], label='$$z_3$$ values')
plt.plot(t, states[:, 3], label='$$z_4$$ values')
plt.plot(t, states[:, 4], label='$$z_5$$ values')
plt.plot(t, states[:, 5], label='$$z_6$$ values')
plt.plot(t, states[:, 6], label='$$z_7$$ values')
plt.xlabel('Time')
plt.ylabel('Solution')
plt.ylim(-1,1000)
plt.legend()
plt.show()


Thanks in advance! (And, if anything, I hope it's at least an interesting quasi-linear system that seems to break integrators.)

• Do you have reason to believe this system should be stable for your choice of coefficients? You can compute the eigenvalues to check. – Spencer Bryngelson Oct 31 '20 at 4:34
• Have you tried solve_ivp using different solver methods? – fibonatic Oct 31 '20 at 7:40
• You can shorten the ODE function to return state*lv_mat.dot(state), array multiplication like any other "scalar" array operations are element-wise. – Lutz Lehmann Oct 31 '20 at 13:41
• @SpencerBryngelson I'm not totally sure in what sense you mean stable, but this system should be diffeomorphic to the Lorenz system. So I expect it to be as stable as chaotic systems get I suppose. I expected it to be stable enough for the integrator since 1) directly computing the solution from one of Lorenz's solutions/trajectories works nicely enough (i.e. using the diffeomorphism's map from state space to state space), and 2) the straightforward quadratic form of the individual ODEs looked harmless enough (to my non-numerics savvy eyes, which I realize isn't a particularly good reason). – Dupin Oct 31 '20 at 16:45
• @fibonatic Other than also trying it in Mathematica + toying arguments to the solver, no I haven't yet. Do you have any suggestions? I'll try taking a stab with other solvers just to see what happens (certainly can't hurt). – Dupin Oct 31 '20 at 16:46

Shortly after reaching $$t=0.4$$ the system has a pole, in component $$z_4$$ with component $$z_3$$ following. Such behavior is to be expected for systems with positive quadratic terms, similar to $$y'=y^2$$. This leads to overflow and the interruption of the integration. The remaining array elements are/remain filled with the random garbage that they contained at the allocation of the memory segment to the array, thus these can be different in different runs.
The error message 'Excess work done on this call (perhaps wrong Dfun type).' is misleading, I would expect that the step size controller falls below the minimal feasible step size at the pole. There is an error message of this type that also displays the last valid time, i have no idea why this does not occur here.
• If you use full_output=True, then in the info structure there is an array tcur where one can see what the last sensible integration time was. // If you expect the solution to exist for a longer time, then the transformation might contain an error, a flipped sign or something. – Lutz Lehmann Oct 31 '20 at 17:06
• The "constants of integration" are the initial conditions, or in general there is a (locally) bijective relationship. Could it be that the constants $z_2,z_5,z_7$ should have larger values, like the shift value of $100$? – Lutz Lehmann Oct 31 '20 at 19:18