Numerical flux and limits of symmetry

Gauss's law says ∮E⋅dA = Q(enc)/ε₀ for every closed surface, whatever the charges do. Nothing in it selects a shape. The step you learned for spheres, cylinders and sheets is separate and comes after: choose a surface on which |E| is constant and everywhere along the outward normal, so the integral becomes |E|A and the law can be rearranged into E = Q(enc)/(ε₀A). That rearrangement is what symmetry licenses, and it is the only thing symmetry licenses. Put a point charge at the centre of a cube and the law still fixes the total flux at q/ε₀, but no rearrangement follows, because the normal component on that surface is not constant: at a face corner it is 1/(3√3) = 0.192 of its value at the face centre, a spread of 5.2. Dividing flux by area gives q/(24ε₀a²) for a cube of side 2a, which is 48% below the face-centre field and 2.7 times the corner field — an average that matches the field at neither the face centre nor a corner — it is attained only on a ring between them, which is no place a formula would ever hand you.

The flux integral is a sum, so evaluate it as one. Partition the closed surface into N patches of area ΔAᵢ, take the centroid rᵢ of each, get the outward unit normal n̂ᵢ there, and add up E(rᵢ)⋅n̂ᵢ ΔAᵢ. That is the midpoint rule applied to a surface, and the two things that break it are a normal accidentally taken inward — which flips the sign of that patch — and a field evaluated somewhere other than the centroid. Try the crudest version on the cube: one patch per face, sampled at the face centre where E = q/(4πε₀a²) along the normal, times the face area 4a², gives q/(πε₀) per face. The exact value is q/(6ε₀), so the estimate is high by 6/π = 1.91, a 91% error from a sum that used the correct field at every point it sampled. The field it did not sample was weaker, and there was much more of it.

Refine and the error obeys a law. For an integrand the mesh resolves, the midpoint error scales as h² in the patch width, so on a fixed surface where N ∝ h⁻² it falls as 1/N. Take the same cube face: 2 × 2 patches give +3.96%, 4 × 4 give +1.53%, 8 × 8 give +0.383%, 16 × 16 give +0.0957%, and 32 × 32 give +0.0239%. From 8 × 8 onward each halving of h divides the error by almost exactly four, and p = log₂[(Φ(h) − Φ(h/2))/(Φ(h/2) − Φ(h/4))] returns 2.00. The two coarsest meshes do not obey it — the step from 2 × 2 to 4 × 4 gains only a factor 2.6 — because they are not yet resolving the integrand. So a convergence order is something you measure and confirm, not something you assume; quoting one from two coarse meshes is how a mesh error gets mistaken for physics.

Now break the symmetry outright. Around a cube of side 2.00 m centred on the origin, put +4.0 nC at (0.30, −0.20, 0.50) m, −1.5 nC at (−0.60, 0.40, −0.10) m, and +6.0 nC outside at (2.50, 0, 0) m. Enclosed charge is +2.5 nC, so the flux must be 2.5 nC/ε₀ = 282.358 N m² C⁻¹. The sum delivers it: 283.12 at 384 patches (+0.27%), 282.41 at 6144 (+0.017%), 282.361 at 98 304 (+0.0011%). The external charge is the interesting part. Its own contribution converges to zero — −0.34, then −0.021, then −0.0013 N m² C⁻¹ — but patch by patch it is large, entering through the near face and leaving through the far one. Its field is on the surface at full strength the whole time. Only the total cancels, which is exactly what the law claims and no more.

Uniform refinement is not always the fix, because the integrand has its own width. Put 1.0 nC in the same 2.00 m cube but only 0.05 m below one face. The normal component there is concentrated in a spot of about that size, so a mesh with h = 0.25 m steps straight over it and returns 84.0 against 112.94 N m² C⁻¹ — 25.6% low. Halving to h = 0.125 m still gives −10.4%. Only once h = 0.0625 m, comparable with the 0.05 m gap, does the error drop to −1.14%, and h = 0.031 m reaches −0.006%. The working rule is that h must be well below the distance from the surface to the nearest charge, which is an argument for adaptive refinement: subdivide the patches where the integrand is peaked and leave the quiet faces coarse. A charge sitting exactly on the surface has no finite answer at all; move the surface.

Treat the result as a test, not a proof. The flux is a single scalar collapsed out of the whole surface, so agreement to 10⁻⁵ does not mean the field is right patch by patch — opposite errors cancel inside the sum and leave no trace. What it does catch is cheap and worth catching: a normal taken inward, a charge left out of Q(enc), a charge wrongly counted as enclosed, a stale ε₀. And notice what the whole procedure never does. It never produces E. You must already have the field, from superposition or from a Poisson solver, before the first patch can be evaluated. Gauss's law is a constraint the solution has to satisfy, not an algorithm for finding it. Symmetry is the special case where the constraint is tight enough to be inverted; the rest of the time the law is exactly as true and simply will not do the work for you.