Implementing a loop for a system of ordinary differential equations

Question

I am new here and with Python. I am working on a code for numerical solutions of ordinary differential equation systems. My problem: I need to implement a loop inside a function. After that, I need to integrate that function (odeint), but it doesn't work.

Here is the (working) code without the loop:

from scipy.integrate import odeint
from pylab import * # for plotting commands
def deriv(y,t): # return derivatives of the array y
    a = 2.0
    return array([ y*a*(1-y/10)]) 
time = linspace(0,500,1000)
yinit = 1 
y = odeint(deriv,yinit,time)
figure()
plot(time,y)
xlabel('t')
ylabel('y')
show()

And here is the non working code:

from scipy.integrate import odeint
from pylab import * # for plotting commands
def deriv(y,t): # return derivatives of the array y
    a = 2.0
    for i in range(0,10):
        return array([ y[i]*a*(1-y[i]/10) ]) 
time = linspace(0,500,1000)
yinit = linspace(1,1,10) 
y = odeint(deriv,yinit,time)
figure()
plot(time,y[:,0],time,y[:,1])
xlabel('t')
ylabel('y')
show()

Can Someone help me?

Answer 1

The documentation and examples for using odeint with more than 4 dimensions or complicated functions are really rather sparse. A working implementation for 5 (independent) Voltera-Lotka systems with a total dimension of 10 could look like this:

from scipy import zeros_like
from scipy.integrate import odeint
from pylab import * # for plotting commands

def deriv(y,t): # return derivatives of the array y
    a = 2.0
    b = 0.5 #0.1
    c = 0.1
    doty = zeros_like(y)
    for i in range(0,5):
        j=2*i   # this has no deep reason
        k=2*i+1 # just keeps notation short
        doty[j] = y[j]*a*(1-y[j]/10)-b*y[j]*y[k]
        doty[k] = b*y[j]*y[k]-c*y[k]

    return doty

time = linspace(0,500,1000)
yinit = linspace(1,1,10) 
y = odeint(deriv,yinit,time)
figure()
plot(time,y[:,0],time,y[:,1],time,y[:,5],time,y[:,8])
xlabel('t')
ylabel('y')
show()

Answer 2

Here is the edited code, your syntaxing was off:

from scipy.integrate import odeint
from pylab import * # for plotting commands
def deriv(y,t): # return derivatives of the array y
    a = 2.0
    b = 0.5 #0.1
    c = 0.1
    for i in range(0,10): #INSERT A COLON HERE
        return array([ y[i]*a*(1-y[i]/10)-b*y[i]*y[i], b*y[i]*y[i]-c*y[i] ]) #INDENT THIS
time = linspace(0,500,1000)
yinit = linspace(1,1,10)
y = odeint(deriv,yinit,time)
figure()
plot(time,y[:,0],time,y[:,1])
xlabel('t')
ylabel('y')
show()

Your error was here:

for i in range(0,10)
return array([ y[i]*a*(1-y[i]/10)-b*y[i]*y[i], b*y[i]*y[i]-c*y[i] ])

You are supposed to put a colon after the for statement, and indent the return:

for i in range(0,10): #INSERT A COLON HERE
    return array([ y[i]*a*(1-y[i]/10)-b*y[i]*y[i], b*y[i]*y[i]-c*y[i] ]) #INDENT THIS