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的数组

看到问题就在眼前会有帮助