Electric Field Potential and the Relaxation Method

1. Laplace and Poisson Equations

Electric potential φ is related to the field E = −∇φ. From Gauss's law ∇·E = ρ/ε₀ we get:

Laplace's equation describes the potential in a charge-free region — between conductors with fixed potentials. Poisson's equation is the general case (charges present). Both are of elliptic type: the solution is "smooth" and depends on boundary conditions across the entire boundary.

Physical Interpretation

2. Finite-Difference Discretisation

We divide the domain into a uniform 2D grid with step h. The Laplacian is approximated by the five-point stencil:

Substituting into Poisson's equation ∇²φ = f (right-hand side = −ρ/ε): we get a linear system A·φ = b, where A is a sparse matrix (5 non-zero elements per row in 2D). Dimension N × N for an N×N grid. Direct solution (LU) is expensive — O(N³). Iterative methods — O(N²) per iteration, converging in O(N) – O(N log N) iterations.

3. Jacobi Method

The simplest iterative method: express φ_{i,j} through its neighbours using old values:

Convergence: spectral radius ρ_J ≈ 1 − π²/(2N²) for a square N×N grid. Iterations for error ε: k ≈ N² ln(1/ε) / π² ∝ N². For N=100 and ε=10⁻⁶ — ~28 000 iterations. Slow!

Why Jacobi Maps Naturally to WebGL

Each texture pixel = one grid node. Jacobi is perfectly parallel — new values are independent of each other. One pass = one draw call in a ping-pong framebuffer.

// GLSL fragment shader — one Jacobi iteration
precision highp float;
uniform sampler2D uPhi;   // current φ approximation
uniform sampler2D uRho;   // charge density ρ (normalised)
uniform vec2 uTexelSize;  // 1/width, 1/height
uniform float uH2;        // h² (grid step squared)
varying vec2 vUV;

void main() {
  float n = texture2D(uPhi, vUV + vec2(0., uTexelSize.y)).r;
  float s = texture2D(uPhi, vUV - vec2(0., uTexelSize.y)).r;
  float e = texture2D(uPhi, vUV + vec2(uTexelSize.x, 0.)).r;
  float w = texture2D(uPhi, vUV - vec2(uTexelSize.x, 0.)).r;
  float rho = texture2D(uRho, vUV).r;
  gl_FragColor = vec4((n + s + e + w + uH2 * rho) / 4.0, 0., 0., 1.);
}

4. Gauss-Seidel Method

Unlike Jacobi: immediately use new neighbour values as soon as they are computed (left and lower neighbours are already updated):

Convergence is twice as fast as Jacobi: ρ_GS ≈ ρ_J². In practice, half the iterations for the same error. GPU challenge: pixel dependencies — requires a checkerboard update order ("red-black Gauss-Seidel").

Red-Black Ordering

Split the grid into "red" (i+j even) and "black" (i+j odd) nodes. Red nodes depend only on black and vice versa. Two passes per iteration:

function redBlackGS(phi, rho, N, h2, iterations) {
  for (let it = 0; it < iterations; it++) {
    for (let color = 0; color < 2; color++) {  // 0=red, 1=black
      for (let i = 1; i < N - 1; i++) {
        for (let j = 1; j < N - 1; j++) {
          if ((i + j) % 2 !== color) continue;
          const idx = i * N + j;
          phi[idx] = (
            phi[(i+1) * N + j] + phi[(i-1) * N + j] +
            phi[i * N + j+1] + phi[i * N + j-1] +
            h2 * rho[idx]
          ) / 4;
        }
      }
    }
  }
}

5. Successive Over-Relaxation (SOR)

SOR (Successive Over-Relaxation) improves Gauss-Seidel by introducing a relaxation parameter ω ∈ (1, 2):

With ω = ω_opt the iteration count drops from O(N²) to O(N) — a fundamentally different order! For N=100 this means ~28 000 → ~300 iterations.

class PoissonSolverSOR {
  constructor(N, h = 1.0 / N) {
    this.N = N;
    this.h2 = h * h;
    this.omega = 2 / (1 + Math.sin(Math.PI / N)); // ω_opt
    this.phi = new Float64Array(N * N);  // potential
    this.rho = new Float64Array(N * N);  // ρ/ε₀
    this.fixed = new Uint8Array(N * N);  // Dirichlet BC mask
  }

  setCharge(i, j, value) { this.rho[i * this.N + j] = -value; }
  setFixed(i, j, value) {
    const idx = i * this.N + j;
    this.phi[idx] = value;
    this.fixed[idx] = 1;
  }

  step(iterations = 1) {
    const { N, h2, omega, phi, rho, fixed } = this;
    for (let it = 0; it < iterations; it++) {
      for (let color = 0; color < 2; color++) {
        for (let i = 1; i < N - 1; i++) {
          for (let j = 1; j < N - 1; j++) {
            if ((i + j) % 2 !== color) continue;
            const idx = i * N + j;
            if (fixed[idx]) continue;
            const gs = (
              phi[(i+1)*N+j] + phi[(i-1)*N+j] +
              phi[i*N+j+1] + phi[i*N+j-1] - h2 * rho[idx]
            ) / 4;
            phi[idx] += omega * (gs - phi[idx]);
          }
        }
      }
    }
  }

  residual() {
    const { N, h2, phi, rho, fixed } = this;
    let max = 0;
    for (let i = 1; i < N - 1; i++) {
      for (let j = 1; j < N - 1; j++) {
        const idx = i * N + j;
        if (fixed[idx]) continue;
        const lap = (
          phi[(i+1)*N+j] + phi[(i-1)*N+j] +
          phi[i*N+j+1] + phi[i*N+j-1] - 4*phi[idx]
        ) / h2;
        max = Math.max(max, Math.abs(lap + rho[idx]));
      }
    }
    return max;
  }
}

6. Dirichlet and Neumann Boundary Conditions

Ghost Cells for Neumann Condition

7. Equipotential Lines and Field Lines: JS Implementation

Equipotential Lines

Isoline at level c: the set of points where φ(x,y) = c. The Marching Squares algorithm iterates over grid squares and determines the level-crossing configuration:

function marchingSquaresIsolines(grid, N, levels) {
  const segments = [];
  for (const level of levels) {
    for (let i = 0; i < N - 1; i++) {
      for (let j = 0; j < N - 1; j++) {
        const v = [
          grid[i*N+j], grid[i*N+j+1],
          grid[(i+1)*N+j], grid[(i+1)*N+j+1]
        ];
        const code = (+(v[0]>level)) | (+(v[1]>level))<<1
                    | (+(v[2]>level))<<2 | (+(v[3]>level))<<3;
        if (code === 0 || code === 15) continue; // inside or outside
        const pts = msLookup(code, v, level, i, j); // lookup table
        if (pts) segments.push(...pts);
      }
    }
  }
  return segments;
}

Electric Field Lines

A field line is a curve tangent to the vector E = −∇φ at every point. Integrate using RK4 from a starting point (or a set of points around charges):

function traceFieldLine(phi, N, startI, startJ, stepSize = 0.5, maxSteps = 2000) {
  const path = [[startI, startJ]];
  let [ci, cj] = [startI, startJ];

  const grad = (i, j) => {
    const gx = (phi[Math.round(i+1)*N + Math.round(j)]
               - phi[Math.round(i-1)*N + Math.round(j)]) / 2;
    const gy = (phi[Math.round(i)*N + Math.round(j+1)]
               - phi[Math.round(i)*N + Math.round(j-1)]) / 2;
    return [gx, gy];
  };

  for (let s = 0; s < maxSteps; s++) {
    const [gx, gy] = grad(ci, cj);
    const mag = Math.hypot(gx, gy);
    if (mag < 1e-6) break;
    ci -= stepSize * gx / mag; // E = -∇φ → negate (line from +)
    cj -= stepSize * gy / mag;
    if (ci < 0 || ci >= N || cj < 0 || cj >= N) break;
    path.push([ci, cj]);
  }
  return path;
}

8. Multigrid and Applications

Multigrid

The reason SOR is slow: it only eliminates "short-wave" errors per iteration, while "long-wave" (smooth) errors persist. Multigrid fixes this: it moves to coarser grids where long-wave components become short-wave and are quickly suppressed (V-cycle or W-cycle):