在python中使用forloop求解耦合微分方程

Question

我正在试着解一组微分方程，但我一直很难做到这一点。我的微分方程包含一个表示从1到n的数字的"i“下标。我尝试实现一个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，并且显式方法不使用雅可比矩阵。因此，不需要对偏导数的差商进行多点评估。

实际上，坐标数组的大小始终为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中一样