在 Python 中使用 odeint 实现积分数学方程

Question

我正在尝试使用 scipy.odeint 函数在 python 中求解以下方程。

目前我能够实现这种形式的方程

在 python 中使用以下脚本：

def dY(y1, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 1
    dC = C/N
    b1 = 0
    return (a/dC)*(yin-y1)+b1*dC

x = np.linspace(0,20,1000)
y0 = 0
res = odeint(dY, y0, x)
plt.plot(t,res, '-')
plt.show()

我对第一个等式的问题是“i”。我不知道如何整合方程并且仍然能够提供当前和以前的 'y'(yi-1 和 yi) 值。 'i' 只是一个 0..100 范围内的序列号。

编辑 1：

原方程为：

我用 y,x,a,b 和 C 重写了它

编辑2： 我编辑了 Pierre de Buyl 的代码并更改了 N 值。幸运的是，我有一个验证表来验证结果。不幸的是，结果并不相同。

这是我的验证表：

这里是 numpy 输出：

使用代码：

def dY(y, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 3
    dC = C/N
    b1 = 0.01
    y_diff = -np.copy(y)
    y_diff[0] += yin
    y_diff[1:] += y[:-1]
    return (a/dC)*(y_diff)+b1*dC

x = np.linspace(0,20,11)
y0 = np.zeros(3)
res = odeint(dY, y0, x)
plt.plot(x,res, '-')

如您所见，这些值相差 0.02.. 的偏移量。

我是否遗漏了导致此偏移的某些内容？

Answer 1

该方程是一个“耦合”常微分方程（请参阅Wikipedia 上的“ODE 系统”。

该变量是一个包含y[0]、y[1] 等的向量。要求解 ODE，您必须输入一个向量作为初始条件，并且函数 dY 也必须返回一个向量。

我已经修改了你的代码来实现这个结果：

def dY(y, x):
    a = 0.001
    yin = 1
    C = 0.01
    N = 1
    dC = C/N
    b1 = 0
    y_diff = -np.copy(y)
    y_diff[0] += yin
    y_diff[1:] += y[:-1]
    return (a/dC)*y_diff+b1*dC

我已将y[i-1] - y[i] 部分编写为 NumPy 向量运算，并将坐标 y[0] 特例化（即您的符号中的 y1，但在 Python 中数组从 0 开始）。

解决方案，对所有 yi 使用初始值为 0 是

x = np.linspace(0,20,1000)
y0 = np.zeros(4)
res = odeint(dY, y0, x)
plt.plot(x,res, '-')