python - Missing factor of 2 when numerically integrating a 1d function on 2d grid - Stack Overflow

admin管理员组

文章数量:1394600

I have a 1d function p(r) = 1-r when r<1 and p(r) = 0 otherwise. Let I1(r) = ∫pdr where integration starts from 0. Let's set the constant of integration so that I1(1)=0. Now, assume r represents radius from origin on a 2d x-y grid so that r(x,y) = √(x²+y²). If I want to instead compute the integral of p on the 2d x-y grid, I can do it in the following manner: I2(x,y) = ∫p (∂r/∂x) dx + ∫p (∂r/∂y) dy because dr = (∂r/∂x) dx + (∂r/∂y) dy. The x,y integral limits can start at any point less than -1. Set the integration constant so that I2=0 at the contour r=1. Since I1(r) is a 1d function and I2(x,y) is a 2d function, they cannot be directly compared but if I plot them, I should see that their extrema are matching. However, when I wrote a simple code to test it, I noticed that extremum of I2 is about 2 times the extremum of I1. Following is the Python-3 code to test it

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)    # setting p=0 outside r=1
I1 = scipy.integrate.cumulative_simpson(p, x=r, initial=0)
I1 -= I1[-1]

nx, ny = 151, 301
x = np.linspace(-1.5, 1.5, nx)
y = np.linspace(-2.0, 2.0, ny)
ix0 = round((0.0-x[0])/(x[-1]-x[0])*(nx-1))    # index where x=0
ix1 = round((1.0-x[0])/(x[-1]-x[0])*(nx-1))    # index where x=1
iy0 = round((0.0-y[0])/(y[-1]-y[0])*(ny-1))    # index where y=0
x2, y2 = np.meshgrid(x, y, indexing='xy')
r2 = np.sqrt(x2**2 + y2**2)
dr2dx = x2 / r2    # ∂r/∂x
dr2dx[iy0,ix0] = 0.0    # to get rid of nan at (0,0)
dr2dy = y2 / r2    # ∂r/∂y
dr2dy[iy0,ix0] = 0.0    # to get rid of nan at (0,0)
p2 = (1 - r2)    # function p on x-y grid
p2 *= (p2 > 0)    # setting p=0 outside radius 1
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]
I2 = I2 * (I2 < 0)    # removes some numerical residues

print(f'Minima of I1: {np.min(I1)}\nMinima of I2: {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(r'I2')
fig.colorbar(c1)
fig.colorbar(c2)
plt.tight_layout(); plt.show(); plt.close()

Have I made a mistake somewhere? If not, what is the actual reason behind this and how to properly rectify it. It might appear that 2d integration is somehow happening twice but I don't see that I am doing anything mathematically wrong here. One can also say that integration is happening twice in both x and y directions, so why not a factor of 4? I also tried the 2d integration in a quadrant, upper half and right half but result was same. For these cases, I modified computation of I2 with the following:

I2x = scipy.integrate.cumulative_simpson(p2 * dr2dx, x=x, axis=1, initial=0)
for j in range(ny):
    I2x[j,:] -= I2x[j,-1]    # to set integration constant at each x
I2y = scipy.integrate.cumulative_simpson(p2 * dr2dy, x=y, axis=0, initial=0)
for i in range(nx):
    I2y[:,i] -= I2y[-1,i]    # to set integration constant at each y
I2 = I2x + I2y
I2 = I2 * (I2 < 0)    # removes some numerical residues

I also tried to do the same in 3d and the answer tripled. So it probably has something to do with dr. Is it wrong to write dr = (∂r/∂x) dx + (∂r/∂y) dy?

本文标签： pythonMissing factor of 2 when numerically integrating a 1d function on 2d gridStack Overflow