Skip to content

Instantly share code, notes, and snippets.

@mbillingr

mbillingr/fieldmodel.py

Last active May 24, 2022 10:51
Show Gist options
  • Save mbillingr/2bfc1b399daf146f22c05128f21ca493 to your computer and use it in GitHub Desktop.
Save mbillingr/2bfc1b399daf146f22c05128f21ca493 to your computer and use it in GitHub Desktop.
Download ZIP
Model the electrical potential on a finite grid of resistors... and use it to naively solve xkcd356
Raw
fieldmodel.py
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.image import imread
try:
import seaborn
except ImportError:
pass
class FieldModel:
""" Model of a regular resistor grid. All nodes in the grid are connected to their direct neighbors with resitors.
"""
def __init__(self, shape=(128, 128), extent=(1, 1)):
self._shape = np.asarray(shape)
self._extent = np.asarray(extent)
self.reset()
def reset(self):
self._step = self._extent / self._shape
self._potential_field = np.zeros(self._shape)
self._temporary_field = np.zeros(self._shape)
self._sx = np.ones(self._shape - [0, 1])
self._sy = np.ones(self._shape - [1, 0])
self._clamps = []
self._inductions = []
@property
def shape(self):
return self._shape
@shape.setter
def shape(self, s):
self._shape = np.asarray(s)
self.reset()
@property
def extent(self):
return self._extent
@extent.setter
def extent(self, e):
self._extent = np.asarray(e)
self.reset()
@property
def dx(self):
return self._step[0]
@property
def dy(self):
return self._step[1]
@property
def potential(self):
return self._potential_field
@property
def current(self):
return np.diff(self._potential_field, axis=1) * self._sx, np.diff(self._potential_field, axis=0) * self._sy
def get_netcurrent(self, x, y):
u0 = self._potential_field[x, y]
u1 = self._potential_field[x-1, y] if x > 0 else 0
u2 = self._potential_field[x+1, y] if x < self.shape[0]-1 else 0
u3 = self._potential_field[x, y-1] if y > 0 else 0
u4 = self._potential_field[x, y+1] if y < self.shape[1]-1 else 0
s1 = self._sy[x-1, y] if x > 0 else 0
s2 = self._sy[x, y] if x < self.shape[0]-1 else 0
s3 = self._sx[x, y-1] if y > 0 else 0
s4 = self._sx[x, y] if y < self.shape[1]-1 else 0
return (u1 - u0) * s1 + (u2 - u0) * s2 + (u3 - u0) * s3 + (u4 - u0) * s4
def add_clamp(self, x, y, u):
self._clamps.append((x, y, u))
def add_induction(self, x, y, d, i):
if d == 'r':
u = i / self._sx[x, y]
self._inductions.append((x, y, x, y + 1, u))
elif d == 'd':
u = i / self._sy[x, y]
self._inductions.append((x, y, x + 1, y, u))
elif d == 'l':
u = i / self._sx[x, y - 1]
self._inductions.append((x, y, x, y - 1, u))
elif d == 'u':
u = i / self._sy[x - 1, y]
self._inductions.append((x, y, x - 1, y, u))
else:
raise AssertionError("direction must be one of l/r/u/d.")
def set_conductivity(self, i, j, d, s):
if d == 'r':
self._sx[i, j] = s
elif d == 'd':
self._sy[i, j] = s
elif d == 'l':
self._sx[i, j - 1] = s
elif d == 'u':
self._sy[i - 1, j] = s
else:
raise AssertionError("direction must be one of l/r/u/d.")
def iterate(self, n_iter):
ws = self._sy[1:, 1:-1] + self._sy[:-1, 1:-1] + self._sx[1:-1, 1:] + self._sx[1:-1, :-1]
for _ in range(n_iter):
self._temporary_field[1:-1, 1:-1] = (self._potential_field[2:, 1:-1] * self._sy[1:, 1:-1] +
self._potential_field[:-2, 1:-1] * self._sy[:-1, 1:-1] +
self._potential_field[1:-1, 2:] * self._sx[1:-1, 1:] +
self._potential_field[1:-1, :-2] * self._sx[1:-1, :-1]) / ws
self.apply_boundary_closed(self._potential_field, self._temporary_field, self._sx, self._sy)
na = np.isnan(self._temporary_field)
self._temporary_field[na] = self._potential_field[na]
self._potential_field, self._temporary_field = self._temporary_field, self._potential_field
for i, j, u in self._clamps:
self._potential_field[i, j] = u
for i1, j1, i2, j2, du in self._inductions:
u = self._potential_field[i2, j2] - self._potential_field[i1, j1]
r = (du - u) / 2
self._potential_field[i1, j1] -= r
self._potential_field[i2, j2] += r
@staticmethod
def apply_boundary_zero(U, V, sx, sy):
V[0, 1:-1] =0
V[-1, 1:-1] = 0
V[1:-1, 0] = 0
V[1:-1, -1] = 0
V[0, 0] = 0
V[-1, 0] = 0
V[0, -1] = 0
V[-1, -1] = 0
return V
@staticmethod
def apply_boundary_open(U, V, sx, sy):
V[0, 1:-1] = U[1, 1:-1]
V[-1, 1:-1] = U[-2, 1:-1]
V[1:-1, 0] = U[1:-1, 1]
V[1:-1, -1] = U[1:-1, -2]
V[0, 0] = (U[0, 1] + U[1, 0]) * 0.5
V[-1, 0] = (U[-1, 1] + U[-2, 0]) * 0.5
V[0, -1] = (U[0, -2] + U[1, -1]) * 0.5
V[-1, -1] = (U[-1, -2] + U[-2, -1]) * 0.5
return V
@staticmethod
def apply_boundary_closed(U, V, sx, sy):
V[0, 1:-1] = (U[1, 1:-1] * sy[0, 1:-1] + U[0, 2:] * sx[0, 1:] + U[0, :-2] * sx[0, :-1]) / (sy[0, 1:-1] + sx[0, 1:] + sx[0, :-1])
V[-1, 1:-1] = (U[-2, 1:-1] * sy[-1, 1:-1] + U[-1, 2:] * sx[-1, 1:] + U[-1, :-2] * sx[-1, :-1]) / (sy[-1, 1:-1] + sx[-1, 1:] + sx[-1, :-1])
V[1:-1, 0] = (U[2:, 0] * sy[1:, 0] + U[:-2, 0] * sy[:-1, 0] + U[1:-1, 1] * sx[1:-1, 0]) / (sy[1:, 0] + sy[:-1, 0] + sx[1:-1, 0])
V[1:-1, -1] = (U[2:, -1] * sy[1:, -1] + U[:-2, -1] * sy[:-1, -1] + U[1:-1, -2] * sx[1:-1, -1]) / (sy[1:, -1] + sy[:-1, -1] + sx[1:-1, -1])
V[0, 0] = (U[1, 0] * sy[0, 0] + U[0, 1] * sx[0, 0]) / (sy[0, 0] + sx[0, 0])
V[-1, 0] = (U[-2, 0] * sy[-1, 0] + U[-1, 1] * sx[-1, 0]) / (sy[-1, 0] + sx[-1, 0])
V[0, -1] = (U[1, -1] * sy[0, -1] + U[0, -2] * sx[0, -1]) / (sy[0, -1] + sx[0, -1])
V[-1, -1] = (U[-2, -1] * sy[-1, -1] + U[-1, -2] * sx[-1, -1]) / (sy[-1, -1] + sx[-1, -1])
return V
def load_model(self, filename):
img = imread(filename)
r, g, b, a = np.rollaxis(img, axis=-1)
self.shape = img.shape[:2]
self.extent = self.shape * 1e-4 # 1 pixel is 0.1mm
# Conductivity values taken from
# http://www.itis.ethz.ch/virtual-population/tissue-properties/database/low-frequency-conductivity/
cond = np.zeros(r.shape)
cond[(r == 1) & (g == 0) & (b == 0)] = 3.69e-1 # White Matter
cond[(r == 1) & (g == 0) & (b == 1)] = 1.85e-1 # Grey Matter
cond[(r == 0) & (g == 0) & (b == 1)] = 1.79e0 # Cerebrospinal Fluid
cond[(r == 0) & (g == 1) & (b == 0)] = 9.5e-2 # Bone
cond[(r == 1) & (g == 1) & (b == 0)] = 1.2e-3 # Skin (wet)
self._sx = (cond[:, :-1] + cond[:, 1:]) * 0.5 / self.dx
self._sy = (cond[:-1, :] + cond[1:, :]) * 0.5 / self.dy
def __str__(self):
s = ''
for i in range(self.shape[0]):
s += '....' + ' '.join('{:+2.3f}V'.format(k) for k in self.potential[i])
s += '\n O'
for k, cond in zip(self.current[0][i], self._sx[i]):
if cond > 0:
s += '--------{}---------O'.format({-1: '<', 0: '=', 1: '>'}[np.sign(k)])
else:
s += ' O'
#s + ('--------{}---------'.join(['O'] * self.shape[1])).format(*[{-1: '<', 0: '=', 1: '>'}[np.sign(k)] for k in self.current[0][i]])
if i < self.shape[0] - 1:
s += '\n '
s += ' ' if self._sy[i][0] == 0 else '|'
for k, sx, sy in zip(abs(self.current[0][i]), self._sx[i], self._sy[i][1:]):
v_line = ' ' if sy == 0 else '|'
if sx == 0:
s += ' ' + v_line
else:
s += ' {:2.3f}A '.format(k) + v_line
s += '\n ' + (' '.join([' ' if sy == 0 else '|' for sy in self._sy[i]]))
s += '\n ' + (' '.join(['{} {:2.3f}A'] * self.shape[1])).format(*[b for a in zip([{-1: '^', 0: '=', 1: 'v'}[np.sign(k)] for k in self.current[1][i]],
np.abs(self.current[1][i])) for b in a])
s += '\n ' + (' '.join([' ' if sy == 0 else '|' for sy in self._sy[i]]))
else:
s += '\n '
for k, cond in zip(abs(self.current[0][i]), self._sx[i]):
if cond == 0:
s += ' '.format(k) + ' '
else:
s += ' {:2.3f}A '.format(k) + ' '
#s += '\n ' + ' '.join(' {:2.3f}A '.format(k) for k in abs(self.current[0][i])) + ' '
s += '\n'
return s
if __name__ == '__main__':
model = FieldModel(shape=(64, 64), extent=(4, 4))
model.add_clamp(32, 20, 1)
model.add_clamp(32, 44, 1)
l, r = 21, 43
t, b = 31, 33
for i in range(l, r):
for j in range(t, b):
model.set_conductivity(j, i, 'r', 100000)
model.set_conductivity(j, i, 'd', 100000)
model.set_conductivity(b, i, 'r', 100000)
model.iterate(10000)
U = model._potential_field
#levels = [-1, -0.5, -1e-1, -5e-2, -1e-2, -0.75e-2, -0.5e-2, -0.25e-2, -1e-3,
# 1e-3, 0.25e-2, 0.5e-2, 0.75e-2, 1e-2, 5e-2, 1e-1, 0.5, 1]
ls = np.logspace(-3, 0, 15)
levels = np.hstack([-ls[::-1], ls])
colors = seaborn.color_palette("coolwarm", len(levels)-1)
plt.contourf(U, vmin=-1, vmax=1, levels=levels, colors=colors, linewidth=0, zorder=1)
plt.plot([l, l, r, r, l], [t, b, b, t, t], 'k')
plt.show()
Raw
xkcd356.py
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from fieldmodel import FieldModel
# https://xkcd.com/356/
# Simulate a finite grid of resistors for increasing grid sizes and see where
# the value of the equivalent resistance converges...
# Solution: somewhere around 0.77, probably
r0 = 1
rs = []
xr = [0, 1, 2, 5, 10, 25, 50, 100, 250, 500, 1000]
for k in xr:
n, m = 2+k*2, 3+k*2
cy = (n-1) // 2
cx = (m-2) // 2
clamp_points = [(cy, cx), (cy+1, cx+2)]
model = FieldModel(shape=(n, m))
model._sx[:] = 1/r0
model._sy[:] = 1/r0
model.add_clamp(*clamp_points[0], 1)
model.add_clamp(*clamp_points[1], 0)
r_last = 1
r_total = 1e-9
while np.abs(r_total - r_last) >= 1e-9:
r_last = r_total
model.iterate(2)
r_total = (model.potential[clamp_points[1]] - model.potential[clamp_points[0]]) / model.get_netcurrent(cy, cx)
rs.append(r_total)
print(k, r_total)
plt.gca().axhline(rs[-1], ls='--', c='grey')
plt.plot(xr, rs)
plt.text(xr[-1]//2, rs[-1]-0.01, rs[-1], va='top', ha='right')
plt.ylim(0.7, 1.4)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment