ValueError: операнды не могут транслироваться вместе с фигурами (3,) (3000,)

Question

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

Я мог понять, что x и xAcc не имеют одинаковых размеров, но я подумал, что исправил это со строкой xPositions[i, :] = x.

import numpy as np

np.seterr(invalid="ignore")

m = 1

x = np.array([1, 5, 9])
y = np.array([16, 20, 24])

def GetLJForce(r, epsilon, sigma):
    return 48 * epsilon * np.power(sigma, 12) / np.power(r, 13) - 24 * epsilon * np.power(sigma, 6) / np.power(r, 7)

def GetAcc(xPositions, yPositions):
    global xAcc
    global yAcc
    xAcc = np.zeros((xPositions.size, xPositions.size), dtype=object)
    yAcc = np.zeros((xPositions.size, xPositions.size), dtype=object)

    for i in range(0, xPositions.shape[0]-1):
        for j in range(i+1, xPositions.shape[0]-1):
            r_x = xPositions[j] - xPositions[i]
            r_y = yPositions[j] - yPositions[i]

            rmag = np.sqrt(r_x*r_x + r_y*r_y)

            if(rmag[0]==0 or rmag[1]==0 or rmag[2]==0):
                rmag += 1

                force_scalar = GetLJForce(rmag, 0.84, 2.56)
            
                force_x = force_scalar * r_x / rmag
                force_y = force_scalar * r_y / rmag

                xAcc[i,j] = force_x / m
                xAcc[j,i] = - force_x / m
            
                yAcc[i,j] = force_y / m
                yAcc[j,i] = - force_y / m
            
            else:
                force_scalar = GetLJForce(rmag, 0.84, 2.56)

                force_x = force_scalar * r_x / rmag
                force_y = force_scalar * r_y / rmag

                xAcc[i,j] = force_x / m
                xAcc[j,i] = - force_x / m

                yAcc[i,j] = force_y / m
                yAcc[j,i] = - force_y / m

    return np.sum(xAcc, axis=0), np.sum(yAcc, axis=0)

def UpdatexPos(x, v_x, a_x, dt):
    return x + v_x*dt + 0.5*a_x*dt*dt

def UpdateyPos(y, v_y, a_y, dt):
    return y + v_y*dt + 0.5*a_y*dt*dt

def UpdatexVel(v_x, a_x, a1_x, dt):
    return v_x + 0.5*(a_x + a1_x)*dt

def UpdateyVel(v_y, a_y, a1_y, dt):
    return v_y + 0.5*(a_y + a1_y)*dt

def RunMD(dt, number_of_steps, x, y):
    xPositions = np.zeros((number_of_steps, 3))
    yPositions = np.zeros((number_of_steps, 3))

    v_x = 0
    v_y = 0

    a_x = GetAcc(xPositions, yPositions)[0]
    a_y = GetAcc(xPositions, yPositions)[1]

    for i in range(number_of_steps):
        x = UpdatexPos(x, v_x, a_x, dt)
        y = UpdateyPos(y, v_y, a_y, dt)

        a1_x = GetAcc(xPositions, yPositions)[0]
        a1_y = GetAcc(xPositions, yPositions)[1]

        v_x = UpdatexVel(v_x, a_x, a1_x, dt)
        v_y = UpdateyVel(v_y, a_y, a1_y, dt)

        a_x = np.array(a1_x)
        a_y = np.array(a1_y)

        xPositions[i, :] = x
        yPositions[i, :] = y

    return xPositions, yPositions

sim_xpos = RunMD(0.1, 1000, x, y)[0]
sim_ypos = RunMD(0.1, 1000, x, y)[1]

np.savetxt("atomsx1.txt", sim_xpos)
np.savetxt("atomsy1.txt", sim_ypos)

Answer 1

Исходя из вашего кода, я не могу понять, почему нужно объявлять global xAcc таким образом. Кажется, вы не используете его нигде за пределами функции.

В любом случае строка xPositions[i, :] = x не поможет вам изменить ваш xAcc размер, потому что он был объявлен ранее в xPositions = np.zeros((number_of_steps, 3)) .

Я вижу, что xAcc = np.zeros((xPositions.size, xPositions.size), dtype=object) объявлен как квадратная матрица (3000, 3000), и в этой строке return np.sum(xAcc, axis=0), np.sum(yAcc, axis=0) вы уменьшили его только по одной оси, поэтому вы получаете xAcc как вектор длиной 3000. Возможно, в вашем случае поможет уменьшение на другой оси np.sum(xAcc), таким образом вы получите массив длины 3.

Было бы полезно увидеть проблему под рукой.