Reverse gradient: Integrate two 2d arrays

Question

I have two 2d arrays of electric field components:

E_x with (256, 512) and E_z with (256, 512) sizes.

I need to calculate scalar potential. Gradient of scalar potential is the electric field vector. So I need the reverse of the gradient.

I tried using the following loop to calculate the potential Φ (phi):

nx = 256
nz = 512
phi = np.zeros([nx, nz])

for i in range(1, nz):
    phi[:, i] = -dx * Ex[:, i-1] + phi[:, i-1]
    for j in range(1, nx):
        phi[j, :] = -dz * Ez[j-1, :] + phi[j-1, :]

This gives me an ok but weird Φ (these vertical lines are suspicious):

I expected to receive something like this (a colleague calculated the same value in Fortran and his profile is different):

This is their Fortran code:

allocate(phi(nx, nz))

phi(1, 1) = 0.0

do i = 2, nx
    phi(i, 1) = -dx * ex(i-1, 1) + phi(i-1, 1)
enddo

do i = 1, nx
    do j = 2, nz
        phi(i, j) = -dz * ez(i, j-1) + phi(i, j-1)
    enddo
enddo

Where did I make a mistake? Or is there a better way to calculate the scalar potential from an electric field?

Answer 1

What you are doing now is integrating the first row and then all the columns one by one (as shown on the top of the image below). The problem with this approach is that the accumulation of the error gets greater and greater the further you are from the starting point and the columns get dissociated because they don't have the same error. You can also see that you basically only use half of your data (only one component of the gradient).

To solve this you could do as suggested on the bottom of the image : integrating by going diagonally. Like that, when you are not on the border of the array, you can use two already calculated values (on the left and below) instead of one and taking the average of those should add some cohesion between the columns. If I'm correct, it would even divide the error by two.

I read your code again and I don't know how I didn't see it the first time but you are using slices. I think that's actually were your problem comes from since because of that your code isn't equivalent to your colleague's one. Your are not affecting one value at the time but whole rows/columns. Another thing I noticed is that you used Ex and dx while moving along the z axis and vice-versa.

Try this modified version of your code:

nx = 256
nz = 512
phi = np.zeros([nx, nz])

for i in range(1, nx):
    phi[i, 0] = -dx * Ex[i-1, 0] + phi[i-1, 0]
    for j in range(1, nz):
        phi[i, j] = -dz * Ez[i, j-1] + phi[i, j-1]