球极坐标雅可比除法的解决方法

Question

我遇到了数据可视化问题。 简短的总结，我正在研究一个涉及极球坐标网格的项目，并试图解决每个细胞的 ODE（化学反应）耦合系统。 出于特定原因，我需要我的状态向量的形式为 (r^2*sin(theta)*n_i), i={1,2,3...}。

我在下面重写了我的问题的一个简单示例，您可以按原样运行它。 为什么我认为等于 np.ones(a[0].shape) 的 cst2 没有显示统一的 pcolormesh。 更令我惊讶的是，为什么添加颜色条会使这个问题消失？

我最好的猜测是除以 r^2*sin(theta) 会导致数值问题，但我该如何解决这个问题？ （我需要查看没有曲率项的数据才能对其进行解释 -> 雅可比除法对我来说似乎是强制性的）。

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

### Edges
r = np.logspace(np.log10(1), np.log10(4.6), num=14) #cell edges
theta = np.linspace(0+0.001,np.pi-0.001,num=10) 
b = np.meshgrid(r,theta)

### Center
r_c = r[0:-1] + np.ediff1d(r)/2 #get the cell center
theta_c = theta[0:-1] + np.ediff1d(theta)/2
a = np.meshgrid(r_c,theta_c)

### The jacobian division
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')
# clb = fig.colorbar(pcm, ax=ax, orientation='horizontal') 

Answer 1

只是缺少括号的问题。这样做：

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