Research and numerical validation

Generated and checked on 18 October 2026. The baseline charges are fixed at (+1, (1,0,0)), (-1, (0,1,0)), (+1, (-1,0,0)), (-1, (0,-1,0)).

Explore the results

• Interactive 3D viewer
• 2D field and potential figure
• Measured validation data (JSON)
• Independent baseline line audit (JSON)
• Charge-variation audits (JSON)
• Python source, tests, and instructions (ZIP)

Drag to orbit in 3D, scroll to zoom, and use the legend to toggle layers. Select a charge setup and use the local equipotentials button for potential surfaces. Source markers show point locations, not physically sized spheres.

Physical assumptions and units

This is isolated, static, vacuum 3D point-charge electrostatics, with potential zero at infinity. The 2D figure samples that same field on z=0. It does not use a two-dimensional logarithmic Green function or infinite line charges. Charges do not move: field lines show the field direction, not particle trajectories.

For target r and source positions r_i, we evaluate Coulomb superposition directly:

d_i = r - r_i
V(r) = k sum_i q_i / |d_i|
E(r) = k sum_i q_i d_i / |d_i|^3 = -gradient(V)
k = 1 / (4 pi epsilon_0)

These laws follow the textbook treatments in Feynman II, Chapter 4 and Fitzpatrick, Electric Field. The electrostatic differential equations are described in Feynman II, Section 6-1.

The calculation uses dimensionless coordinates r/a and charges q/q0 with k=q0=a=1. If the coordinate scale is a metres and the unit charge is q0 coulombs, multiply plotted potentials by k*q0/a volts and fields by k*q0/a^2 N/C. No physical charge magnitude or distance unit was assumed from the user's ±1.

The exact field is undefined at a source; the solver returns NaNs there. The default 3D tracer resolves approaches to charges down to radius 0.03a; the static 2D drawing excludes radius 0.07a disks. The default local 3D potential- surface mesh excludes radius 0.03a. No softening parameter changes the Coulomb kernel. Source markers are larger visual symbols; they do not represent charged bodies with finite physical radii. The static 2D image is a finite zoomed window of the same unbounded field.

Symmetric textbook sampling and unbounded curves

The baseline displays 24 ordinary charge-to-charge curves: twelve per positive charge, twelve per negative charge, and six for each positive/negative pairing. Seed directions are equal-solid-angle cell centers with x-component u = -2/3, 0, 2/3 and four azimuth centers separated by 90 degrees. Their transverse components stay away from zero. Exact null separatrices and the slow, unusually long excursions near them are deliberately absent from this presentation. They remain available in the separate research sampler.

Complete baseline paths preserve x, y, and z reflections exactly. Swapping x and y exchanges opposite charges, so E(Sr)=-S E(r). The tracer integrates three canonical half-curves to x=y, then constructs the remaining halves by that swap and direction reversal. At the midpoint E is proportional to (1,-1,0), so both halves join with the same tangent. Reflections generate the full family. An independent, stricter integration of each complete curve checks this symmetry construction; it does not replace Coulomb's law.

The planar figure similarly shows sixteen symmetry-completed ordinary paths. Its image window is fitted to these paths; it is not an integration boundary.

Selectable charge variations

The Setup dropdown switches between eight separate Python runs. Coordinates are unchanged; multiplying a negative charge means increasing its magnitude while retaining its negative sign. P1=(1,0,0), P2=(-1,0,0), N1=(0,1,0), and N2=(0,-1,0). The baseline 2D figure remains the equal-charge configuration.

Every run uses the actual Coulomb sum for its charge array, including the potential surfaces, arrows, colors, and path directions. Increasing a charge can break symmetries: a variation preserves only symmetries that its actual charge distribution still has. The equal-charge null-axis formula is not applied to unequal charges.

The variation sampler traces candidate symmetry orbits and selects whole orbits with an integer constraint: twelve displayed lines per unit magnitude at every positive and negative charge, with twelve per unit net charge continuing to/from infinity. Selection uses SciPy's documented mixed-integer linear optimizer. It chooses among already integrated physical curves; it never redirects a curve to satisfy a quota. Counts follow the textbook convention, while local line spacing is an illustrative sample, not a calibrated flux-density map.

Field lines can connect to infinity when charge is unbalanced, as explained in OpenStax's field-line rules. The general far-field interpretation follows the multipole expansion: nonzero net charge contributes a leading monopole, whereas the balanced, unequal adjacent pair has a nonzero leading dipole. The baseline's quadrupole asymptotics below do not apply to every variation.

Continuation and finite display

The baseline/research tracer has no fixed coordinate boundary and no global |E| stopping threshold. Adaptive integration continues through successively longer arclength intervals until a negative charge is reached or a known zero-field limiting point is resolved. A finite computational work budget is still necessary: exhausting it would be reported as unresolved, never as a physical endpoint. All 24 baseline curves end at negative-charge approach spheres; none has a boundary, null, or unresolved-work ending.

Unequal-charge runs can have paths to or from infinity. A finite viewer cannot contain an infinite curve: its computed tail stops only after a global radial escape certificate, and the marker/label explicitly denotes continuation. Near-charge and full-tail views are display choices, not new physical walls. Negative-net-charge runs integrate backward from negative charges and then reverse their arrays, so the displayed arrows still follow E.

Here is the certificate, derived directly from the Coulomb kernel. All source positions satisfy |ri| <= a, with a=1. Let Q=sum(qi), A=sum(|qi|), and R=|r|>a. The operator norm of the derivative of g(r)=r/|r|^3 is 2/|r|^3. Integrating this bound along each source displacement gives:

|E(r) - Q*r/R^3| <= A * [(R-a)^(-2) - R^(-2)]
B(R) = (A/|Q|) * [(R/(R-a))^2 - 1]

When B<1, the field's radial component has the sign of Q at every point outside that radius, and the bound decreases outward. Therefore the outward tracing direction cannot return to the source region. This certifies escape; it does not approximate the remaining curve by an exact straight ray or claim that a finite endpoint is literally infinity. The implementation and independent audit record the bound at each displayed continuation.

Optional research sampler: zero-field limits

The separate Fibonacci sampler retains the difficult cases for numerical research. With 32 seeds per positive charge, start radius 0.1a and stop radius 0.01a, it yields 62 charge captures and two null limits. Its largest excursion reaches 17.52a from the origin. These 64 curves are audited separately and are not part of the default textbook presentation.

Exactly on the z-axis E=0. The normalized arclength equation has no defined direction there, but its curve can have a well-defined limiting endpoint. For the two exactly symmetric seeds in y=0, that plane is an analytic invariant. Near the axis, the tracer changes coordinates from arclength to x and solves dz/dx = Ez/Ex to x=0, using its removable endpoint derivative of zero. It resolves the approach and includes its limiting point; it does not draw a connection across the null to another branch.

This follows the electrostatic gradient relation in Feynman II, Section 6-1 and the uniqueness of smooth vector-field integral curves, discussed in this university treatment of ODE uniqueness. Under the regular parameterization dr/dlambda=E, a null is an equilibrium approached asymptotically, rather than a point that an ordinary solution crosses. Lambda is a curve parameter, not the time of a particle trajectory.

The cancellation-resistant invariant-plane field is derived directly from the four Coulomb terms. Let A=1+x*x+z*z, t=2*x/A, and define Hminus=(1-t)^(-3/2), Hplus=(1+t)^(-3/2):

even = Hminus + Hplus - 2
odd  = Hplus - Hminus
Ex = A^(-3/2) * (x*even + odd)
Ey = 0
Ez = A^(-3/2) * z*even

The binomial even/odd parts are evaluated without subtracting nearly equal constants. At small x, with S=1+z*z, Ex=-6*x/S^(5/2)+O(x^3) and Ez=15*z*x*x/S^(7/2)+O(x^4). Thus dz/dx=-(5/2)*z*x/S+O(x^3) tends smoothly to zero. The calculation preserves the exact y=0 symmetry only for the exact requested charge geometry; an approximately symmetric input is not forced onto that invariant plane. Independent high-precision Coulomb evaluations and a halved coordinate-chart switch check this treatment.

Literature review and algorithm choices

The review covers the governing laws, visual interpretation, numerical tracing, floating-point errors, GPU execution, and scalable alternatives.

For four sources, a mesh PDE solve adds discretization error and boundary choices without helping compute the requested field. Direct analytic summation is the reference. For N sources and M targets it costs O(MN). The CPU implementation vectorizes targets in chunks; the GPU kernel loops over sources independently for each target. Neither allocates an M-by-N displacement tensor. Full host input/output arrays still require O(M+N) storage.

Both implementations use float64 and compensated (Kahan) accumulation. The quadrupole's net field is much smaller than its individual terms far away, so relative error deteriorates with cancellation. Compensated summation cannot remove all cancellation or subtraction error. Absolute errors and scale-aware vector errors accompany relative asymptotic errors. Work with normalized coordinates: this implementation does not guarantee scale-safe intermediates throughout the entire float64 exponent range.

The field values are evaluated analytically up to floating-point rounding. Potential contours and 3D isosurfaces interpolate a sampled grid, and curves have ODE integration error. Those displayed surfaces have not been subjected to a separate mesh-refinement study; they are illustrations of the validated field. Plotly's isosurface extraction requires finite grid values. Omitted source- core grid samples, if present at the selected resolution, use signed display sentinels; an explicit 26-neighbor halo check ensures those cores and their immediate neighbors are beyond the displayed potential levels, so the substituted values cannot create a displayed surface. The underlying field/potential evaluator retains its exact singularity handling.

Derived predictions and measured validation

The following configuration-specific formulas are our algebraic specialization of the sourced Coulomb and multipole laws, rather than values quoted from a paper:

Total charge = 0; dipole moment = (0,0,0)
Q_ij = sum_n q_n (3 r_ni r_nj - |r_n|^2 delta_ij) = diag(6,-6,0)
V_far = 3 (x^2-y^2) / r^5
On z=0: V_far = 3 cos(2 theta)/r^3
         E_r = 9 cos(2 theta)/r^4; E_theta = 6 sin(2 theta)/r^4
Near the origin: V = 3(x^2-y^2) + O(r^4); E = (-6x,6y,0) + O(r^3)

Thus the planar origin is a saddle, and every point on the z-axis is an exact field null. On the planes y=±x the potential is zero, although the field is generally nonzero. In a generic direction the far potential decays as r^-3 and field magnitude as r^-4. Along special directions a leading coefficient can vanish; relative potential errors on nodal planes are not meaningful.

The automated tests cover single positive/negative charges, finite dipoles, exact axis/diagonal values, reflections, source/zero-charge handling, the center Jacobian, the full null line, quadrupole asymptotics, potential gradients, full 3D divergence/curl, Gauss surface flux, chunking, and actual GPU parity.

Measured results in validation.json:

The Gauss checks follow OpenStax, Applying Gauss's Law. They integrate over actual 3D spheres, not 2D circles. In a planar slice, dEx/dx+dEy/dy generally is not zero: the omitted dEz/dz completes the three-dimensional divergence. Ez being zero on z=0 does not make its z derivative zero.

The baseline's largest excursion reaches 2.17a from the origin. This value results from the chosen noncritical seed directions, rather than a spatial cutoff. The retained research run demonstrates complete excursions beyond the former 3a display box.

The independent line audit in reports/field_line_audit.json uses a separate scalar Coulomb implementation, compensated scalar sums, and high-precision arithmetic near cancellations. It checks endpoint signs, every displayed segment's tangent, proximity to other positive charges, invariant-plane preservation, null limiting-point refinement, and rendered arrow directions. Individual plot traces remain separate, so two projected curves cannot be joined by a serialized connector. A 2D projection of the 3D scene can still visually overlap distinct curves; the single-curve selector and top views make their actual connections inspectable.

These are analytic/software validation results, not experimental measurements.

The full suite passes 43 tests, including actual CUDA execution. The separate reports/scenario_audit.json reconstructs all 396 displayed curves from the eight actual HTML files. It checks their charge values, endpoint quotas, arrows, scalar potential decrease, every segment's tangent, remaining source symmetries, and infinity certificates. Independent Gauss integrals and far-field checks use each run's actual charges; CPU/GPU comparisons include the changed magnitudes. No selected path is an unresolved or null case. Across the seven changed setups, a stricter representative integration agrees within 7.69e-10a; the largest displayed chord-versus-field angle is 0.0448 degrees. These are measured sampled checks, not formal bounds on every point of a continuous curve.

GPU performance and larger systems

The local RTX 5060 Ti executed the CUDA kernel successfully, using CuPy 14.1.1 and CUDA runtime 12.9. Timings below are medians of three repetitions after warm-up, including input/output transfers, with both backends in float64:

These compare this vectorized Python/NumPy source loop with this CUDA kernel, not an optimized compiled or multithreaded CPU implementation. CPU/GPU maximum pointwise field differences scaled by max(|E_cpu|,1) were below 2.70e-15 in these benchmarks; larger absolute differences occur near large fields.

Curve tracing uses SciPy on CPU; browser rendering uses WebGL. Generating the 24 baseline curves took about 0.35 seconds in this run. GPU sampling speed does not describe total viewer generation or tracing speed. For larger N, first benchmark shared-memory tiling and retaining targets/results on-device; then compare a mature free-space 3D FMM against this direct reference. Tree/FMM, particle-mesh methods, moving particles, conductors, and dielectric interfaces were reviewed but are not implemented in this starting model.

Run and extend

Fresh environment (Windows PowerShell, compatible Python installed):

python -m venv .venv
.\.venv\Scripts\python.exe -m pip install -e ".[test]"
# Optional NVIDIA/CUDA 12 backend, using official PyPI distributions:
.\.venv\Scripts\python.exe -m pip install -e ".[test,gpu]"
.\run.ps1 -Task all -Backend gpu
# CPU-only path:
.\run.ps1 -Task all -Backend cpu

run.ps1 keeps temporary CUDA compilation and Matplotlib files inside .cache and fails on a failed command. Available tasks are test, plot, viewer, validate, and audit. The viewer task generates all eight setups. Explicit GPU requests do not silently fall back; auto reports a fallback warning. The default output viewer embeds its Plotly runtime.

Generate just one alternative, or reproduce the retained research sample:

.\.venv\Scripts\python.exe viewer3d.py --scenario p1x4 --backend gpu
.\.venv\Scripts\python.exe audit_scenarios.py

from field_tracing import TraceConfig, trace_field_lines
curves, metadata = trace_field_lines(config=TraceConfig(
    sampling="fibonacci", seeds_per_positive=32,
    start_radius=0.1, stop_radius=0.01))

Reusable numerical API:

import numpy as np
from electrostatics import evaluate

positions = np.array([[1., 0., 0.], [0., 1., 0.],
                      [-1., 0., 0.], [0., -1., 0.]])
charges = np.array([1., -1., 1., -1.])
targets = np.array([[2., 0., 0.], [.3, .2, .8]])
potential, field = evaluate(targets, positions, charges, backend="gpu")

evaluate accepts arbitrary source positions, signed charges, and target arrays whose last dimension is three. The numerical core generalizes to additional charges. The default 2D illustration and its analytic tests describe this specific quadrupole; updating just its constants is not a complete workflow for an arbitrary new physical model. The 3D Python functions accept custom source arrays. The viewer CLI also offers the eight documented precomputed setups.

Before adding another physical model, record its governing reference, boundary conditions, units, exact cases, and convergence study as required by AGENTS.md.

Static publication

package_site.py assembles the public viewer, 2D figure, research notes, measured validation JSON, and an allowlisted Python-source ZIP into site/dist. It does not copy the virtual environment, caches, credentials, or local deployment logs. The explicit custom-domain configuration is in wrangler.jsonc; it creates a separate Worker named electricfield-simulation for electricfield.santerihukari.com.

.\.venv\Scripts\python.exe package_site.py
# With the official Cloudflare Wrangler CLI installed and authenticated:
npx wrangler deploy --dry-run
npx wrangler deploy

The numerical model is computed by Python before publication. The public static viewer displays those precomputed results locally using browser WebGL; it does not run CUDA or Python on the visitor's device.