Electrostatic Simulator: Mathematical & Physical Foundations

1. Physics Foundations

Governing laws, fields, potentials, and fundamental physical constants

1.1 Coulomb's Law

Coulomb's Law quantifies the electrostatic force $\mathbf{F}_{12}$ exerted by a point charge $q_1$ located at $\mathbf{r}_1$ upon a second charge $q_2$ located at $\mathbf{r}_2$ in a vacuum.

$$\mathbf{F}_{12} = \frac{1}{4\pi\varepsilon_0} \frac{q_1 q_2}{|\mathbf{r}_2 - \mathbf{r}_1|^2} \hat{\mathbf{r}}_{12} = \frac{1}{4\pi\varepsilon_0} \frac{q_1 q_2 (\mathbf{r}_2 - \mathbf{r}_1)}{|\mathbf{r}_2 - \mathbf{r}_1|^3}$$

Here, $\varepsilon_0$ denotes the permittivity of free space, defining the capability of a vacuum to permit electric field line formation. The force vector obeys Newton's third law ($\mathbf{F}_{12} = -\mathbf{F}_{21}$) and follows an inverse-square distance dependence ($r^{-2}$).

1.2 Electric Field ($\mathbf{E}$) and Electric Potential ($V$)

The electric field $\mathbf{E}(\mathbf{r})$ is a vector field representing force per unit test charge $q_0$. For a continuous spatial charge distribution with volumetric charge density $\rho(\mathbf{r}')$:

$$\mathbf{E}(\mathbf{r}) = \frac{1}{4\pi\varepsilon_0} \int_V \frac{\rho(\mathbf{r}') (\mathbf{r} - \mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|^3} \, d^3\mathbf{r}'$$

Because static electric fields are conservative ($\nabla \times \mathbf{E} = \mathbf{0}$), $\mathbf{E}$ can be uniquely expressed as the negative gradient of a scalar Electric Potential $V(\mathbf{r})$:

$$\mathbf{E} = -\nabla V \quad \iff \quad V(\mathbf{r}_b) - V(\mathbf{r}_a) = -\int_{\mathbf{r}_a}^{\mathbf{r}_b} \mathbf{E} \cdot d\mathbf{l}$$

Physical meaning: Potential $V$ measures work required per unit charge to move from a reference point (typically infinity, where $V(\infty)=0$) to location $\mathbf{r}$. The negative sign dictates that the field points directionally along the steepest decrease in potential energy.

1.3 Related Electrostatic Maxwell Equations

Electrostatics is governed by two fundamental Maxwell equations in static limit ($\frac{\partial \mathbf{B}}{\partial t} = \mathbf{0}$):

Gauss's Law

Differential Form: $$\nabla \cdot \mathbf{E} = \frac{\rho}{\varepsilon_0}$$

Integral Form: $$\oint_{\partial \Omega} \mathbf{E} \cdot d\mathbf{A} = \frac{Q_{\text{enclosed}}}{\varepsilon_0}$$

Electrostatic Condition

Differential Form: $$\nabla \times \mathbf{E} = \mathbf{0}$$

Integral Form: $$\oint_{\mathcal{C}} \mathbf{E} \cdot d\mathbf{l} = 0$$

Fundamental Physical Constants

Vacuum Permittivity ($\varepsilon_0$) $\varepsilon_0 \approx 8.854187 \times 10^{-12} \text{ F/m}$

Coulomb Constant ($k_e = \frac{1}{4\pi\varepsilon_0}$) $k_e \approx 8.987551 \times 10^{9} \text{ N}\cdot\text{m}^2/\text{C}^2$

Concept of Scalar Heat Maps

Numerical simulators represent the continuous scalar field $V(x,y)$ on discrete 2D spatial grids. Mapping raw numerical potential values to visual spectrums provides intuitive field visualization:

Positive High Potential ($V \gg 0$) $\rightarrow$ Warm Red/Orange

Neutral Ground Zero ($V = 0$) $\rightarrow$ Dark Charcoal/Black

Negative Low Potential ($V \ll 0$) $\rightarrow$ Cool Blue/Cyan

2. Vector Calculus

Differential operators, gradients, field flux divergence, and irrotational conditions

Scalar Field $\rightarrow$ Vector Field

Gradient ($\nabla V$)

The gradient operator $\nabla$ operating on scalar potential $V(x,y,z)$ produces a vector field pointing toward the maximal rate of spatial increase:

$$\nabla V = \frac{\partial V}{\partial x}\hat{\mathbf{i}} + \frac{\partial V}{\partial y}\hat{\mathbf{j}} + \frac{\partial V}{\partial z}\hat{\mathbf{k}}$$

Directional derivative along unit vector $\hat{\mathbf{u}}$ is $D_{\hat{\mathbf{u}}}V = \nabla V \cdot \hat{\mathbf{u}}$. Thus $\mathbf{E} = -\nabla V$ naturally directs charge flow down steepest potential gradients.

Vector Field $\rightarrow$ Scalar Field

Divergence ($\nabla \cdot \mathbf{E}$)

Computes net outward vector flux per unit volume at an infinitesimal volume element:

$$\nabla \cdot \mathbf{E} = \frac{\partial E_x}{\partial x} + \frac{\partial E_y}{\partial y} + \frac{\partial E_z}{\partial z}$$

By Divergence Theorem, $\iiint_V (\nabla \cdot \mathbf{E}) dV = \oiint_{\partial V} \mathbf{E} \cdot d\mathbf{A}$. $\nabla \cdot \mathbf{E} > 0$ denotes positive charge sources ($\rho > 0$).

Vector Field $\rightarrow$ Vector Field

Curl ($\nabla \times \mathbf{A}$)

Measures localized circulation/rotation density of vector field $\mathbf{A}$:

$$\nabla \times \mathbf{A} = \begin{vmatrix} \hat{\mathbf{i}} & \hat{\mathbf{j}} & \hat{\mathbf{k}} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ A_x & A_y & A_z \end{vmatrix}$$

Electrostatic Proof: Since $\mathbf{E} = -\nabla V$, then $\nabla \times \mathbf{E} = -\nabla \times (\nabla V) \equiv \mathbf{0}$ identically (the curl of any gradient field is identically zero).

3. Numerical Methods & Discretization

Taylor expansion finite differences, 5-point Poisson stencil, and iterative solver schemes

3.1 First & Second-Degree Finite Differences

Consider a 1D grid with spatial step spacing $h = \Delta x$. The first-order spatial derivatives of potential $V(x)$ are approximated via forward, backward, and central difference stencils:

Forward Difference $\mathcal{O}(h)$

$$\frac{\partial V}{\partial x} \approx \frac{V_{i+1,j} - V_{i,j}}{h}$$

Backward Difference $\mathcal{O}(h)$

$$\frac{\partial V}{\partial x} \approx \frac{V_{i,j} - V_{i-1,j}}{h}$$

Central Difference $\mathcal{O}(h^2)$

$$\frac{\partial V}{\partial x} \approx \frac{V_{i+1,j} - V_{i-1,j}}{2h}$$

Step-by-Step Derivation of Second-Degree Central Derivative

To derive the second-degree partial derivative $\frac{\partial^2 V}{\partial x^2}$, expand $V(x+h, y)$ and $V(x-h, y)$ using Taylor series expansions around node $(x,y)$:

(Eq 1 - Forward): $$V(x+h, y) = V(x,y) + h \frac{\partial V}{\partial x} + \frac{h^2}{2!} \frac{\partial^2 V}{\partial x^2} + \frac{h^3}{3!} \frac{\partial^3 V}{\partial x^3} + \frac{h^4}{4!} \frac{\partial^4 V}{\partial x^4} + \mathcal{O}(h^5)$$

(Eq 2 - Backward): $$V(x-h, y) = V(x,y) - h \frac{\partial V}{\partial x} + \frac{h^2}{2!} \frac{\partial^2 V}{\partial x^2} - \frac{h^3}{3!} \frac{\partial^3 V}{\partial x^3} + \frac{h^4}{4!} \frac{\partial^4 V}{\partial x^4} - \mathcal{O}(h^5)$$

Adding Equation (1) and Equation (2) cancels all odd-order derivative terms ($h, h^3, \dots$):

$$V(x+h, y) + V(x-h, y) = 2V(x,y) + h^2 \frac{\partial^2 V}{\partial x^2} + \mathcal{O}(h^4)$$

Isolating $\frac{\partial^2 V}{\partial x^2}$ and switching to discrete grid indexing $(i, j)$ yields the second-order accurate central difference formula:

$$\frac{\partial^2 V}{\partial x^2} \Bigg|_{i,j} = \frac{V_{i+1,j} - 2V_{i,j} + V_{i-1,j}}{h^2} + \mathcal{O}(h^2)$$

3.2 Derivation of the 5-Point Stencil for Poisson's Equation

Combining Gauss's Law ($\nabla \cdot \mathbf{E} = \frac{\rho}{\varepsilon_0}$) with conservative gradient relation ($\mathbf{E} = -\nabla V$) yields Poisson's Equation:

$$\nabla \cdot (-\nabla V) = -\nabla^2 V = \frac{\rho}{\varepsilon_0} \quad \implies \quad \nabla^2 V(x,y) = -\frac{\rho(x,y)}{\varepsilon_0}$$

Expand the 2D Laplacian operator $\nabla^2 V = \frac{\partial^2 V}{\partial x^2} + \frac{\partial^2 V}{\partial y^2}$ by substituting second-order central difference formulas along both grid dimensions $x$ and $y$ with uniform mesh size $h = \Delta x = \Delta y$:

$$\frac{V_{i+1,j} - 2V_{i,j} + V_{i-1,j}}{h^2} + \frac{V_{i,j+1} - 2V_{i,j} + V_{i,j-1}}{h^2} = -\frac{\rho_{i,j}}{\varepsilon_0}$$

Group terms with common denominator $h^2$:

$$\frac{V_{i+1,j} + V_{i-1,j} + V_{i,j+1} + V_{i,j-1} - 4V_{i,j}}{h^2} = -\frac{\rho_{i,j}}{\varepsilon_0}$$

Multiply by $h^2$, rearrange terms, and solve for node potential $V_{i,j}$ to arrive at the explicit discrete update formula:

$$V_{i,j} = \frac{1}{4} \left( V_{i+1,j} + V_{i-1,j} + V_{i,j+1} + V_{i,j-1} + \frac{h^2 \rho_{i,j}}{\varepsilon_0} \right)$$

5-Point Stencil Weights

NORTH$+1/4$

WEST$+1/4$

CENTER$V_{i,j}$

EAST$+1/4$

SOUTH$+1/4$

Charge Source: $+ \frac{h^2 \rho_{i,j}}{4\varepsilon_0}$

3.3 Iterative Solvers for Poisson's Grid Matrix

Jacobi Method

Simultaneous

Computes step $(k+1)$ strictly using previous iteration values $k$. Requires keeping dual array buffers in memory.

$$V_{i,j}^{(k+1)} = \frac{1}{4}\left(V_{i+1,j}^{(k)} + V_{i-1,j}^{(k)} + V_{i,j+1}^{(k)} + V_{i,j-1}^{(k)} + \frac{h^2 \rho_{i,j}}{\varepsilon_0}\right)$$

Slow convergence rate ($\mathcal{O}(N^2)$ iterations)

Gauss-Seidel

In-Place

Updates array values in-place. Uses newly computed $(k+1)$ values for calculated neighbors $(i-1,j)$ and $(i,j-1)$ immediately.

$$V_{i,j}^{(k+1)} = \frac{1}{4}\left(V_{i+1,j}^{(k)} + V_{i-1,j}^{(k+1)} + V_{i,j+1}^{(k)} + V_{i,j-1}^{(k+1)} + \frac{h^2 \rho_{i,j}}{\varepsilon_0}\right)$$

Single buffer array (~2x faster convergence)

SOR Method

Accelerated

Over-relaxes Gauss-Seidel step via relaxation parameter $\omega \in (1, 2)$ to dramatically accelerate convergence:

$$V_{i,j}^{(k+1)} = (1-\omega)V_{i,j}^{(k)} + \omega V_{i,j}^{\text{GS}}$$

Optimal $\omega = \frac{2}{1 + \sin(\pi/N)}$ reduces iterations to $\mathcal{O}(N)$

4. Connecting the Dots: Complete Pipeline

Mapping continuous electrodynamics to computational grid routines

The block diagram below outlines the full algorithmic flow executed inside the electrostatic engine:

1

Spatial Grid & Boundary Initialization

Discretize spatial domain $\Omega$ into $N \times N$ matrix mesh nodes $V[N][N]$. Set boundary condition potentials (e.g. $V_{\text{boundary}} = 0$ for grounded walls or Dirichlet conductor plates).

2

Charge Density Mapping ($\rho$)

Map discrete point charges $q_k$ onto neighboring spatial mesh points $(i, j)$ to build source charge matrix $\rho[i][j]$.

3

Iterative Poisson Stencil Solving (SOR)

Loop over internal non-boundary grid points applying 5-point Poisson SOR stencil updates until residual norm drops below threshold $\epsilon$:

$$\max_{i,j} \left| V_{i,j}^{(k+1)} - V_{i,j}^{(k)} \right| < 10^{-5}$$

4

Electric Vector Field Gradient Calculation ($\mathbf{E} = -\nabla V$)

Compute discrete vector components $E_x(i,j)$ and $E_y(i,j)$ across converged potential array $V$ using central differences:

$$E_x(i,j) = -\frac{V_{i+1,j} - V_{i-1,j}}{2h}, \quad E_y(i,j) = -\frac{V_{i,j+1} - V_{i,j-1}}{2h}$$

5

Color Heatmap Rendering & Vector Quiver Overlay

Map potential array $V(i,j)$ into RGB color spectrum buffer while overlaying scaled directional arrows for vector field $\mathbf{E}(i,j)$.

Explore Code Implementation

View complete source code, matrix algorithms, and build instructions on GitHub.

anuz-912/electrostatic-simulator

Mathematical & Physical Foundations of Electrostatic Simulation