为常微分方程组实现循环答案

【问题标题】：Implementing a loop for a system of ordinary differential equations为常微分方程组实现循环
【发布时间】：2014-04-03 04:07:17
【问题描述】：

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

有人可以帮帮我吗？

【问题讨论】：

能否请您记录下第二种情况的微分方程是什么？您应该返回一个长度为 10 的数组，您在更正版本的 i=0 循环中返回一个长度为 2 的数组。没有明显的修复方法。
我只是想用“n”个方程求解 ODE 系统。我是生物学家。我正在模拟一个具有“n”个种群的系统，每个种群都指一个方程。
是的，但是方程式是什么？在第一种情况下，n=2，两个分量有 2 个方程。但在第二种情况下，当 n=10 时，您尝试组装 10 对函数，总共有 20 个。除了第一个返回只返回一对。如果没有进一步的知识，就无法以有意义的方式修复它。
好的，我将发布一组更简单的方程式，以便将注意力集中在代码问题上。
我去寻找 odeint 的非平凡系统示例，但一无所获。此外，如果需要 (numpy/scipy) 数组对象或只需要裸 [1,2,3] 数组，两者都可以找到，这是不一致的。我在下面的答案中放置了 5 个相同的 Voltera-Lotka 系统，这应该很容易扩展到真正耦合的系统，无论有无循环。

标签： python scipy ode differential-equations

【解决方案1】：

使用 odeint 超过 4 维或复杂函数的文档和示例确实很少。总维度为 10 的 5 个（独立）Voltera-Lotka 系统的工作实现可能如下所示：

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

【讨论】：

【解决方案2】：

这是编辑的代码，你的语法关闭：

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

【讨论】：

真的很有用吗？函数在i = 0处的第一个循环中返回，从未达到i = 1。 span>
它使代码工作没有错误消息（至少是！）。但是，它只解决了一个方程。似乎循环不起作用！什么是错的？ span>