Нестабильность массивной пружинной системы с 2 узлами и RK45

Question

Я хочу смоделировать деформацию системы массовой пружины с множеством узлов.Я сталкиваюсь с проблемой нестабильности в определенных условиях.

Для этого примера я моделирую влияние гравитации на узел весом 1 кг.

Когда мой узел 0 находится в положении 0., 0., 0. и мой узел 1 в положении 1., 1., 1. все работает нормально и эффект пружины согласован.

Когда мой узел 0 находится в позиции 0., 1., 0, а мой узел 1 в позиции 0., 0., 0. результат расходится.

Я не понимаю, почему в аналогичных условиях я получаю разные результаты только из-за положения узлов.Я использую решатель неправильно?Или мой расчет силы пружины неверен?Я пробовал как классический способ с циклом в цикле, так и векторизацию, и результаты одинаковы.

import numpy as np
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import scipy.integrate

X = np.array([[0., 0., 0.], # 0
             [1., 1., 1.]]) # 1

edges = np.array([[0, 1]])

W = np.array([[0, 1],
              [1, 0]])

d = cdist(X, X, "euclidean")

mass = 1 * np.ones(X.shape[0])
mass[0] = np.inf
mass[1] = 1
Force = np.zeros((X.shape[0], 3))
Force[1, 1] = 9.81
N = X.shape[0]

dt = 0.1
T = int(1e4)

k = 5


def compute_F_spring_vec(X, W, di):
    d_inv = np.reciprocal(cdist(X, X, "euclidean"))
    d_inv[np.arange(np.shape(d_inv)[0]),np.arange(np.shape(d_inv)[0])] = 0
    op1 = np.dot(W, X)
    op2 = X * np.repeat(np.sum(W, axis=1), 3).reshape(X.shape)
    op3 = np.dot(W * di * d_inv, X)
    op4 = X * np.repeat(np.sum(W * di * d_inv, axis=1), 3).reshape(X.shape)
    F_spring = - k * (-op1 + op2 + op3 - op4)

    return F_spring


def compute_F_spring(X, W, d):
    F_spring = []
    for i in range(X.shape[0]):
        F_i = np.zeros(3)
        for j in range(X.shape[0]):
            if W[i,j]:
                F_i += -k * (np.linalg.norm(X[i, :] - X[j, :]) - d[i, j]) * (X[i, :] - X[j, :])/np.linalg.norm(X[i, :] - X[j, :])

        F_spring.append(F_i)

    return F_spring

def fun(t, X):
    X = X.reshape((int(X.shape[0]/3), 3))
    F_spring = compute_F_spring_vec(X, W, d)
    F = F_spring + Force
    a = F / mass[:, np.newaxis]

    X = X + a*t
    return X.ravel()

solver = scipy.integrate.RK45(fun, 0, X.ravel(), t_bound=100, rtol=0.001, vectorized=True)

Ylast = np.inf
X1 = []

while solver.status == 'running' and solver.t < 10:
    error = solver.step()
    X1.append([solver.y[4], solver.t, solver.y[3], solver.y[5]])

    if (np.abs(X1[-1][0] - Ylast))/np.abs(Ylast) < 1e-6:
        break
    else:
        Ylast = X1[-1][0]

Еще немного точности:

X содержит позиции x, y, z всехузлы.Края - это связь между каждой парой узлов.Из этого я могу получить матрицу связности W. 1 - связь, 0 - нет.

d представляет евклидово расстояние между всеми узлами.