为系统的一般微分方程实现循环

Question

我刚来这里，对Python也不太熟悉。我正在写一段代码，用来解决普通微分方程组的数值解。

我遇到的问题是：我需要在一个函数里面实现一个循环。然后，我想用odeint来对这个函数进行积分，但它不工作。

下面是没有循环的（可以正常运行的）代码：

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

这是不工作的代码：

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

有人能帮我吗？

Answer 1

这是修改过的代码，你的语法写错了：

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

你的错误在这里：

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

你应该在 for 语句后面加一个冒号，并且把 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

Answer 2

关于如何使用odeint处理超过4个维度或者复杂函数的文档和示例其实很少。下面是一个可以工作的例子，展示了如何实现5个（独立的）Voltera-Lotka系统，总维度为10：

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