0
$\begingroup$

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.)

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ Do you have reason to believe this system should be stable for your choice of coefficients? You can compute the eigenvalues to check. $\endgroup$
    – user20857
    Oct 31, 2020 at 4:34
  • $\begingroup$ Have you tried solve_ivp using different solver methods? $\endgroup$
    – fibonatic
    Oct 31, 2020 at 7:40
  • $\begingroup$ You can shorten the ODE function to return state*lv_mat.dot(state), array multiplication like any other "scalar" array operations are element-wise. $\endgroup$
    – Lutz Lehmann
    Oct 31, 2020 at 13:41
  • $\begingroup$ @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). $\endgroup$
    – Dupin
    Oct 31, 2020 at 16:45
  • $\begingroup$ @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). $\endgroup$
    – Dupin
    Oct 31, 2020 at 16:46

1 Answer 1

Reset to default
3
$\begingroup$

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.

relevant part of the plot

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.

Improve this answer
$\endgroup$
5
  • $\begingroup$ That makes sense, thanks! Do you have any suggestions about what I can do to solve this problem? Is there a more appropriate integrator? If nothing comes to mind for you, do you happen to have a good source for me to read about what sorts of problems this pole is giving the integrator? (I'm currently googling around, but not finding too much practical/actionable information.) Also, as a side note, glad to hear I wasn't losing my mind/being dumb by not knowing how to solve this with that error's suggestion and also thinking that the erratic behaviour looked "almost like an int overflow". $\endgroup$
    – Dupin
    Oct 31, 2020 at 16:57
  • $\begingroup$ 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. $\endgroup$
    – Lutz Lehmann
    Oct 31, 2020 at 17:06
  • 1
    $\begingroup$ Thanks again! After toying around a bit longer I realized that the problem might actually be quite a silly mistake on my end. Namely, this system was derived from a chain of diffeomorphisms that start with the classic Lorenz equations. However, the first step in this chain requires simply "shifting" trajectories of the Lorenz system by a constant like 100 (enough to guarantee the strange attractor is all positive by simply using the diffeomorphic map f(x)=x+C). Obviously this map doesn't change the derivatives and, when I transform the derivatives from Lorenz, this shift isn't accounted for. $\endgroup$
    – Dupin
    Oct 31, 2020 at 18:51
  • $\begingroup$ With all that said, is there a simple way to add a "constant of integration" to python's integrator? It doesn't seem to be in the documentation. (I believe that initial shift equates to simply changing the constant of integration in the original Lorenz system, so hopefully this will do the trick.) $\endgroup$
    – Dupin
    Oct 31, 2020 at 18:52
  • 1
    $\begingroup$ 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$? $\endgroup$
    – Lutz Lehmann
    Oct 31, 2020 at 19:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.