I am working on a loop to calculate the elements of a 3D numpy array. The code looks like this:
A = np.zeros((N - 2, N - 2, N - 2))
for i in range(1, N - 1):
for j in range(1, N - 1):
for k in range(1, N - 1):
A[i - 1, j - 1, k - 1] = (B[i,j,k] * dy[j] * dz[k] + B[i,j-1,k] * dy[j-1] * dz[k]
+ B[i,j,k-1] * dy[j] * dz[k-1] + B[i,j-1,k-1] * dy[j-1] * dz[k-1]) / (4 * dx[i])
The above method is becoming slow for large values of N, so I attempted a different approach:
A = (B[1:-1, 1:-1, 1:-1] * dy[1:] * dz[1:] + B[1:-1, :-2, 1:-1] * dy[:-1] * dz[1:]
+ B[1:-1, 1:-1, :-2] * dy[1:] * dz[:-1] + B[1:-1, :-2, :-2] * dy[:-1] * dz[:-1]) / (4 * dx[1:])
However, the alternate method fails when dx, dy, and dz are not constant. I also experimented with the following approach:
dX, dY, dZ = np.meshgrid(dx, dy, dz)
A = (B[1:-1, 1:-1, 1:-1] * dY[1:, 1:, 1:] * dZ[1:, 1:, 1:] + B[1:-1, :-2, 1:-1] * dY[:-1, :-1, :-1] * dZ[1:, 1:, 1:]
+ B[1:-1, 1:-1, :-2] * dY[1:, 1:, 1:] * dZ[:-1, :-1, :-1] + B[1:-1, :-2, :-2] * dY[:-1, :-1, :-1] * dZ[:-1, :-1, :-1]) / (4 * dX[:-1, :-1, :-1])
Unfortunately, this approach also failed to produce the desired result.
If anyone has any ideas or suggestions on how to improve this calculation, I would greatly appreciate it.