使用Numpy实现绳索模拟出现振荡爆炸，疑存数据竞争及向量化困惑

我尝试用Numpy编写基础绳索模拟程序，但程序运行时出现振荡甚至崩溃的情况，猜测可能存在数据竞争问题。

边被存储为(2,N)的整数数组，对应顶点的索引。代码中的一些矩阵操作可能导致同一顶点被多次引用，例如constraintpositions[edges[:,1]] -= edgevertexdeltasweighted，这会引发数据竞争吗？

我的完整脚本如下：

import bpy
import numpy as np
import time

mesh_object = bpy.context.active_object

# Initialise timestep (seconds), vertex weight (kg), gravity
timestep = 0.5
weight = 0.05
gravity = np.array([0,0,-9.8])

# Ensure it's a mesh and convert to numpy arrays
if mesh_object and mesh_object.type == 'MESH':
    mesh = mesh_object.data
    
    # Extract vertex coordinates as NumPy array
    length = len(mesh.vertices)
    positions = np.empty(length*3, dtype=np.float64)
    mesh.vertices.foreach_get('co', positions)
    positions.shape = (length, 3)
    
    # Extract edges as NumPy array
    edgelength = len(mesh.edges)
    edges = np.empty(edgelength*2, dtype=np.uint8)
    mesh.edges.foreach_get('vertices', edges)
    edges.shape = (edgelength, 2)

    # get edge lengths
    edgelengthsstart = np.empty(edgelength, dtype = np.float64)
    edgevectors = positions[edges[:,1]] - positions[edges[:,0]]
    edgelengthsstart = np.linalg.norm(edgevectors, axis= 1)

    #initialise arrays for distance constraint
    edgevertexdeltas = np.empty((edgelength, 2), dtype = np.float64)
    edgelengthsnew = np.empty(edgelength, dtype = np.float64)
    constraintpositions = np.zeros_like(positions)


else:
    print("Please select a valid mesh object.")


#initialise velocities
velocities = np.zeros_like(positions)
#print(velocities)

def update(velocities, positions, timestep, weight, gravity):

    #Update velocities and positions
    newpositions=np.copy(positions)
    velocities+= timestep*weight*gravity
    newpositions += timestep*velocities

    # calculate constraints

    for i in range (1,4):
        #v1, v2 = 
        edgevectorsnew = newpositions[edges[:,1]] - newpositions[edges[:,0]]
        edgelengthsnew = np.linalg.norm(edgevectorsnew, axis= 1)
        edgedifferences = edgelengthsnew - edgelengthsstart
        edgeratios = edgedifferences/edgelengthsnew
        #print(edgeratios)
        edgevertexdeltasweighted = 0.5 * edgevectorsnew * edgeratios[:, np.newaxis]

        #Update constraint positions
        constraintpositions = np.copy(newpositions)
        constraintpositions[edges[:,1]] -= edgevertexdeltasweighted
        constraintpositions[edges[:,0]] += edgevertexdeltasweighted
        newpositions = np.copy(constraintpositions)
        newpositions[0] = [0,0,0]


    #update velocity and positions
    velocities = (1/timestep)*(newpositions-positions)
    positions = np.copy(newpositions)
    #print(positions)

    # Put it back into mesh
    vertexpositions = np.copy(positions)
    vertexpositions.shape = length*3
    mesh.vertices.foreach_set("co", vertexpositions)
    mesh.update()


    return velocities, positions

for i in range(0,2):   
    velocities, positions = update(velocities, positions, timestep, weight, gravity)
    bpy.ops.wm.redraw_timer(type='DRAW_WIN_SWAP', iterations=1)
    time.sleep(0.1) 

另外，使用for循环逐个遍历边的版本运行正常，但弹性过大。我不清楚如何对边约束进行向量化处理。

内容的提问来源于stack exchange，提问作者ollie d