# Interpolating the gradient of a cylindrically symmetric potential field that's 'supposed to' obey the Laplace equation?

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

It's hot-off-the-press and I've just started to debug and look for things to test it against.

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? I don't need script, just a short explanation, with a link if necessary. 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]

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()

• Use local Lagrange interpolation using piecewise Lagrange polynomials to interpolate the solution and use their derivatives to calculate the gradient: en.wikipedia.org/wiki/Lagrange_polynomial – smh Jun 3 '19 at 15:36
• @smh Do Lagrange polynomials directly address either of the two numbered points in the question? – uhoh Jun 3 '19 at 16:10
• No, but if you're going to be calculating forces acting on moving particles, it's going to be quite expensive to project the solution onto the eigenmodes of the Laplace operator and evaluate their gradients at each particle location at each time step. A carefully implemented Lagrange interpolator should work just fine. – smh Jun 3 '19 at 17:49
• @smh I've done this in the past; for what I need to do speed won't be a problem at all. This time however I need to focus (no pun intended) on accuracy, thus it's time to take a closer look at the interpolation. Perhaps you would like to expand on "project the solution onto the eigenmodes of the Laplace operator" in an answer? – uhoh Jun 3 '19 at 17:53

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.
• Okay I see what you mean. My field is cylindrically symmetric so only $n=0$ terms are used? I was imagining a polynomial-like interpolator of local grid points only, not requiring calls to transcendentals for each evaluation, but this looks really interesting so I'll give it a try today. Thanks! – uhoh Jun 3 '19 at 22:04