在 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 值。幸运的是，我有一个验证 table 来验证结果。不幸的是，结果并不相同。

这是我的验证table:

这是 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

该方程是 "coupled" 常微分方程（参见 Wikipedia 上的 "System of ODEs"。

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

我修改了您的代码以实现此结果：

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, '-')

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

implement an integration math equation using odeint in Python

python

math

integration

numpy

scipy