球极坐标雅可比除法pcolormesh可视化异常咨询

问题背景 #

• 项目基于球极坐标网格开展计算，需要为每个网格单元求解对应化学反应过程的耦合ODE方程组，状态向量要求采用r²*sin(theta)*n_i（i取1,2,3...）的形式。
• 按照理论推导，经雅可比除法得到的cst2应为和a[0]形状一致的全1数组，但实际绘制pcolormesh伪彩图时并未呈现均匀纯色，而添加colorbar后该显示异常会直接消失。
• 初步猜测异常由除以r²*sin(theta)的操作引入数值误差导致，由于需要去除曲率项才能完成结果解读，雅可比除法是必要处理步骤，需要找到可落地的规避方案。

最小复现代码 #

import numpy as np
import matplotlib.pylab as plt
fig, ax = plt.subplots()

### 网格边界定义
r = np.logspace(np.log10(1), np.log10(4.6), num=14) # 径向网格边坐标
theta = np.linspace(0+0.001,np.pi-0.001,num=10) # 极角网格边坐标
b = np.meshgrid(r,theta)

### 网格中心计算
r_c = r[0:-1] + np.ediff1d(r)/2 # 径向单元中心坐标
theta_c = theta[0:-1] + np.ediff1d(theta)/2 # 极角单元中心坐标
a = np.meshgrid(r_c,theta_c)

### 雅可比除法操作
cst = pow(a[0],2)*np.sin(a[1])
cst2 = np.copy(cst)/pow(a[0],2)/np.sin(a[1])
pcm = ax.pcolormesh(b[0]*np.cos(b[1]),
                    b[0]*np.sin(b[1]), 
                    cst2,cmap='seismic',edgecolor='black')
# 取消下方注释添加colorbar后，显示异常会消失
# clb = fig.colorbar(pcm, ax=ax, orientation='horizontal')

问题根因 #

• cst2并非严格的全1数组：双精度浮点数的乘除运算存在固有机器精度误差，实际cst2的取值在1±1e-15区间内波动，偏差属于浮点运算的正常范围，不属于计算逻辑错误。
• 未添加colorbar时，pcolormesh默认以传入数据的全局最小值、最大值作为色阶映射的上下限，此时数据极差仅为1e-15量级，极微小的数值波动会被拉伸到seismic色图从蓝到红的全映射范围，叠加黑色网格边线的视觉干扰，就会呈现出颜色不均的异常效果。
• 添加colorbar会触发matplotlib的色阶重计算逻辑，由于数据和1的偏差极小，重算时会自动将色阶范围钳位到包含1的合理区间，浮点误差带来的颜色差异被压缩到人眼无法识别的程度，异常视觉效果就消失了。

可行规避方案 #

• 方案1：显式固定颜色映射范围，从根源避免微小误差被色阶拉伸。绘制伪彩图时手动传入vmin、vmax参数，根据实际物理量的合理区间设置色阶上下限，对于本应全1的测试数据可参考如下写法：

  pcm = ax.pcolormesh(b[0]*np.cos(b[1]),
                      b[0]*np.sin(b[1]), 
                      cst2,cmap='seismic',edgecolor='black',
                      vmin=0.99, vmax=1.01)
  

  该方案是最稳妥的通用方案，不会修改原始计算数据，仅通过可视化参数调整解决显示问题。
• 方案2：对雅可比除法后的结果做精度截断，过滤机器精度级别的浮点误差。使用np.round将结果保留12位左右有效数字，即可消除1e-15量级的数值波动：

  cst2 = cst / pow(a[0],2) / np.sin(a[1])
  cst2 = np.round(cst2, decimals=12)
  

• 方案3：针对去除曲率项后的物理量做校准，若理论上某部分结果应为常数，可在计算完成后通过归一化操作将常数项校准为精确值，避免多步计算累积的浮点误差影响后续可视化。

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