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:
- it doesn't leverage the fact that the field should be a pretty good solution of the Laplace equation
- 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()