Python绘制含除零场景数学函数的报错与绘图异常排查

现有代码存在的错误 #

• 变量缺失：计算非零角度函数值时用到的sin_the没有定义，实际应为极角THETA的正弦值np.sin(THETA)，属于笔误。
• 冗余死代码：dB函数中return语句之后的频率数组定义永远不会执行；开头提前定义的4个角度值的theta变量，后续被全角度theta直接覆盖，属于无效代码。
• 警告触发原因：θ=0位置sinθ=0，直接计算时分母为0触发除零警告；除零得到的inf值和后续正弦项计算得到的0相乘得到nan，触发无效值乘法警告。虽然已经通过洛必达法则算出了θ=0处的极限，但代码没有对该位置做特殊赋值，直接裸跑计算必然出警告。
• 极坐标图不从0°起始的原因：绘图逻辑错误把频率拆成两组，前4个频率只绘制0_{179°区间，后4个频率从181°位置开始绘制，直接跳过了0°位置的点；且`theta[181:361]`属于越界切片——theta数组总共360个点（对应0}359°），索引最大值为359，这个切片取不到180°位置的点，进一步导致角度轴数据不连续。
• 数值逻辑错误：θ=0位置的比值Γ理论值为1（分子分母同取极限，比值恒为1），但代码中该位置值为nan，会导致绘图在0°位置断线。

修正后可运行代码 #

import numpy as np
import matplotlib.pyplot as plt

def find_nearest(array, value):
    array = np.asarray(array)
    idx = (np.abs(array - value)).argmin()
    return idx

def dB(x):
    return 20*np.log10(np.abs(x))

# 基础参数配置
c = 340  # 空气中声速，单位m/s
target_freqs = [125, 250, 500, 1000, 2000, 4000, 8000, 16000]
f = np.arange(start=125, stop=16001, step=1)
theta = np.deg2rad(np.arange(start=0, stop=360, step=1))
F, THETA = np.meshgrid(f, theta, indexing='ij')

k = 2 * np.pi * F / c  # 波数计算
sin_the = np.sin(THETA)  # 补全缺失的角度正弦变量

# 临时屏蔽计算过程中已知会出现的除零、无效值警告
with np.errstate(divide='ignore', invalid='ignore'):
    f_neq0 = (-5) / (k * sin_the) * np.sin(2.5 * k * sin_the)
f_0 = (-2.5 * 1j) / (np.pi * F)
Gamma = f_neq0 / f_0
# 手动修正θ=0位置的比值，该位置极限值为1
theta0_idx = find_nearest(theta, 0)
Gamma[:, theta0_idx] = 1

# 极坐标绘图
plt.figure(figsize=(10, 10))
ax = plt.subplot(111, polar=True)
ax.set_theta_zero_location('N')  # 设置0°位置在北向
ax.set_theta_direction(-1)  # 角度顺时针递增，符合声学指向性图常规习惯，不需要可删除
ax.set_rlabel_position(0)  # 径向刻度标签放在0°位置

# 遍历所有目标频率绘图，不再拆分角度段
for freq in target_freqs:
    f_idx = find_nearest(f, freq)
    ax.plot(theta, dB(Gamma[f_idx, :]), label=f'{freq}Hz')

ax.set_ylim(-40, 0)  # 可根据实际计算的dB范围调整纵轴上下限
plt.legend(loc='lower right')
plt.show()

关键修正点说明 #

• 补全所有缺失变量，删除所有无效死代码，解决基础语法问题。
• 用np.errstate上下文管理器屏蔽已知的计算警告，计算完成后手动给θ=0位置赋理论值1，从根源消除两类RuntimeWarning。
• 废除原来拆分角度段绘图的错误逻辑，所有频率统一绘制全角度范围的曲线，修复0°位置缺线、索引越界的问题，保证极坐标图从0°位置连续起始。
• 增加图例、轴范围配置，提升绘图可读性，可根据作业要求自行调整样式参数。

内容的提问来源于stack exchange，提问作者Sebastian N