Solve a linear system with NumPy: the currents in a resistor network

The problem

You have a resistor network, a voltage source, and a ground, and want every node voltage and branch current. Kirchhoff's current law at each node \(i\) of unknown voltage, \(\sum_j (V_i - V_j)/R_{ij} = 0\), is one linear equation. Known voltages move to the right-hand side, and the unknown nodes, in the order of the list free, become the rows and columns of a linear system \(A\mathbf x = \mathbf b\), solved as in the eigenvalue tutorial.

The example is a Wheatstone bridge fed with 12 V through 10 Ω (source terminal 1, top 2, sides 3 and 4, ground 0). The bridge resistor carries current because 100/330 is not 220/150; at 726 Ω instead of 150 Ω it would carry none, and that is the code's test. Swap in your own network by editing the two lists.

The code

The two functions build and solve the system; the rest checks the result, prints it, and plots the currents.

Copy

import numpy as np
import matplotlib.pyplot as plt

# ---- network: replace these two lists with your own
resistors = [(1, 2, 10), (2, 3, 100), (2, 4, 220), (3, 0, 330), (4, 0, 150), (3, 4, 470)]  # (node, node, ohms)
fixed = {0: 0.0, 1: 12.0}                                    # known node voltages in V; node 0 is ground

# ---- Kirchhoff's current law at every node whose voltage is unknown
def node_equations(resistors, fixed):
    n = 1 + max(max(i, j) for i, j, _ in resistors)
    free = [k for k in range(n) if k not in fixed]           # unknown k is row and column free.index(k)
    A, b = np.zeros((len(free), len(free))), np.zeros(len(free))
    for i, j, R in resistors:
        for p, q in [(i, j), (j, i)]:                        # the current leaving p through R
            if p in fixed:
                continue
            A[free.index(p), free.index(p)] += 1 / R
            if q in fixed:
                b[free.index(p)] += fixed[q] / R             # a known voltage moves to the right-hand side
            else:
                A[free.index(p), free.index(q)] -= 1 / R
    return A, b, free

def node_voltages(resistors, fixed):
    A, b, free = node_equations(resistors, fixed)
    x = np.linalg.solve(A, b)
    V = np.zeros(len(free) + len(fixed))
    V[free] = x
    V[list(fixed)] = list(fixed.values())
    return V, np.abs(A @ x - b).max()

# ---- solve and check
V, residual = node_voltages(resistors, fixed)
print(f"residual of A x = b: {residual:.0e} A")
balanced = [(1, 2, 10), (2, 3, 100), (2, 4, 220), (3, 0, 330), (4, 0, 726), (3, 4, 470)]  # 100/330 = 220/726
V_test, _ = node_voltages(balanced, {0: 0.0, 1: 12.0})
R_arms = 1 / (1 / (100 + 330) + 1 / (220 + 726))           # no current in 470 Ω: two arms in parallel
V_hand = 12.0 * R_arms / (10 + R_arms) * 330 / (100 + 330)
print(f"balanced bridge: V3 = {V_test[3]:.3f} V, V4 = {V_test[4]:.3f} V, by hand {V_hand:.3f} V")

# ---- currents from Ohm's law, report and plot
for k, v in enumerate(V):
    print(f"V{k} = {v:7.3f} V")
labels, currents = [], []
for i, j, R in resistors:
    I = (V[i] - V[j]) / R
    p, q = (i, j) if I >= 0 else (j, i)                      # name the resistor in the direction of flow
    labels.append(f"{p}→{q}, {R:g} Ω")
    currents.append(1000 * abs(I))                           # A to mA
    print(f"{labels[-1]:>12s}: {currents[-1]:6.2f} mA")

fig, ax = plt.subplots(figsize=(7, 3.6), dpi=110)
y = np.arange(len(labels))
ax.barh(y, currents, color="#c8553d", height=0.6)
for yk, c in zip(y, currents):
    ax.text(c + 0.6, yk, f"{c:.2f} mA", va="center", color="#1f2a44")
ax.set(yticks=y, yticklabels=labels, xlabel="current / mA", xlim=(0, 1.18 * max(currents)))
ax.invert_yaxis()
ax.spines[["top", "right"]].set_visible(False)
plt.show()

residual of A x = b: 2e-17 A
balanced bridge: V3 = 8.908 V, V4 = 8.908 V, by hand 8.908 V
V0 =   0.000 V
V1 =  12.000 V
V2 =  11.403 V
V3 =   8.253 V
V4 =   5.202 V
   1→2, 10 Ω:  59.69 mA
  2→3, 100 Ω:  31.50 mA
  2→4, 220 Ω:  28.19 mA
  3→0, 330 Ω:  25.01 mA
  4→0, 150 Ω:  34.68 mA
  3→4, 470 Ω:   6.49 mA

The knobs

The list resistors takes any number of entries, with two resistors in parallel as two entries and the nodes numbered from 0 without gaps. Ground is only the zero of the scale: add 5 V to every entry of fixed, and every voltage rises by 5 V while no current changes. A second voltage source to ground is one more entry in fixed. Leave ground out while node 1 stays at 12 V, and nothing complains: every node sits at 12 V and every current is zero, because no path returns current to the source. The lists cannot hold two things. A current source adds its current to b in the row of the node it feeds and subtracts it in the row of the node it drains, and when that node is ground there is nothing to subtract, because ground has no row; the code has no list for it. A voltage source between two nodes that are both unknown needs its own current as an extra unknown, an extension called modified nodal analysis, which this recipe does not do.

The residual, below 10⁻¹⁵ A against 1.2 A, the largest entry of b, says that the computed voltages satisfy the equations to rounding. It says nothing about whether the equations are right. The balanced bridge does. With no current in the 470 Ω, both sides are plain voltage dividers, and V₃ and V₄ must come out at the 8.908 V computed by hand on the same line. Compare them with that value, not only with each other. Flip the sign of b and both read −8.908 V. Skip the (j, i) pass of the loop and both read 0 V. That known case tests the function, not your network. Neither check knows whether the list you typed is the circuit on your bench, and no computation can: read the six printed resistor lines against the schematic. Use np.linalg.solve, never inv(A) @ b, for the reason given in the eigenvalue tutorial: one factorization serves the one right-hand side you have. Dense storage stops paying off at a few thousand nodes. A row of \(A\) has nonzeros only for its node and that node's neighbors, yet a dense array stores every entry, and at 10,000 nodes that is 10⁸ numbers, 800 MB in double precision. Build \(A\) as a scipy.sparse matrix there and solve with scipy.sparse.linalg.spsolve.

Pitfalls

A node with nothing attached. Renumber node 4 as 5 throughout the list, the simplest way to leave a gap, and node 4 still exists but touches nothing. Its row of \(A\) is all zeros:

skipped = [(1, 2, 10), (2, 3, 100), (2, 5, 220), (3, 0, 330), (5, 0, 150), (3, 5, 470)]  # node 4 renumbered as 5
try:
    node_voltages(skipped, fixed)
except np.linalg.LinAlgError as err:
    print(f"LinAlgError: {err}")

LinAlgError: Singular matrix

The current law fixes only differences of voltages, so a node with no path to a fixed voltage has no level, and the matrix is singular. A whole part of the network cut off from every fixed node fails the same way. Number the nodes from 0 to the largest without gaps, and make sure every part of the network reaches a fixed node.

Mixed units. The report assumes ohms: 1000 * abs(I) turns amperes into the milliamperes it prints. A list written entirely in kilohms solves to the same voltages and prints every current exactly 1000 times too large, which you will notice. One value in kΩ among ohms is worse. Type the 470 Ω bridge resistor as 0.47 and it acts as a near short: V₃ and V₄ fall to 6.810 V and 6.798 V, the bridge current reads 24.67 mA instead of 6.49 mA, and the other five currents move by no recognizable factor. Write every value in ohms.

A wire typed as a zero resistance. Enter a wire as (3, 4, 0) and node_equations stops with ZeroDivisionError: division by zero. An ideal wire has infinite conductance, and that has no place in \(A\). A wire makes its two ends one node: give them one number and delete the entry.