假设我想解一列方程，但它们并不耦合，只是由numpy数组定义的不同。
第一个月
其中c只是numpy数组中给定的一个已知常数，例如c = np.array([1, 2])。这意味着我需要求解两个方程，可能还需要不同的初始值。当然，这可以很容易地通过解析求解。
y(t) = t^2 + c*t + C个
其中C是另一个依赖于y(0)的初始值的常数。例如，第一个方程有C = y(0) = 0，第二个方程有C = y(0) = 1。然后我可以得到精确解：

y =  t^2 + 1*t
y =  t^2 + 2*t + 1

字符串
如果我想使用scipy.integrate.solve_ivp来获取t=1处的值，那么我可以“假装”它们是耦合的。

def test_fun(t, y):
    c = np.array([1, 2])
    dy1 = 2*t + c
    return dy1
sol = solve_ivp(test_fun, [0, 1], np.array([0, 1]))
y = sol.y[:, -1]

型
这实际上是给我的y = np.array([2., 4.])，它与给予的解析值相匹配。
但是现在假设我有一个带有常数的方程列表，不同的是一个numpy数组，但是每个方程都由3个耦合的微分方程组成。让我们仍然使用一个玩具示例，它们不是真正的“耦合”，而是只需要在函数中返回多个值，就像你在求解耦合微分方程时所做的那样。

def test_fun(t, y):
    c = np.array([1, 2])
    dy1 = 2*t + c
    dy2 = 3*t + c
    dy3 = 4*t + c
    return dy1, dy2, dy3

sol = solve_ivp(test_fun, [0, 1], np.array([[1, 2], [2, 2], [3, 2]]))

型
但是现在它会抱怨ValueError: y0 must be 1-dimensional. Flatten y0似乎没有意义，因为我希望函数的输入是一个大小为2的数组，但是对于每个元素，它将返回一个由3个耦合方程组成的系统。如果我只是flatten，它将假设它是一个由6个耦合方程组成的系统的输入，然后抱怨ValueError: could not broadcast input array from shape (3,2) into shape (6,)
当然，我可以只使用for循环，但一般的想法是避免它，因为它在python中非常慢，而且我的c数组非常大。