在python中使用for循环求解耦合微分方程

Question

我正在尝试求解一组微分方程，但我一直难以完成这项工作。我的微分方程包含一个“i”下标，表示从 1 到 n 的数字。我尝试按如下方式实现forloop，但我一直收到此索引错误（错误消息如下）。我尝试更改初始条件 (y0) 和其他值，但似乎没有任何效果。在这段代码中，我使用的是solve_ivp。代码如下：

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import solve_ivp

def testmodel(t, y):
    X = y[0]
    Y = y[1]
    J = y[2]
    Q = y[3]
    a = 3
    S = 0.4
    K = 0.8
    L = 2.3
    n = 100 
    for i in range(1,n+1):
        dXdt[i] = K**a+(Q[i]**a) - S*X[i]
        dYdt[i] = (K*X[i])-(L*Y[i])
        dJdt[i] = S*Y[i]-(K*Q[i])
        dQdt[i] = K*X[i]/L+J[i]
        return dXdt, dYdt, dJdt, dQdt
t_span= np.array([0, 120])
times = np.linspace(t_span[0], t_span[1], 1000) 
y0 = 0,0,0,0
soln = solve_ivp(testmodel, t_span, y0, t_eval=times, 
vectorized=True)
t = soln.t
X = soln.y[0]
Y = soln.y[1]
J = soln.y[2]
Q = soln.y[3]

plt.plot(t, X,linewidth=2, color='red')
plt.show()    

我得到的错误是

IndexError                Traceback (most recent call last)
<ipython-input-107-3a0cfa6e42ed> in testmodel(t, y)
     15     n = 100
     16     for i in range(1,n+1):
 --> 17     dXdt[i] = K**a+(Q[i]**a) - S*X[i]
   
 IndexError: index 1 is out of bounds for axis 0 with size 1

我已经分散了网络来解决这个问题，但我无法对这个问题应用任何解决方案。我不确定我做错了什么以及实际要改变什么。

我试图删除“vectorized=True”参数，但随后我收到一个错误，指出我无法索引标量变量。这很令人困惑，因为我认为这些值不应该是标量的。我该如何解决这个问题，我的最终目标是绘制这些微分方程。提前谢谢你。

Answer 1

很高兴您为标准求解器提供了用于多点评估的矢量化 ODE 函数。但默认方法是显式 RK45，显式方法不使用 Jacobi 矩阵。所以不需要对偏导数的差商进行多点求值。

本质上，坐标数组的大小始终为 1，因为评估是在一个点上进行的，所以例如 Q 是一个长度为 1 的数组，唯一有效的索引是 0。记住，在所有“真”编程语言，数组索引从 0 开始。只有一些 CAS 脚本语言使用“更数学”的 1 作为索引开始。 （设置n=100并忽略求解器提供的数组长度也是错误的。）

考虑到标准算术运算是按元素应用于 numpy 数组的，您可以避免所有这些并缩短您的例程，所以

def testmodel(t, y):
    X,Y,J,Q = y
    a = 3; S = 0.4; K = 0.8; L = 2.3
    dXdt = K**a + Q**a - S*X
    dYdt = K*X - L*Y
    dJdt = S*Y - K*Q
    dQdt = K*X/L + J
    return dXdt, dYdt, dJdt, dQdt

为具有相同动态的多个隔间修改代码

您需要向求解器传递状态的平面向量。第一个设计决策是隔间及其组件如何排列在平面矢量中。与现有代码最兼容的一种变体是将相同的组件聚集在一起。那么在ODE函数中，第一个操作就是分离出这些簇。

    X,Y,J,Q = y.reshape([4,-1])

这会将输入向量分成 4 个长度相等的部分。最后，您需要反转此拆分，以便导数再次位于平面向量中。

    return np.concatenate([dXdt, dYdt, dJdt, dQdt])

其他一切都保持不变。除了初始向量，它需要有 4 个长度为 N 的段，其中包含隔间的数据。这里可能就是

    y0 = np.zeros(4*N)

如果初始数据来自任何其他来源，并且在每个隔间的记录中给出，您可能必须在展平之前转置结果数组。

请注意，此构造不是矢量化的，因此请不要在其默认的False 中设置该选项。

对于像圆圈这样的统一交互模式，我建议使用numpy.roll 来继续避免使用显式循环。对于看起来像网络的交互模式，可以使用连接矩阵和掩码，如Using python built-in functions for coupled ODEs