My goal is to plot the following set of coupled ODEs:
The odeint method was giving me problems (such as some of the x_i's taking on negative values which analytically does not happen), so I decided to use a fourth-order Runge-Kutta solver. I was following the examples here to use the dopri5 method (Using adaptive step sizes with scipy.integrate.ode):
import scipy as sp
import pylab as plt
import numpy as np
import scipy.integrate as spi
#Constants
c13 = 6.2
c14 = 1.0
c21 = 7.3
c32 = 2.4
c34 = 12.7
c42 = 5.7
c43 = 5.0
e12 = 1.5
e23 = 2.5
e24 = 2.0
e31 = 3.0
e41 = 4.8
#Time
t_end = 700
t_start = 0
t_step = 1
t_interval = sp.arange(t_start, t_end, t_step)
#Initial Condition
ic = [0.2,0.3,0.3,0.5]
def model(t,ic):
Eqs= np.zeros((4))
Eqs[0] = (ic[0]*(1-ic[0]*ic[0]-ic[1]*ic[1]-ic[2]*ic[2]-ic[3]*ic[3])-c21*((ic[1]*ic[1])*ic[0])+e31*((ic[2]*ic[2])*ic[0])+e41*((ic[3]*ic[3])*ic[0]))
Eqs[1] = (ic[1]*(1-ic[0]*ic[0]-ic[1]*ic[1]-ic[2]*ic[2]-ic[3]*ic[3])+e12*((ic[0]*ic[0])*ic[1])-c32*((ic[2]*ic[2])*ic[1])-c42*((ic[3]*ic[3])*ic[1]))
Eqs[2] = (ic[2]*(1-ic[0]*ic[0]-ic[1]*ic[1]-ic[2]*ic[2]-ic[3]*ic[3])-c13*((ic[0]*ic[0])*ic[2])+e23*((ic[1]*ic[1])*ic[2])-c43*((ic[3]*ic[3])*ic[2]))
Eqs[3] = (ic[3]*(1-ic[0]*ic[0]-ic[1]*ic[1]-ic[2]*ic[2]-ic[3]*ic[3])-c14*((ic[0]*ic[0])*ic[3])+e24*((ic[1]*ic[1])*ic[3])-c34*((ic[2]*ic[2])*ic[3]))
return Eqs
ode = spi.ode(model)
ode.set_integrator('dopri5')
ode.set_initial_value(ic,t_start)
ts = []
ys = []
while ode.successful() and ode.t < t_end:
ode.integrate(ode.t + t_step)
ts.append(ode.t)
ys.append(ode.y)
t = np.vstack(ts)
x1,x2,x3,x4 = np.vstack(ys).T
plt.subplot(1, 1, 1)
plt.plot(t, x1, 'r', label = 'x1')
plt.plot(t, x2, 'b', label = 'x2')
plt.plot(t, x3, 'g', label = 'x3')
plt.plot(t, x4, 'purple', label = 'x4')
plt.xlim([0,t_end])
plt.legend()
plt.ylim([-0.2,1.3])
plt.show()
which produces the following plot:
However, my plots have something odd about them -- the x1 value seems to randomly spike above 1. Because (1,0,0,0) is a fixed point of the dynamical system, it doesn't make much sense to me that it spikes above one (there might've been a reason if the spikes had shown some sort of repetitive pattern, but they look random in the plot, so I'm wondering if that has more to do with just the numerical integration rather than the actual dynamics). Changing parameters also changes that spike (some parameter values I tried still have that spiking but it's substantially less). I therefore have 2 questions:
1) What is causing these spikes to appear? My original thought was that it was the "t_step" being large, so I tried playing around with it (which leads me to my second question)
2) I tried changing the step size from 1 to 0.1. However it produced this plot, which actually looks much different than if step size was 1 (which I thought lower step size would produce more accurate plots?) In this plot, it appears that the time each x_i spends near 1 is longer than if the step size was 1 instead of 0.1, and the heights of the spikes are all equal -- how do I know which plot is more accurate of the true system considering the plots have some substantial differences?
The real crux of the problem I'm working on is looking up the interaction between x3 and x4, but I want to make sure the simulations are set up correctly so that I can get accurate results about x3 and x4!
