Interpolating the gradient of a cylindrically symmetric potential field expected obey the Laplace equation, especially near/across r=0

Question

The script below tries to implement a Jacobi iterative relaxation of a potential field for an electrostatic lens.

In order to plot electric field lines and calculate trajectories for charged particles, I need to write a function that calculates an interpolated gradient at a specified point between grid points. I can find the coefficients of a quadratic function $f(x, y) = ax^2 + bxy + cy^2 + dx + ey + f$ that fits the 3x3 array nearest the point of interest then calculate the derivatives in $x$ and $y$, but that worries me because:

1. it doesn't leverage the fact that the field should be a pretty good solution of the Laplace equation
2. it doesn't "know" that this is a solution to the Laplace equation in cylindrical coordinates; for example, there's no $1/r$-like term in it like there is in the relaxation (the one_over_8ir term).

That second one becomes increasingly important if/when particle trajectories cross the symmetry axis, and the interpolator should honor that symmetry.

note: I'll be calculating trajectories in 3D. They will have a velocity component in $\mathbf{\hat{z}}$, $\mathbf{\hat{r}}$ and $\mathbf{\hat{\theta}}$ directions, though I don't think this matters for the purposes of this question.

Question: What's a good 2nd (or perhaps 3rd) order method to interpolate the gradient of my numerically calculated potential field in cylindrical coordinates, especially near/across r=0.

Python script:

class Ring(object):
    def __init__(self, iz0, iz1, ir, phi):
        self.ir = ir
        self.iz0 = iz0
        self.iz1 = iz1
        self.phi = float(phi)

def do_it(N):
    # https://math.stackexchange.com/questions/2067439/numerical-solution-to-laplace-equation-using-a-centred-difference-approach-in-cy
    # phi (i, k) = (1/4) * ((i+1, k) + (i-1, k) + (i, k+1) + (i, ki1)) + (1/8k) * ((i+1, k) - (i-1, k))  # r > 0
    # phi (i, 0) = (2/3) * (i, 1) + (1/6) * ((i+1, 0) + (i-1, 0))  # r = 0
    for i in range(N):
        phi_new = (quarter      * (np.roll(phi, -1, axis=0) + np.roll(phi, +1, axis=0) +
                                   np.roll(phi, -1, axis=1) + np.roll(phi, +1, axis=1)) +
                   one_over_8ir * (np.roll(phi, +1, axis=0) - (np.roll(phi, +1, axis=0))))

        phi_new[:, 0]  = (two_thirds * phi[:, 1] +
                          one_sixth  * (np.roll(phi[:, 0], -1, axis=0) + np.roll(phi[:, 0], +1, axis=0)))

        phi[do_me]     = phi_new[do_me]
                                 
# see also https://scicomp.stackexchange.com/questions/30839/python-finite-difference-schemes-for-1d-heat-equation-how-to-express-for-loop-u

import numpy as np
import matplotlib.pyplot as plt

from scipy.integrate import odeint as ODEint

nz, nr = 60, 20

quarter, two_thirds, one_sixth = 1./4, 2./3, 1./6
one_over_8ir = np.hstack(([0.0], 1./(8*np.arange(1, nr, dtype=float))))[None, :]

lens_def = ((6, 13, 16, -1), (19, 41, 16, +1), (47, 54, 16, -1))

rings = []
for thing in lens_def:
    ring = Ring(*thing)
    rings.append(ring)

phi0  = np.zeros((nz, nr)) # zero everythwere

do_me = np.ones_like(phi0, dtype = bool)

do_me[ 0,  :] = False
do_me[-1,  :] = False
do_me[ :, -1] = False

for ring in rings:
    do_me[ring.iz0:ring.iz1+1, ring.ir] = False
    phi0[ ring.iz0:ring.iz1+1, ring.ir] = ring.phi

phi = phi0.copy()

do_it(10000)

phit = phi.T.copy()

phii = np.vstack((phit[1:][::-1], phit))

if True:
    plt.figure()
    vv = np.abs(phii).max()
    plt.imshow(phii, vmin=-vv, vmax=+vv, cmap='RdBu',
               origin='lower', interpolation='nearest')
    plt.colorbar()
    plt.show()

if False:
    plt.figure()
    for (i, phi) in enumerate(phiz):
        vv = np.abs(phi).max()
        plt.subplot(2, 2, i+1)
        plt.imshow(phi.T, vmin=-vv, vmax=+vv, cmap='RdBu',
                   origin='lower', interpolation='nearest')
        plt.colorbar()
    plt.show()

Answer 1

Provided that your geometry conforms to the cylindrical coordinate system, the separation of variables solution should look something like $$ \Phi(\mathbf{r}) = \sum_{n,m} a_{nm} \left\lbrace \begin{split}J_n(k_{nm} r)\\ N_n(k_{nm}r)\end{split} \right\rbrace e^{in\theta} e^{\pm k_{nm}z}. $$ You can figure out which functions are appropriate from your boundary conditions. You can find the coefficients $a_{nm}$ by setting your solution as a function of $\mathbf{r}$ equal to the sum on the right-hand side above and numerically integrating your solution against each eigenfunction (you should exploit the orthogonality of complex exponentials and Bessel functions).

Once you have the coefficients, you can just evaluate the first sum (or its gradient) above at any point $(r,\theta,z)$ in your solution domain to obtain your desired quantity.