How to deal with numerical errors in electrostatic field calculations
I want to trace electrostatic field lines emerging from 2D surfaces in 3D space. Eventually I want to find their intersection with an (uncharged) mesh.
The charge distribution $\sigma(x), x \in \Gamma \subset \mathbb{R}^3$ is defined by requiring a constant potential $\phi(x) = \phi_\Gamma\ \forall x\in\Gamma$.
For the potential field $\phi$ and the electrostatic field $E$ we have:
$$ \phi(x) = \int_\Gamma \sigma(y) \frac{1}{\|x-y\|} dy,\quad x\in\mathbb{R}^3 \ E(x) = \nabla \phi(x) = \int_\Gamma \sigma(y) \frac{x-y}{\|x-y\|^3} dy,\quad x\in\mathbb{R}^3 $$
The potential function $\phi$ is harmonic ($\Delta\phi=0$), so the points on the charged surfaces are the maxima of $\phi$ and the electric field is $0$ there. To trace the field lines I place two particles with a small offset in/against normal direction. The particles have in general a non-zero $E$ field, which can be used to trace the field lines with a time stepping scheme
$$ x_{i+1} = x_i + E(x_i) $$
For a surface like a plane you get two field lines in opposite directions and the field vanishes at infinity. 402
When the surface belongs to a solid object, e.g a sphere, the potential field is constant and the electric field is $0$ on the inside, because $\phi$ is harmonic.
To calculate the field line, I use an approximation $\tilde{E}$ by replacing the integral over $\Gamma$ with a finite sum over triangles $T_i$:
$$ E(x) \approx \sum_j q_j \int_{T_j} \frac{1}{\|x-y\|} dy \approx \sum_j q_j \sum_k w_k \frac{1}{\|x-y_k\|} = \tilde{E}(x) $$
where $q_j = \int_{T_j} \sigma(y) dy$ and the last sum is a gauss quadrature of the integral over the triangle $T_j$. The potential $\tilde{\phi}$ can be approximated in the same way.
Now the problem is that as the functions are approximations, we have on the inside $\Omega \subset \mathbb{R}^3$ of the solid object
$$ \|\tilde{E}(x)\| > 0 \ \tilde{\phi}(x) \not\equiv \text{const}, x \in \Omega $$
and the field $\tilde{E}$ often points outside the object.
While $\|\tilde{E}\|$ inside the object is typically small, the field lines still eventually cross the surface at a step $x_i \to x_{i+1}$.
The field has a discontinuity on the surface, so $\|E\|$ is zero on the inside and large on the outside near the surface. When tracing a particle in the approximated field, you have near the surface $\tilde{E}(x_i) \ll \tilde{E}(x_{i+1}), x_i \in \Omega, \ x_{i+1} \not\in \Omega$.
So the error between the position $x_k$ on an analytic field line and the position $\tilde{x}_k$ is small for $k \le i$ and huge for $k > i$.
The problem does not only affect solids, but creates local extrema in the space between non-closed charged surfaces as well, which leads to wrong field lines.
I thought of different ways to fix this numeric problems, but all ways seem to affect valid field lines as well:
* Rounding small values to $0$ can fix the field inside solids. But this will make local minima even worse.
* Disallowing field lines to cross surfaces. This fixes the solid problem but not the local extrema.
* Adding some inertia term could fix field lines getting stuck in local minima, but will worsen the problem of the non-zero field inside solids. 903