两个多变量耦合 ODE 的系统

Question

我正在尝试使用 scipy 中的 odeint() 解决耦合 ODE 的以下问题。系统看起来像这样：

X'_k = mean(Y_k) + F

Y'_{k,j} = X_k - Y_{k,j}

这是一个有3个X变量的系统，对于每个X变量，还有另外3个Y变量。

根据我从文档中阅读的内容和示例 here and here，我可以将方程组作为列表传递。这就是我在以下示例中尝试的方法：

import numpy as np
from scipy.integrate import odeint

def dZdt(Z, t):

    X = Z[0]
    Y = Z[1]

    F = 4

    d_x = np.zeros(3)
    d_y = np.zeros(3*3).reshape(3,3)

    # Compute the Y values
    for k in range(3):
        for j in range(3):
            d_y[k][j] = X[k] - Y[k][j]

        # X values
        d_x[k] = Y[k].mean() + F

    d = [d_x, d_y]

    return d

# Initial conditions
X0 = np.random.uniform(size=3)
Y0 = np.random.uniform(size = 3*3).reshape(3,3)
Z0 = [X0, Y0]

t = range(20)

Z = odeint(dZdt, Z0, t)

其中 k、j = (1,2,3) 且 Z = [X,Y]

但我恐怕会遇到以下错误：

ValueError: could not broadcast input array from shape (3,3) into shape (3)

我的实际问题更复杂，因为 j 和 k 可以大于 3（它们分别从 1 到 j_max 和 K_max）所以我不能写 12变量一个一个。

我的猜测是，在代码的某处，Y 变量试图填充 X 形状...但不知道在哪里。

知道我做错了什么吗？

Answer 1

您正试图用列表中的两个数组表示未知函数。它必须是一维数组。因此，不是 3 个 X 变量和 9 个 Y 变量，它必须是 12 个变量的平面列表。像这样：

def dZdt(Z, t):

    X = Z[:3]
    Y = Z[3:].reshape(3, 3)
    F = 4

    d_x = np.zeros(3)
    d_y = np.zeros((3, 3))

    # Compute the Y values
    for k in range(3):
        for j in range(3):
            d_y[k, j] = X[k] - Y[k, j]

        # X values
        d_x[k] = Y[k].mean() + F

    d = np.concatenate((d_x.ravel(), d_y.ravel()))

    return d

# Initial conditions
X0 = np.random.uniform(size=3)
Y0 = np.random.uniform(size=(3, 3))
Z0 = np.concatenate((X0.ravel(), Y0.ravel()))

t = range(20)

Z = odeint(dZdt, Z0, t)

NumPy 数组的索引为 Y[k, j]，而不是 Y[k][j]。并且有大量的矢量化机会可以消除 dZdt 计算中的循环。像这样：

def dZdt(Z, t):

    X = Z[:3]
    Y = Z[3:].reshape(3, 3)
    F = 4

    d_y = X[:, None] - Y 
    d_x = Y.mean(axis=1) + F

    d = np.concatenate((d_x.ravel(), d_y.ravel()))

    return d

两个多变量耦合 ODE 的系统

A system of two multivariable coupled ODEs

python

scipy

odeint