Python复刻Matlab质子吸引电子模拟的numpy数组操作错误排查

问题背景 #

通过复刻本科时期的Matlab质子吸引电子模拟项目提升Python技能，项目采用多维数组存储各时间步的粒子位置与速度。使用numpy数组实现时，受Matlab数组思维影响，在循环中调用np.append更新电子位置（re）和速度（v）数组时，先后出现axis=1轴越界、axis=0维度不匹配的错误。

原始代码 #

#This code is for plotting the trajectory of an electron travelling through space
#The particle is in the vacinity of a central force (A proton at the origin)

import numpy as np
import matplotlib.pyplot as plt

re = np.array([[1e-11],[1e-11]]) #let re denote the trajectory of the electron with x = r[0] and y = r[1]
m = 9.11e-31 #mass of electron
k = 8.99e9 #Coulomb's constant [N m^2/C^2]
q = 1.6e-19 #charge of electron and proton [C]
rp = [0,0] #rp is the position of the proton
dt = 0.001 #time differential [s]

v = np.array([[-3e12],[0]]) #the electron has initial velocity of v = (-3, 0) [m/s]
phi = np.arctan2(re[1][0], re[0][0]) #starting angle

for i in range(1,10):
   # nrex = (re[0][i-1])+v[0][i-1]*dt #nuew position in x
   # nrey = (re[1][i-1])+v[1][i-1]*dt #new position in y
    re[0] = np.append(re[0], ((re[0][i-1])+v[1][i-1]*dt), axis=1) #, axis=1) #for each timestep move the velocity in x
    re[1] = np.append(re[1], ((re[1][i-1])+v[1][i-1]*dt), axis=1) #for each timestep mobe the velocity in y

    phi = np.arctan2(re[1][i],re[0][i]) #update the angle
    rho = np.sqrt(re[0][i]**2 + re[1][i]**2) #update separation from proton
    v[0] = np.append(n[0], (v[0][i-1]+((k*(q**2)/(rho**2))/m)*np.cos(phi)*dt), axis=1) #update velocity in x
    v[1] = np.append(v[1], (v[1][i-1]+((k*(q**2)/(rho**2))/m)*np.sin(phi)*dt), axis=1) #update velocity in y
    
plt.scatter(re[0][:], re[1][:], s=2, c='b') #Plot electron's trajectory
plt.scatter(rp[0],rp[1], s=3, c='r') #Show proton's position
plt.show() #Show

错误信息 #

1. 使用axis=1时的错误：

Traceback (most recent call last):
  File "c:\Users\hecto\OneDrive\Documentos\ITESM\8vo semestre\Repaso Python\ParticleInElectricField.py", line 20, in <module>    
    re[0] = np.append(re[0], ((re[0][i-1])+v[1][i-1]*dt), axis=1) #, axis=1) #for each timestep move the velocity in x
  File "<__array_function__ internals>", line 200, in append
  File "C:\Users\hecto\AppData\Local\Programs\Python\Python38\lib\site-packages\numpy\lib\function_base.py", line 5499, in append
    return concatenate((arr, values), axis=axis)
  File "<__array_function__ internals>", line 200, in concatenate
numpy.AxisError: axis 1 is out of bounds for array of dimension 1

2. 改用axis=0时的错误：

Traceback (most recent call last):
  File "c:\Users\hecto\OneDrive\Documentos\ITESM\8vo semestre\Repaso Python\ParticleInElectricField.py", line 20, in <module>
    re[0] = np.append(re[0], ((re[0][i-1])+v[1][i-1]*dt), axis=0) #, axis=1) #for each timestep move the velocity in x
  File "<__array_function__ internals>", line 200, in append
  File "C:\Users\hecto\AppData\Local\Programs\Python\Python38\lib\site-packages\numpy\lib\function_base.py", line 5499, in append        
    return concatenate((arr, values), axis=axis)
  File "<__array_function__ internals>", line 200, in concatenate
ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 1 dimension(s) and the array at index 1 has 0 dimension(s)

错误原因分析 #

• axis=1越界：re[0]是从初始2D数组re（形状(2,1)）中提取的行，本质是1维数组（形状(1,)），而axis=1仅适用于2D及以上数组，因此触发轴越界错误。
• 维度不匹配：np.append要求两个输入数组维度一致，代码中第二个参数是计算得到的标量（0维），而re[0]是1维数组，维度不匹配导致错误。
• 核心思维差异：Matlab支持动态扩展数组且维度处理更宽松，但Numpy数组是固定大小的，频繁使用np.append不仅效率低，还容易触发维度问题。

正确实现方案 #

方案1：预分配数组（推荐，高效符合Numpy最佳实践） #

先确定总时间步数，预创建固定大小的2D数组，直接通过索引赋值更新每一步的位置和速度，避免动态扩展的性能损耗和维度问题。

方案2：动态扩展（仅需时使用） #

若必须动态扩展，需保证np.append的两个输入维度一致：

• 将新增的标量转为与原数组同维度的数组（如np.array([new_value])）
• 操作整个2D数组而非单独行，保持维度统一

修正后的代码（预分配数组版本） #

import numpy as np
import matplotlib.pyplot as plt

# 模拟参数
total_steps = 10  # 总时间步数
re = np.zeros((2, total_steps))  # 预分配位置数组：行=坐标(x/y)，列=时间步
re[:, 0] = [1e-11, 1e-11]  # 初始位置
m = 9.11e-31
k = 8.99e9
q = 1.6e-19
rp = [0, 0]
dt = 0.001

v = np.zeros((2, total_steps))  # 预分配速度数组
v[:, 0] = [-3e12, 0]  # 初始速度

for i in range(1, total_steps):
    # 更新位置：基于上一步速度
    re[0, i] = re[0, i-1] + v[0, i-1] * dt
    re[1, i] = re[1, i-1] + v[1, i-1] * dt

    # 计算受力角度和距离
    phi = np.arctan2(re[1, i], re[0, i])
    rho = np.sqrt(re[0, i]**2 + re[1, i]**2)
    # 库仑力加速度
    a = (k * q**2) / (m * rho**2)

    # 更新速度：基于库仑力
    v[0, i] = v[0, i-1] + a * np.cos(phi) * dt
    v[1, i] = v[1, i-1] + a * np.sin(phi) * dt

# 绘制轨迹
plt.scatter(re[0, :], re[1, :], s=2, c='b')
plt.scatter(rp[0], rp[1], s=30, c='r')  # 放大质子点便于观察
plt.axis('equal')  # 保持坐标轴比例一致，避免轨迹变形
plt.show()

关键改进点 #

• 预分配(2, total_steps)的数组，直接通过索引赋值，效率远高于动态np.append
• 修正了原始代码中位置更新的错误（原代码误用v[1]计算x方向位置）
• 添加plt.axis('equal')保证轨迹显示比例正确
• 简化索引操作，更符合Numpy的数组访问习惯

内容的提问来源于stack exchange，提问作者Hector Viruega