在solve_ivp()过程中调整求解变量

Question

我正在使用solve_ivp()函数来求解一些运动方程。更具体地说，我正在对社会力模型进行建模，其中solve_ivp（）以数值方式计算模拟中每个人的速度和位置。现在我遇到了一个奇怪的问题，我试图在求解过程中实现一个相对简单的操作，但结果与我的预期相反。

此代码片段显示了用于solve_ivp()函数的相关函数。

def motion(t, X):
        
        all_forces = Forces(X)
    
        for i in range(N):
            
            # allocate vector dX
            dX[i*MMM] = X[i*MMM+2]
            dX[i*MMM+1] = X[i*MMM+3]
            dX[i*MMM+2] = (1/p[i].mass)*all_forces[i,0]
            dX[i*MMM+3] = (1/p[i].mass)*all_forces[i,1]
        
        return dX

然后在“Forces(X)”函数中，我使用以下代码来检查人们是否在污染区域内，如果是，我想将他们的速度降低一个因子：

# update pedestrian variables
    for j in range(N):
        
        # check if the pedestrians are inside the contaminated area
        if np.linalg.norm(c[0].pos-p[j].pos[0]) <= c[0].radius:
            p[j].pos[0,0] = X[j*MMM]
            p[j].pos[0,1] = X[j*MMM+1]
            p[j].velocity[0] = X[j*MMM+2]*c[0].weight
            p[j].velocity[1] = X[j*MMM+3]*c[0].weight
        
        else:
            p[j].pos[0,0] = X[j*MMM]
            p[j].pos[0,1] = X[j*MMM+1]
            p[j].velocity[0] = X[j*MMM+2]
            p[j].velocity[1] = X[j*MMM+3]

此代码位于“Forces(X)”函数的顶部，然后根据人员的位置和速度值再次计算人员所受的力。

结果应该相当简单，因为速度只是乘以一个数字，但结果并没有表明这一点。问题是我使用的权重似乎是相反使用的，低权重（0.4）会增加它们的速度，而较高的权重（2）会降低它们的速度。我添加了 2 个 gif 来展示这一点。有谁明白为什么会发生这种情况？我正在寻找力的计算中可能导致此问题的原因，但想知道它是否符合我正在监督的代码的逻辑。

此 gif 显示了将权重设置为 0.4 的结果，您可以看到速度增加。

来自

的评论

c1_pos = np.array([18,5,7.5])
c1_radius = 3.5
c1_weight = 1.4

# initialize contaminant
c = []
c.append(contaminant(c1_pos, c1_radius, c1_weight))

Answer 1

我发现在Forces(X)函数中写入条件并没有得到想要的结果。我认为发生的事情是它会改变速度，然后根据改变的速度计算新的速度并更新它。由于 Forces(X) 中的其中一个力是通过当前速度与“所需速度”之间的差来计算的，“所需速度”是一个常数，这意味着计算出的力会变大。解决方案是将检查语句放入 Motion(t, X) 函数中，如下所示。这段代码只是添加了一些额外的 if 语句来检查各种半径。

    def motion(t, X):
        
        all_forces = Forces(X)
    
        for i in range(N):
            # check if the pedestrians are inside the contaminated area
            if np.linalg.norm(c[0].pos-p[i].pos[0]) <= c[0].radius:
                dX[i*MMM] = X[i*MMM+2]*c[0].weight
                dX[i*MMM+1] = X[i*MMM+3]*c[0].weight
                dX[i*MMM+2] = (1/p[i].mass)*all_forces[i,0]
                dX[i*MMM+3] = (1/p[i].mass)*all_forces[i,1]
            elif np.linalg.norm(c[1].pos-p[i].pos[0]) <= c[1].radius:
                dX[i*MMM] = X[i*MMM+2]*c[1].weight
                dX[i*MMM+1] = X[i*MMM+3]*c[1].weight
                dX[i*MMM+2] = (1/p[i].mass)*all_forces[i,0]
                dX[i*MMM+3] = (1/p[i].mass)*all_forces[i,1]
            elif np.linalg.norm(c[2].pos-p[i].pos[0]) <= c[2].radius:
                dX[i*MMM] = X[i*MMM+2]*c[2].weight
                dX[i*MMM+1] = X[i*MMM+3]*c[2].weight
                dX[i*MMM+2] = (1/p[i].mass)*all_forces[i,0]
                dX[i*MMM+3] = (1/p[i].mass)*all_forces[i,1]
            else:
                dX[i*MMM] = X[i*MMM+2]
                dX[i*MMM+1] = X[i*MMM+3]
                dX[i*MMM+2] = (1/p[i].mass)*all_forces[i,0]
                dX[i*MMM+3] = (1/p[i].mass)*all_forces[i,1]
            

        return dX