python odeint 给出奇怪的结果

Question

我想了解 scipy.odeint 的工作原理，但我遇到了一些问题。例如：

import numpy as np
from scipy.integrate import  odeint
from scipy.integrate import  odeint

def Diffeq(v,t, lam, gam,c, a):
    vdot = v
    for i in range(0,len(v)):
        if i == 0:
            vdot[0] = c[1]*v[1]- lam*a[0]*v[0]
        elif i == (len(v)-1):
            vdot[i] =  lam*a[i-1]*v[i-1] - (lam*a[i]+c[i])*v[i]
        else:
            vdot[i]=  lam*a[i-1]*v[i-1] - (lam*a[i]+c[i])*v[i]+ c[i+1]*v[i+1]
    print vdot
    return vdot



incond=np.array([0]*900)
incond[1] =1
t = np.linspace(0.0, 2, 1000)
ak = [2*i for i in range(0,900)]
lamma =2
gamma =1
c=[i*gamma for i in range(0,900)]

y = odeint(Diffeq, incond, t, args=(lamma,gamma,c,ak) )

此代码应计算以下形式的微分方程组：

其中

xdot_0 = -(a_0 + c_0)*x_0(t) + c_1*x_1(t)
xdot_899 = a_898*x_898(t) -(a_899 + c_899)*x_900(t)

初始条件 x(0) = (0,1,0...0) 当我尝试分析结果时，我注意到我的函数爆炸到 + 无穷大。如果我使用常量 ak、lamma 和 gamma，我可以得到类似

的结果

x_k(t) = [0,-21,21, 0, 0,...,0]

每次操作。因此，我认为我在代码中犯了一些错误，但我看不到哪里。

Answer 1

在Python中，当你执行行

    vdot = v

vdot 不是 v 的副本。这两个名称现在指的是同一个对象。所以当你在函数Diffeq中修改vdot时，你也修改了输入参数。如果你修改了vdot[0]然后尝试使用v[0]，你实际上得到的是修改后的值vdot[0]，所以你的计算是不正确的。

将该行更改为

    vdot = np.empty_like(v)

当我这样做时（我删除了 print 语句，因此函数在合理的时间内完成），odeint returns 成功。这是解决方案组件的图表：

python odeint 给出奇怪的结果

python odeint gives strange results

python

scipy

differential-equations