仿真报错ValueError：Operands无法广播（形状不匹配）求助

错误根源分析 #

报错ValueError: operands could not be broadcast together with shapes (9999,10000) (10000,10000) (9999,10000)来自两处核心问题：

• 电场更新时，切片后的E与补全后的Hx/Hy形状不匹配
• 孔径函数的平方运算符写错，且未适配numpy数组的向量化计算逻辑

具体修复步骤 #

1. 修正电场更新的形状匹配问题 #

补全后的Hx/Hy是(num_points, num_points)形状，而E[:-1,:]和E[:,:-1]是(num_points-1, num_points)和(num_points, num_points-1)，需要对Hx/Hy做对应切片：

# Update electric field
E[:-1,:] -= time_step * Hx[:-1,:]  # 取Hx前num_points-1行匹配E的切片
E[:,:-1] += time_step * Hy[:,:-1]  # 取Hy前num_points-1列匹配E的切片

2. 修复孔径函数的错误并向量化 #

原函数中x*2/y*2是平方运算符笔误，应为x**2/y**2；同时不能用标量if判断numpy数组，改用向量化的np.where实现：

# Define aperture function
def aperture(x, y):
    r = np.sqrt(x**2 + y**2)
    return np.where(r <= aperture_diameter / 2, 1, 0)

3. 优化孔径掩码的计算位置 #

x/y/xx/yy是固定不变的网格坐标，无需在每次时间循环内重复生成，将其移到循环外减少冗余计算：

# Precompute aperture mask once (outside loop)
x = np.linspace(-screen_distance/2, screen_distance/2, num_points)
y = np.linspace(-screen_distance/2, screen_distance/2, num_points)
xx, yy = np.meshgrid(x, y, indexing='ij')
aper = aperture(xx, yy)

# Main simulation loop
for t in range(num_time_steps):
    # ... 其他代码 ...
    # Apply aperture function
    E *= aper

修改后的完整代码 #

import numpy as np
import matplotlib.pyplot as plt

# Define simulation parameters
wavelength = 1e-3   # meters
aperture_diameter = 0.5   # meters
screen_distance = 5   # meters
grid_spacing = wavelength / 2  # meters
num_points = int(screen_distance / grid_spacing)
# 遵循FDTD的CFL稳定性条件（简化模型假设波速为1，实际应用可替换为光速c=3e8）
time_step = grid_spacing / (2 * np.sqrt(2))   # seconds
# 计算波传播到屏幕的时间对应步数
num_time_steps = int(screen_distance / (time_step * 1))  

# Define aperture function
def aperture(x, y):
    r = np.sqrt(x**2 + y**2)
    return np.where(r <= aperture_diameter / 2, 1, 0)

# Initialize electric field array
E = np.zeros((num_points, num_points), dtype=np.complex128)

# Precompute aperture mask once to save computation
x = np.linspace(-screen_distance/2, screen_distance/2, num_points)
y = np.linspace(-screen_distance/2, screen_distance/2, num_points)
xx, yy = np.meshgrid(x, y, indexing='ij')
aper = aperture(xx, yy)

# Main simulation loop
for t in range(num_time_steps):
    # Update magnetic field
    Hx = np.diff(E, axis=0) / grid_spacing
    Hy = np.diff(E, axis=1) / grid_spacing
    Hx = np.pad(Hx, ((0, 1), (0, 0)), mode='constant')
    Hy = np.pad(Hy, ((0, 0), (0, 1)), mode='constant')
    
    # Update electric field with shape-matched slices
    E[:-1,:] -= time_step * Hx[:-1,:]
    E[:,:-1] += time_step * Hy[:,:-1]
    
    # Apply aperture function
    E *= aper

# Calculate intensity at screen
I = np.abs(E)**2

# Plot results
plt.imshow(I, cmap='gray', extent=[-screen_distance/2, screen_distance/2, -screen_distance/2, screen_distance/2])
plt.xlabel('Position (m)')
plt.ylabel('Position (m)')
plt.title('Diffraction pattern for circular aperture')
plt.show()

内容的提问来源于stack exchange，提问作者Sophie Er