你的问题出在numpy的广播规则和MATLAB的运算逻辑细微差异上，我来帮你拆解和解决：

问题原因 #

在你的循环代码中：

for n in range(31):
    b[:, n] = a[:, n] ** c

• a[:, n]是形状为(40,)的一维数组，numpy广播时会将其视为"行向量"
• c是形状为(40,1)的二维列数组

执行幂运算时，numpy的广播规则会把这两个数组扩展成(40,40)的二维矩阵，得到的结果是40×40的数组，但b[:, n]只能接受40个元素的一维数组，因此抛出ValueError: could not broadcast input array from shape (40,40) into shape (40)。

两种解决方案 #

方案1：调整维度适配循环逻辑

如果你想保留循环，可以修改其中一个数组的维度，让运算结果是一维数组：

• 把c转成一维数组（用c.ravel()或c[:,0]）：

import numpy as np
a = np.arange(1240).reshape(40, 31)
b = np.ones(a.shape)
c = np.ones((40, 1))

for n in range(31):
    b[:, n] = a[:, n] ** c.ravel()  # 将c转为(40,)的一维数组

• 或者把a[:,n]转成二维列数组（用a[:,n, None]）：

for n in range(31):
    b[:, n] = a[:, n, None] ** c  # 将a的列转为(40,1)，运算后结果为(40,1)，赋值时自动降维

方案2：利用numpy广播特性，去掉循环（更高效）

numpy的广播天生支持这种矩阵与列向量的逐列运算，完全不需要写循环，直接一行代码搞定，和MATLAB的向量化风格一致：

import numpy as np
a = np.arange(1240).reshape(40, 31)
c = np.ones((40, 1))
b = a ** c  # 直接对整个矩阵运算，自动广播c到(40,31)的形状

这里a是(40,31)，c是(40,1)，numpy会自动把c的每一行重复31次，扩展成和a相同的形状，然后逐元素执行幂运算，结果就是你想要的(40,31)矩阵，比循环快得多，也更符合numpy的最佳实践。

验证小例子 #

用更小的数组测试效果：

a = np.arange(8).reshape(4,2)  # 形状(4,2)
c = np.array([[2],[2],[2],[2]])  # 形状(4,1)
b = a ** c
print(b)
# 输出：
# [[ 0  1]
#  [ 4  9]
#  [16 25]
#  [36 49]]

这个结果和你逐列循环运算的结果完全一致。

内容的提问来源于stack exchange，提问作者B Furtado