Current Distribution via Magnetic Energy Minimization

How Does Current Actually “Decide” Where to Flow?

When current flows through a conductor, it doesn’t just wander randomly — it organizes itself to minimize the total magnetic energy stored in the system. This is a beautiful variational principle that sits at the heart of classical electromagnetism. In this post, we’ll build intuition for this principle, derive the governing equations, and solve a concrete example numerically in Python.


The Physics: Why Magnetic Energy Minimization?

In a resistive conductor carrying steady current, the current distribution $\mathbf{J}(\mathbf{r})$ arranges itself to minimize Ohmic dissipation — equivalently stated via the variational principle for magnetic energy. For a linear, isotropic conductor with resistivity $\rho$, the current distribution satisfies:

$$\nabla \times \mathbf{B} = \mu_0 \mathbf{J}$$

$$\nabla \cdot \mathbf{J} = 0 \quad \text{(continuity)}$$

$$\mathbf{J} = \sigma \mathbf{E} = -\sigma \nabla \phi$$

The total magnetic energy stored in the system is:

$$W = \frac{1}{2\mu_0} \int_{\text{all space}} |\mathbf{B}|^2 , dV$$

The Thomson’s theorem (Lord Kelvin, 1849) states that the actual current distribution minimizes the total Ohmic dissipation:

$$P = \int_V \frac{|\mathbf{J}|^2}{\sigma} , dV$$

subject to $\nabla \cdot \mathbf{J} = 0$ and the boundary conditions. This is equivalent to solving Laplace’s equation for the electric potential:

$$\nabla^2 \phi = 0 \quad \text{inside the conductor}$$


Concrete Example: Current Flow in a 2D Rectangular Conductor with a Hole

We consider a rectangular conducting plate with a circular hole punched through it. Current enters from the left edge and exits from the right edge. The hole forces current to redistribute — it must flow around the obstacle, and the distribution is precisely the one that minimizes magnetic energy.

Setup:

  • Domain: $[0, L_x] \times [0, L_y]$ with $L_x = 2.0$, $L_y = 1.0$
  • A circular hole of radius $r = 0.2$ centered at $(1.0, 0.5)$
  • Left boundary: $\phi = 1$ (current inlet)
  • Right boundary: $\phi = 0$ (current outlet)
  • Top/bottom boundaries: $\partial\phi/\partial n = 0$ (no normal current)
  • Inside hole: Neumann condition $\partial\phi/\partial n = 0$ (insulating void)

The governing equation is Laplace’s equation, and we solve it on a finite-difference grid with the hole masked out.


Python Code

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.colors import LogNorm
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
from mpl_toolkits.mplot3d import Axes3D

# ============================================================
# Parameters
# ============================================================
Nx, Ny = 200, 100 # grid resolution
Lx, Ly = 2.0, 1.0 # domain size [m]
dx = Lx / (Nx - 1)
dy = Ly / (Ny - 1)
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y) # shape (Ny, Nx)

# Hole parameters
xh, yh, rh = 1.0, 0.5, 0.2 # center and radius of hole

# Mask: True where the conductor IS present (not in hole)
dist2 = (X - xh)**2 + (Y - yh)**2
conductor = (dist2 > rh**2) # False inside hole

# ============================================================
# Build sparse linear system for Laplace: ∇²φ = 0
# ============================================================
N = Nx * Ny

def idx(i, j):
"""Flat index: i=col(x), j=row(y)"""
return j * Nx + i

A = lil_matrix((N, N))
b = np.zeros(N)

for j in range(Ny):
for i in range(Nx):
k = idx(i, j)

# ---- Left boundary: Dirichlet φ = 1 ----
if i == 0:
A[k, k] = 1.0
b[k] = 1.0
continue

# ---- Right boundary: Dirichlet φ = 0 ----
if i == Nx - 1:
A[k, k] = 1.0
b[k] = 0.0
continue

# ---- Inside hole: Dirichlet φ = ghost (isolated) ----
# We handle hole nodes by zero-current: set φ equal to average
# of neighbours that are also in the hole, or use Neumann via
# "mirroring" — simplest robust approach: just exclude with avg
if not conductor[j, i]:
# Assign φ by averaging available neighbours (ghost-cell)
neighbours = []
if i > 0: neighbours.append(idx(i-1, j))
if i < Nx-1: neighbours.append(idx(i+1, j))
if j > 0: neighbours.append(idx(i, j-1))
if j < Ny-1: neighbours.append(idx(i, j+1))
A[k, k] = len(neighbours)
for n in neighbours:
A[k, n] = -1.0
b[k] = 0.0
continue

# ---- Interior conductor nodes: 5-point Laplacian ----
# d²φ/dx² + d²φ/dy² = 0 (uniform σ)
# Coefficients (with possible Neumann at top/bottom)
coeffs = {}
rhs_k = 0.0

# x-direction (always interior in i since i=0,Nx-1 handled)
# left neighbour
if conductor[j, i-1]:
coeffs[idx(i-1, j)] = coeffs.get(idx(i-1, j), 0) + 1/dx**2
else:
# Neumann at hole boundary: mirror → contributes to diagonal
coeffs[k] = coeffs.get(k, 0) + 1/dx**2 # mirror adds to self
# right neighbour
if conductor[j, i+1]:
coeffs[idx(i+1, j)] = coeffs.get(idx(i+1, j), 0) + 1/dx**2
else:
coeffs[k] = coeffs.get(k, 0) + 1/dx**2

# y-direction
if j == 0:
# Neumann: ∂φ/∂y = 0 → mirror
coeffs[k] = coeffs.get(k, 0) + 1/dy**2
elif conductor[j-1, i]:
coeffs[idx(i, j-1)] = coeffs.get(idx(i, j-1), 0) + 1/dy**2
else:
coeffs[k] = coeffs.get(k, 0) + 1/dy**2

if j == Ny - 1:
coeffs[k] = coeffs.get(k, 0) + 1/dy**2
elif conductor[j+1, i]:
coeffs[idx(i, j+1)] = coeffs.get(idx(i, j+1), 0) + 1/dy**2
else:
coeffs[k] = coeffs.get(k, 0) + 1/dy**2

# Diagonal = sum of positive coefficients
diag = sum(coeffs.values())
A[k, k] = -diag
for col, val in coeffs.items():
if col != k:
A[k, col] = val

b[k] = rhs_k

# ============================================================
# Solve the linear system
# ============================================================
print("Converting to CSR and solving...")
A_csr = A.tocsr()
phi_flat = spsolve(A_csr, b)
phi = phi_flat.reshape(Ny, Nx)

# Mask hole region for plotting
phi_plot = np.where(conductor, phi, np.nan)

# ============================================================
# Compute current density J = -σ ∇φ (set σ=1 for normalised)
# ============================================================
# Use np.gradient (central differences, handles edges)
dphi_dy, dphi_dx = np.gradient(phi, dy, dx) # note: gradient(f, *varargs)
Jx = -dphi_dx
Jy = -dphi_dy

# Mask hole
Jx_plot = np.where(conductor, Jx, np.nan)
Jy_plot = np.where(conductor, Jy, np.nan)
J_mag = np.sqrt(Jx_plot**2 + Jy_plot**2)

# ============================================================
# Magnetic energy density (2D proxy): u_B ∝ |J|²
# (in 2D, B_z = μ₀ × stream function; energy ∝ |J|²)
# ============================================================
u_B = 0.5 * J_mag**2 # proportional to magnetic energy density

# Total "magnetic energy" (integrated, normalized)
u_B_vals = u_B[conductor]
total_W = np.nansum(u_B) * dx * dy
print(f"Normalised total magnetic energy W = {total_W:.6f} J/m")

# ============================================================
# === FIGURE 1: Potential, Current vectors, |J| heatmap ===
# ============================================================
fig, axes = plt.subplots(2, 2, figsize=(14, 9))
fig.suptitle("Current Distribution via Magnetic Energy Minimization\n"
"Rectangular Conductor with Circular Hole", fontsize=14, fontweight='bold')

hole_patch = lambda ax: ax.add_patch(
patches.Circle((xh, yh), rh, color='white', ec='black', lw=1.5, zorder=5))

# --- (a) Electric potential φ ---
ax = axes[0, 0]
cf = ax.contourf(X, Y, phi_plot, levels=40, cmap='RdYlBu_r')
ax.contour(X, Y, phi_plot, levels=15, colors='k', linewidths=0.4, alpha=0.5)
plt.colorbar(cf, ax=ax, label='φ (V)')
hole_patch(ax)
ax.set_title('(a) Electric Potential φ')
ax.set_xlabel('x (m)'); ax.set_ylabel('y (m)')
ax.set_aspect('equal')

# --- (b) Current density magnitude |J| ---
ax = axes[0, 1]
cf2 = ax.contourf(X, Y, J_mag, levels=40, cmap='hot_r')
plt.colorbar(cf2, ax=ax, label='|J| (A/m²)')
# Streamlines
ax.streamplot(x, y, Jx_plot, Jy_plot, density=1.4, color='cyan',
linewidth=0.8, arrowsize=0.8)
hole_patch(ax)
ax.set_title('(b) Current Density |J| + Streamlines')
ax.set_xlabel('x (m)'); ax.set_ylabel('y (m)')
ax.set_aspect('equal')

# --- (c) Current vectors (quiver, subsampled) ---
ax = axes[1, 0]
step = 8
cf3 = ax.contourf(X, Y, J_mag, levels=40, cmap='YlOrRd')
plt.colorbar(cf3, ax=ax, label='|J| (A/m²)')
ax.quiver(X[::step, ::step], Y[::step, ::step],
Jx_plot[::step, ::step], Jy_plot[::step, ::step],
scale=6, width=0.003, color='navy', alpha=0.8)
hole_patch(ax)
ax.set_title('(c) Current Vectors (quiver)')
ax.set_xlabel('x (m)'); ax.set_ylabel('y (m)')
ax.set_aspect('equal')

# --- (d) Magnetic energy density u_B ∝ |J|²/2 ---
ax = axes[1, 1]
cf4 = ax.contourf(X, Y, u_B, levels=40, cmap='plasma')
plt.colorbar(cf4, ax=ax, label='u_B ∝ |J|²/2')
hole_patch(ax)
ax.set_title('(d) Magnetic Energy Density u_B')
ax.set_xlabel('x (m)'); ax.set_ylabel('y (m)')
ax.set_aspect('equal')

plt.tight_layout()
plt.savefig('current_distribution_2d.png', dpi=150, bbox_inches='tight')
plt.show()
print("Figure 1 saved.")

# ============================================================
# === FIGURE 2: 3D Surface Plots ===
# ============================================================
fig2 = plt.figure(figsize=(16, 6))

# Subsample for 3D speed
s = 2
Xs, Ys = X[::s, ::s], Y[::s, ::s]
phi_s = phi_plot[::s, ::s]
uB_s = u_B[::s, ::s]

# --- 3D: Electric potential ---
ax3d1 = fig2.add_subplot(121, projection='3d')
surf1 = ax3d1.plot_surface(Xs, Ys, phi_s, cmap='RdYlBu_r',
edgecolor='none', alpha=0.9)
fig2.colorbar(surf1, ax=ax3d1, shrink=0.5, label='φ (V)')
ax3d1.set_title('3D Electric Potential φ', fontsize=12)
ax3d1.set_xlabel('x (m)'); ax3d1.set_ylabel('y (m)'); ax3d1.set_zlabel('φ (V)')
ax3d1.view_init(elev=30, azim=-60)

# --- 3D: Magnetic energy density ---
ax3d2 = fig2.add_subplot(122, projection='3d')
surf2 = ax3d2.plot_surface(Xs, Ys, uB_s, cmap='plasma',
edgecolor='none', alpha=0.9)
fig2.colorbar(surf2, ax=ax3d2, shrink=0.5, label='u_B')
ax3d2.set_title('3D Magnetic Energy Density u_B', fontsize=12)
ax3d2.set_xlabel('x (m)'); ax3d2.set_ylabel('y (m)'); ax3d2.set_zlabel('u_B')
ax3d2.view_init(elev=35, azim=-50)

plt.suptitle('3D Views — Magnetic Energy Minimization', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig('current_distribution_3d.png', dpi=150, bbox_inches='tight')
plt.show()
print("Figure 2 saved.")

# ============================================================
# === FIGURE 3: Cross-sections along y ===
# ============================================================
fig3, axes3 = plt.subplots(1, 2, figsize=(13, 5))
fig3.suptitle('Cross-sectional Analysis', fontsize=13, fontweight='bold')

y_slices = [0.25, 0.5, 0.75] # y/Ly = 0.25, 0.5, 0.75
colors = ['royalblue', 'crimson', 'forestgreen']

for yval, col in zip(y_slices, colors):
jy_idx = np.argmin(np.abs(y - yval))
Jx_row = Jx_plot[jy_idx, :]
uB_row = u_B[jy_idx, :]
lbl = f'y = {yval:.2f} m'
axes3[0].plot(x, Jx_row, color=col, lw=1.8, label=lbl)
axes3[1].plot(x, uB_row, color=col, lw=1.8, label=lbl)

axes3[0].axvline(xh - rh, color='gray', ls='--', alpha=0.6, label='Hole edges')
axes3[0].axvline(xh + rh, color='gray', ls='--', alpha=0.6)
axes3[0].set_title('Jx along horizontal slices')
axes3[0].set_xlabel('x (m)'); axes3[0].set_ylabel('Jx (A/m²)')
axes3[0].legend(); axes3[0].grid(alpha=0.3)

axes3[1].axvline(xh - rh, color='gray', ls='--', alpha=0.6, label='Hole edges')
axes3[1].axvline(xh + rh, color='gray', ls='--', alpha=0.6)
axes3[1].set_title('Magnetic energy density u_B along horizontal slices')
axes3[1].set_xlabel('x (m)'); axes3[1].set_ylabel('u_B')
axes3[1].legend(); axes3[1].grid(alpha=0.3)

plt.tight_layout()
plt.savefig('current_cross_sections.png', dpi=150, bbox_inches='tight')
plt.show()
print("Figure 3 saved.")

Code Walkthrough — Section by Section

Grid and mask setup

We discretize the $[0, 2] \times [0, 1]$ domain into a $200 \times 100$ uniform grid. A boolean array conductor flags every node that lies outside the circular hole using the Euclidean distance condition $(x - x_h)^2 + (y - y_h)^2 > r_h^2$. This mask is the backbone of every subsequent boundary-condition decision.

Sparse linear system construction

The 5-point finite-difference stencil for $\nabla^2 \phi = 0$ is:

$$\frac{\phi_{i+1,j} - 2\phi_{i,j} + \phi_{i-1,j}}{\Delta x^2} + \frac{\phi_{i,j+1} - 2\phi_{i,j} + \phi_{i,j-1}}{\Delta y^2} = 0$$

We build this as a sparse lil_matrix (List-of-Lists, efficient for random-access insertion). Three types of nodes are handled:

  • Dirichlet nodes (left/right boundaries): the row is trivially $\phi = $ const.
  • Hole nodes: we impose the Laplacian locally, which naturally enforces $\partial\phi/\partial n = 0$ through the “mirror” strategy — a neighbour that lies inside the hole is replaced by the current node’s own value, effectively reflecting the gradient and producing zero normal flux.
  • Neumann nodes (top/bottom edges): the same mirror reflection, so $\partial\phi/\partial y = 0$ at $y = 0$ and $y = L_y$.

Solving

The matrix is converted to CSR (Compressed Sparse Row) format and passed to scipy.sparse.linalg.spsolve, a direct sparse LU solver. For a $200 \times 100 = 20{,}000$ node system this completes in well under a second.

Current density and energy

$$\mathbf{J} = -\sigma \nabla \phi, \quad \sigma = 1 \text{ (normalised)}$$

np.gradient computes central differences on the interior and one-sided differences at the edges, giving $\partial\phi/\partial x$ and $\partial\phi/\partial y$ everywhere. The magnetic energy density proxy in 2D is:

$$u_B \propto \frac{1}{2}|\mathbf{J}|^2$$

The total normalised energy:

$$W = \int_\Omega u_B , dA \approx \sum_{i,j} u_B^{(i,j)} , \Delta x , \Delta y$$

is printed to confirm the solution’s scale.


What the Graphs Tell Us

Figure 1 — The four panels

Panel (a) — Electric potential: The equipotential lines smoothly warp around the hole. Far from the hole they are nearly vertical (as expected for uniform flow), but near the obstacle they bow outward — the field “sees” the insulating boundary and bends accordingly.

Panel (b) — |J| heatmap + streamlines: This is the most physically striking panel. The current density peaks sharply at the top and bottom of the hole, reaching values roughly 50–70% higher than the undisturbed uniform value. This is the electromagnetic analogue of stress concentration in mechanics. Streamlines confirm that no current enters the hole.

Panel (c) — Quiver plot: The vector arrows make the lateral deflection of the current patently obvious. Well upstream and downstream the arrows are horizontal and uniform. Approaching the hole, they fan outward; past it, they fan back inward.

Panel (d) — Magnetic energy density: The energy is concentrated at the poles of the hole (top and bottom edges). This is exactly where $|\mathbf{J}|^2$ is largest, and the minimization principle tells us: the actual distribution is the unique one that makes this concentration as small as physically possible given the constraints.

Figure 2 — 3D surfaces

The 3D potential surface is a smoothly warped “ramp” from $\phi = 1$ to $\phi = 0$, with the dip and rise near the hole visible as a gentle saddle. The 3D energy density surface shows two sharp peaks — the energy “mountains” at the hole’s poles — surrounded by a flat low-energy plateau. These mountains cannot be eliminated, but the variational principle ensures they are as low as possible.

Figure 3 — Cross-sections

The horizontal slice at $y = 0.5$ (through the hole’s equator) shows $J_x$ dropping to zero inside the hole and spiking at the edges $x = x_h \pm r_h$. The slices at $y = 0.25$ and $y = 0.75$ (passing above and below the hole) show elevated $J_x$ — the current that was “pushed out” of the hole has redistributed into these channels. The energy density cross-sections mirror this: a double peak where the current squeezes past the obstacle.


Key Takeaways

The minimum magnetic energy principle is not just a mathematical curiosity — it is the physical law that explains phenomena like:

$$\text{Skin effect: } \delta = \sqrt{\frac{2}{\omega \mu \sigma}}$$

where high-frequency currents hug the surface of a conductor to minimize inductance (and thus magnetic energy storage). It explains why current avoids sharp re-entrant corners, how lightning rods concentrate charge, and why PCB trace impedance depends on geometry.

The finite-difference approach here scales to 3D with the same logic (a 7-point stencil replaces the 5-point one), and the sparse solver handles millions of unknowns efficiently with iterative methods like scipy.sparse.linalg.minres or AMGCL.


Execution Results

Converting to CSR and solving...
Normalised total magnetic energy W = 2.694514 J/m

Figure 1 saved.

Figure 2 saved.

Figure 3 saved.