Scipy odeint 给出索引越界错误

Question

我正在尝试使用 Scipy 的 odeint 函数求解 Python 中的微分方程。该方程的形式为dy/dt = w(t)，其中w(t) = w1*(1+A*sin(w2*t)) 用于某些参数w1、w2 和A。我编写的代码适用于某些参数，但对于其他参数，我得到的索引超出范围错误。

这是一些有效的示例代码

import numpy as np
import scipy.integrate as integrate

t = np.arange(1000)

w1 = 2*np.pi
w2 = 0.016*np.pi
A = 1.0

w = w1*(1+A*np.sin(w2*t))

def f(y,t0):
    return w[t0]

y = integrate.odeint(f,0,t)

这是一些不起作用的示例代码

import numpy as np
import scipy.integrate as integrate

t = np.arange(1000)

w1 = 0.3*np.pi
w2 = 0.005*np.pi
A = 0.15

w = w1*(1+A*np.sin(w2*t))

def f(y,t0):
    return w[t0]

y = integrate.odeint(f,0,t)

这两者之间唯一不同的是，w1、w2和A这三个参数在第二个中变小了，但是第二个总是给我以下错误

line 13, in f 
   return w[t0]

IndexError: index 1001 is out of bounds for axis 0 with size 1000

即使在重新启动 python 并首先运行第二个代码后，此错误仍然存​​在。我尝试过使用其他参数，有些似乎有效，但其他参数给了我不同的索引越界错误。有人说 1001 超出范围，有人说 1000，有人说 1008，等等。

更改 y 的初始条件（odeint 的第二个输入，我在上面的代码中将其设为 0）也会更改索引错误上的数字，因此可能是我误解了该放在这里的内容。除了将 y 用作信号的相位之外，我没有被告知初始条件应该是什么，所以我假设它最初是 0。

Answer 1

你想做的是

def w(t):
    return w1*(1+A*np.sin(w2*t))

def f(y,t0):
   return w(t0)

数组索引通常是整数，时间参数和微分方程解的值通常是实数。因此在调用w[t0] 时存在一些概念上的困难。

你也可以尝试直接集成函数w，这个例子没有内在的困难。

---

对于耦合系统，您可以将它们作为耦合系统来解决。

def w(t):
    return w1*(1+A*np.sin(w2*t))

def f(y,t):
   wt = w(t)
   return np.array([ wt, wt*sin(y[1]-y[0]) ])