# 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. 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.
• Normalizing $E$ to $E/\|E\|$. This makes the problem much worse for errors inside solid objects and requires choosing a smaller timestep size for all steps.

The approximations, which influence the accurracy of the solution:

• The the resolution of the triangulation. Can be refined if needed, but makes the calculation much more costly as calculating the field at a single point is in $\mathcal{O}(n)$. So calculating one timestep of the field lines is in $\mathcal{O}(n^2)$.
• The approximation of the charges. I have a piecewise constant charge on each triangle and currently no way to approximate it with higher order.
• The approximation of the integrals in $\int_{T_j} q_j \frac{x-y}{\|x-y\|^3}dy$. I currently use 6th order gauss quadrature.
• The calculation of the charges $q_j$ by requiring $\phi(x)=\phi_\Gamma \ \forall x\in\Gamma$ is approximated by requiring $\phi(x_i)=\phi_\Gamma$, where $x_i \in T_i$ are the centers of the triangles.
• What is your question? Aug 24 '17 at 21:23
• I made a bigger edit. I hope its more clear now. My main problem is not the actual error in the field (which is expected when using an approximation), but handling how the small error leads to a very big error when particles cross the field discontinuity on the surface.
– allo
Aug 25 '17 at 7:53
• I think that what you are looking looking is similar to a boundary integral formulation. And the numerical counterpart would be a boundary element method. Aug 28 '17 at 4:03
• This would mean solving for E at position $x\in R^3$ by discretizing the integral over each triangle by using BEM on the edges of the triangle? If possible I was looking for something easier, especially because I got the charge value piecewise constant on the triangles only.
– allo
Aug 29 '17 at 13:39
• No, it means solving in $R^3$ using the values in your triangular mesh, for example. Intuitively, you could think that you can concentrate the charge of your triangle in each barycenter, and then just superimpose the electric field (Green function) of each one. That will give you the total electric field Aug 29 '17 at 13:51

I will try to answer the part of the question regarding the accuracy of the calculation, which certainly affects all the other things.

Integrals of the type ($T_j$ denoting the $j$th triangle in the discretization and $\vec{x}$ some arbitrary observation point)

$$\int\limits_{T_j}\frac{1}{||\vec{x}-\vec{y}||}dy,\quad \vec{y}\in T_j$$

have a singularity at $\vec{x}=\vec{y}$ which implies that a special care should be taken when $\vec{x}$ is close to $\vec{y}$. As I understood, your special treatment is limited to choosing even/odd Gauss quadrature rules, so that $\vec{x}$ and $\vec{y}$ never collide exactly. While such ad-hoc method will provide some result, it is avoiding the problem of handling the singularity rather than solving it. Moreover, increasing accuracy using this approach is nearly impossible; thus, another technique should be used.

There are special quadrature rules designed to handle such singularities:

But, I would propose handling such singularity analytically. The following paper:

describes how to take those integrals analytically for planar polygons (certainly works for triangles). The theory is not too complicated and the implementation should be faster than using quadrature rules with a large number of quadrature points.

This is a must when $\vec{x}$ and $\vec{y}$ can lie in the same triangle, but is recommended for close enough interactions. You will choose your singularity treatment region yourself based on accuracy consideration. A good rule of thumb is to choose

$$d_\text{sing. treatment}=4\max_\limits{j}{\max_\limits{s=1,2,3}e_{j,s}},$$

where $e_{j,s}$ is the length of the $s$th edge of the $j$th triangle. Then, if $||\vec{x}-\vec{y}||<d_\text{sing. treatment}$, you apply the singularity treatment of your choice (analytic or special quadrature rule); otherwise, proceed with regular cheap Gaussian quadrature.

I have also discussed a bit this paper and the technique in this question for the case of a non-static integral (with $e^{jkR}$).

Using a proper singularity treatment technique should increase the accuracy of your results dramatically which might eliminate a lot of problems with local errors and field lines near the interfaces.

• Thank you, this looks like a possible solution and I accept this answer now. I will probably not be able to test it soon, but I think it will help when I am working at the corresponding code again.
– allo
Jun 19 '18 at 13:29