计算科学 - 对预期的圆柱对称势场的梯度进行插值遵循拉普拉斯方程，尤其是在 r=0 附近/跨过 r=0 - 吾爱随笔录

下面的脚本尝试为静电透镜实现势场的Jacobi 迭代松弛。

为了绘制电场线并计算带电粒子的轨迹，我需要编写一个函数来计算网格点之间指定点处的插值梯度。我可以找到适合最接近兴趣点的 3x3 数组的二次函数的系数，然后计算和的导数，但这让我很担心，因为： $f(x, y) = ax^2 + bxy + cy^2 + dx + ey + f$ $x$ $y$

它没有利用该领域应该是拉普拉斯方程的一个很好的解决方案这一事实
它不“知道”这是柱坐标中拉普拉斯方程的解；例如，其中没有像松弛（术语） $1/r$ one_over_8ir

如果/当粒子轨迹穿过对称轴时，第二个变得越来越重要，并且插值器应该尊重这种对称性。

注意：我将在 3D 中计算轨迹。它们将在、和方向上有一个速度分量，尽管我认为这对于目的而言并不重要这个问题的。 $\mathbf{\hat{z}}$ $\mathbf{\hat{r}}$ $\mathbf{\hat{\theta}}$

问题：有什么好的二阶（或者可能是三阶）方法可以在柱坐标中插值我数值计算的势场梯度，尤其是在 r=0 附近/跨过 r=0。

静电透镜三联体的势场

Python脚本：

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