我最近在做一个积分对比测试，遇到了一个让我困惑的问题，希望能得到解答：

首先我定义了一个一维分段函数：

• 当r < 1时，p(r) = 1 - r
• 否则p(r) = 0

接着我计算了它的一维积分I1(r) = ∫p dr，积分从0开始，并且设置积分常数使得I1(1) = 0。

现在我把r扩展到二维x-y网格中，作为原点到点(x,y)的距离，也就是r(x,y) = √(x² + y²)。我想在这个二维网格上计算p的积分，于是根据微分关系dr = (∂r/∂x)dx + (∂r/∂y)dy，推导出二维积分的形式应该是：
I2(x,y) = ∫p(∂r/∂x)dx + ∫p(∂r/∂y)dy
我设置积分下限为任何小于-1的点，同时让r=1的轮廓上I2=0来确定积分常数。

理论上，虽然I1(r)是一维函数，I2(x,y)是二维函数，不能直接代数比较，但它们的极值应该是匹配的。可当我用Python代码测试时，却发现I2的极值大约是I1极值的2倍；而且我在3D网格上做同样的测试时，这个倍数变成了3倍，明显和维度相关。

以下是我用来测试的Python 3代码：

import numpy as np
import scipy.integrate
import matplotlib.pyplot as plt

nr = 101
r = np.linspace(0, 2, nr)
p = 1 - r
p *= (p > 0)    # 将r≥1区域的p设为0
I1 = scipy.integrate.cumulative_simpson(p, x=r, initial=0)
I1 -= I1[-1]  # 确保I1(1)=0

nx, ny = 151, 301
x = np.linspace(-1.5, 1.5, nx)
y = np.linspace(-2.0, 2.0, ny)
# 计算x=0和x=1对应的索引
ix0 = round((0.0-x[0])/(x[-1]-x[0])*(nx-1))
ix1 = round((1.0-x[0])/(x[-1]-x[0])*(nx-1))
iy0 = round((0.0-y[0])/(y[-1]-y[0])*(ny-1))  # y=0对应的索引

x2, y2 = np.meshgrid(x, y, indexing='xy')
r2 = np.sqrt(x2**2 + y2**2)
# 计算r对x和y的偏导，处理原点的NaN
dr2dx = x2 / r2
dr2dx[iy0, ix0] = 0.0
dr2dy = y2 / r2
dr2dy[iy0, ix0] = 0.0

# 生成二维网格上的p函数
p2 = (1 - r2)
p2 *= (p2 > 0)

# 计算二维积分I2
I2 = scipy.integrate.cumulative_simpson(p2 * dr2dx, x=x, axis=1, initial=0)
I2 += scipy.integrate.cumulative_simpson(p2 * dr2dy, x=y, axis=0, initial=0)
I2 -= I2[iy0, ix1]  # 确保r=1时I2=0
I2 = I2 * (I2 < 0)  # 去除数值误差带来的残留正值

print(f'I1的最小值: {np.min(I1)}\nI2的最小值: {np.min(I2)}')

# 绘图可视化
fig, axs = plt.subplots(1, 3, figsize=(12, 4))
axs[0].plot(r, p, label='p')
axs[0].plot(r, I1, label=r'I1')
axs[0].legend()
axs[0].set_xlabel('r')

c1 = axs[1].contour(x2, y2, p2)
axs[1].set_aspect(1)
axs[1].set_xlabel('x')
axs[1].set_ylabel('y')
axs[1].set_title('二维p函数')

c2 = axs[2].contour(x2, y2, I2)
axs[2].set_aspect(1)
axs[2].set_xlabel('x')
axs[2].set_ylabel('y')
axs[2].set_title('二维积分I2')

fig.colorbar(c1)
fig.colorbar(c2)
plt.tight_layout()
plt.show()
plt.close()

我还尝试过缩小积分区域（比如只在单个象限、上半平面或右半平面计算），修改后的I2计算代码如下，但结果依然显示极值存在2倍的偏差：

I2x = scipy.integrate.cumulative_simpson(p2 * dr2dx, x=x, axis=1, initial=0)
for j in range(ny):
    I2x[j,:] -= I2x[j,-1]    # 为每个y行单独设置积分常数
I2y = scipy.integrate.cumulative_simpson(p2 * dr2dy, x=y, axis=0, initial=0)
for i in range(nx):
    I2y[:,i] -= I2y[-1,i]    # 为每个x列单独设置积分常数
I2 = I2x + I2y
I2 = I2 * (I2 < 0)    # 去除数值残留

我现在的疑问是：

• 是不是我对微分关系dr = (∂r/∂x)dx + (∂r/∂y)dy的应用有误？
• 还是二维积分的计算方式本身就存在错误？
• 为什么极值的倍数刚好等于空间的维度（2倍对应2D，3倍对应3D）？

备注：内容来源于stack exchange，提问作者udy11