ValueError: operands could not be broadcast together with shapes (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的数组。

看到手头的问题会很有帮助。

ValueError: operands could not be broadcast together with shapes (3,) (3000,)

ValueError: operands could not be broadcast together with shapes (3,) (3000,)

python

physics