Resistor mesh
You are encouraged to solve this task according to the task description, using any language you may know.
Given 10 × 10 grid nodes interconnected by 1Ω resistors as shown, find the resistance between point A and B.
See also [[1]]
Ada
<lang Ada>with Ada.Text_IO; use Ada.Text_IO; procedure ResistMesh is
H, W : constant Positive := 10; rowA, colA : constant Positive := 2; -- row/col indexed from 1 rowB : constant Positive := 7; colB : constant Positive := 8;
type Ntype is (A, B, Free); type Vtype is digits 15; type Node is record volt : Vtype := 0.0; name : Ntype := Free; end record; type NodeMesh is array (Positive range <>, Positive range <>) of Node; package IIO is new Ada.Text_IO.Float_IO (Vtype); mesh, dmesh : NodeMesh (1 .. H, 1 .. W); curA, curB, diff : Vtype;
procedure set_AB (mesh : in out NodeMesh) is begin mesh (rowA, colA).volt := 1.0; mesh (rowA, colA).name := A; mesh (rowB, colB).volt := -1.0; mesh (rowB, colB).name := B; end set_AB;
function sides (i : Positive; j : Positive) return Vtype is s : Integer := 0; begin if i /= 1 and i /= H then s := s + 2; else s := s + 1; end if; if j /= 1 and j /= W then s := s + 2; else s := s + 1; end if; return Vtype (s); end sides;
procedure calc_diff (mesh : NodeMesh; dmesh : out NodeMesh; total : out Vtype) is n : Natural; v : Vtype := 0.0; begin total := 0.0; for i in Positive range 1 .. H loop for j in Positive range 1 .. W loop n := 0; v := 0.0; if i /= 1 then v := v + mesh (i - 1, j).volt; n := n + 1; end if; if j /= 1 then v := v + mesh (i, j - 1).volt; n := n + 1; end if; if i < H then v := v + mesh (i + 1, j).volt; n := n + 1; end if; if j < W then v := v + mesh (i, j + 1).volt; n := n + 1; end if; v := mesh (i, j).volt - v / Vtype (n); dmesh (i, j).volt := v; if mesh (i, j).name = Free then total := total + v ** 2; end if; end loop; end loop; end calc_diff;
begin
loop set_AB (mesh); calc_diff (mesh, dmesh, diff); exit when diff < 1.0e-40; for i in Positive range 1 .. H loop for j in Positive range 1 .. W loop mesh (i, j).volt := mesh (i, j).volt - dmesh (i, j).volt; end loop; end loop; end loop;
curA := dmesh (rowA, colA).volt * sides (rowA, colA); curB := dmesh (rowB, colB).volt * sides (rowB, colB); diff := 4.0 / (curA - curB); IIO.Put (diff, Exp => 0); New_Line;
end ResistMesh;</lang>
- Output:
1.60899124173073
C
<lang c>#include <stdio.h>
- include <stdlib.h>
- define S 10
typedef struct { double v; int fixed; } node;
- define each(i, x) for(i = 0; i < x; i++)
node **alloc2(int w, int h) { int i; node **a = calloc(1, sizeof(node*)*h + sizeof(node)*w*h); each(i, h) a[i] = i ? a[i-1] + w : (node*)(a + h); return a; }
void set_boundary(node **m) { m[1][1].fixed = 1; m[1][1].v = 1; m[6][7].fixed = -1; m[6][7].v = -1; }
double calc_diff(node **m, node **d, int w, int h) { int i, j, n; double v, total = 0; each(i, h) each(j, w) { v = 0; n = 0; if (i) v += m[i-1][j].v, n++; if (j) v += m[i][j-1].v, n++; if (i+1 < h) v += m[i+1][j].v, n++; if (j+1 < w) v += m[i][j+1].v, n++;
d[i][j].v = v = m[i][j].v - v / n; if (!m[i][j].fixed) total += v * v; } return total; }
double iter(node **m, int w, int h) { node **d = alloc2(w, h); int i, j; double diff = 1e10; double cur[] = {0, 0, 0};
while (diff > 1e-24) { set_boundary(m); diff = calc_diff(m, d, w, h); each(i,h) each(j, w) m[i][j].v -= d[i][j].v; }
each(i, h) each(j, w) cur[ m[i][j].fixed + 1 ] += d[i][j].v * (!!i + !!j + (i < h-1) + (j < w -1));
free(d); return (cur[2] - cur[0])/2; }
int main() { node **mesh = alloc2(S, S); printf("R = %g\n", 2 / iter(mesh, S, S)); return 0; }</lang>
D
<lang d>import std.stdio, std.traits;
enum Node.FP differenceThreshold = 1e-40;
struct Node {
alias real FP; enum Kind : size_t { free, A, B }
FP voltage = 0.0; private Kind kind = Kind.free;
@property Kind fixed() const pure nothrow { return kind; }
}
Node.FP iter(size_t w, size_t h)(ref Node[w][h] m) pure nothrow {
static void enforceBoundaryConditions(ref Node[w][h] m) pure nothrow { m[1][1].voltage = 1.0; m[6][7].voltage = -1.0; }
static Node.FP calcDifference(in ref Node[w][h] m, ref Node[w][h] d) pure nothrow { Node.FP total = 0.0;
foreach (i; 0 .. h) foreach (j; 0 .. w) { Node.FP v = 0.0; { size_t n = 0; if (i != 0) { v += m[i - 1][j].voltage; n++; } if (j != 0) { v += m[i][j - 1].voltage; n++; } if (i < h-1) { v += m[i + 1][j].voltage; n++; } if (j < w-1) { v += m[i][j + 1].voltage; n++; } v = m[i][j].voltage - v / n; }
d[i][j].voltage = v; if (m[i][j].fixed == Node.Kind.free) total += v ^^ 2; }
return total; }
Node[w][h] difference;
while (true) { enforceBoundaryConditions(m); if (calcDifference(m, difference) < differenceThreshold) break; foreach (i, const(Node[]) di; difference) foreach (j, ref const(Node) dij; di) m[i][j].voltage -= dij.voltage; }
Node.FP[EnumMembers!(Node.Kind).length] cur = 0.0; foreach (i, const(Node[]) di; difference) foreach (j, ref const(Node) dij; di) cur[m[i][j].fixed] += dij.voltage * (!!i + !!j + (i < h-1) + (j < w-1));
return (cur[Node.Kind.A] - cur[Node.Kind.B]) / 2.0;
}
void main() {
enum size_t w = 10, h = w; Node[w][h] mesh;
// Set A and B Nodes. mesh[1][1] = Node( 1.0, Node.Kind.A); mesh[6][7] = Node(-1.0, Node.Kind.B);
writefln("R = %.19f", 2 / iter(mesh));
}</lang>
- Output:
R = 1.6089912417307296556
J
We represent the mesh as a Ybus matrix with B as the reference node and A as the first node and invert it to find the Z bus (which represents resistance). The first element of the first row of this Z matrix is the resistance between nodes A and B.
<lang J>nodes=: 10 10 #: i. 100 nodeA=: 1 1 nodeB=: 6 7
conn =: [: (#~ e.~/@|:~&0 2) ([ ,: +)"1 ref =: ~. nodeA,nodes-.nodeB wiring=: /:~ ref i. ,/ nodes conn"2 1 (,-)=i.2 Yii=: (* =@i.@#) _1 }. #/.~ {."1 wiring Yij=: - _1 _1 }. 1:`(<"1@[)`]}&(+/~ 0*i.1+#ref) wiring Y=: Yii+Yij</lang>
Here, the result of offset conn nodes
represents all pairs of nodes where we can connect the argument nodes to neighboring nodes at the specified offset, and wiring
is a list of index pairs representing all connections made by all resistors (note that each connection is represented twice -- node e is connected to node f AND node f is connected to node e). Yii contains the values for the diagonal elements of the Y bus while Yij contains the values for the off diagonal elements of the Y bus.
So:
<lang J> {.{. %. Y 1.60899</lang>
Or, if we want an exact answer (this is slow), we can assume our resistors are perfect:
<lang J> {.{.%. x:Y 455859137025721r283319837425200</lang>
(here, the letter 'r' separates the numerator from the denominator)
See http://jsoftware.com/pipermail/programming/2013-January/030892.html for a solution to the xkcd problem referenced in the task description.
Maxima
<lang maxima>/* Place a current souce between A and B, providing 1 A. Then we are really looking
for the potential at A and B, since I = R (V(B) - V(A)) where I is given and we want R. Atually, we will compute potential at each node, except A where we assume it's 0. Without with assumption, there would be infinitely many solutions since potential is known up to a constant. For A we will simply write the equation V(A) = 0, to keep the program simple. Hence, for a general grid of p rows and q columns, there are n = p * q nodes, so n unknowns, and n equations. Write Kirchhoff's current law at each node. Be careful with the node A (equation A = 0) and the node B (there is a constant current to add, from the source, that will go in the constant terms of the system). Finally, we have a n x n linear system of equations to solve. Simply use Maxima's builtin LU decomposition. Since all computations are exact, the result will be also exact, written as a fraction. Also, the program can work with any grid, and any two nodes on the grid.
For those who want more speed and less space, notice the system is sparse and necessarily symmetric, so one can use conjugate gradient or any other sparse symmetric solver. */
/* Auxiliary function to get rid of the borders */ ongrid(i, j, p, q) := is(i >= 1 and i <= p and j >= 1 and j <= q)$
grid_resistor(p, q, ai, aj, bi, bj) := block(
[n: p * q, A, B, M, k, c, V], A: zeromatrix(n, n), for i thru p do for j thru q do ( k: (i - 1) * q + j, if i = ai and j = aj then A[k, k]: 1 else ( c: 0, if ongrid(i + 1, j, p, q) then (c: c + 1, A[k, k + q]: -1), if ongrid(i - 1, j, p, q) then (c: c + 1, A[k, k - q]: -1), if ongrid(i, j + 1, p, q) then (c: c + 1, A[k, k + 1]: -1), if ongrid(i, j - 1, p, q) then (c: c + 1, A[k, k - 1]: -1), A[k, k]: c ) ), B: zeromatrix(n, 1), B[k: (bi - 1) * q + bj, 1]: 1, M: lu_factor(A), V: lu_backsub(M, B), V[k, 1]
)$
grid_resistor(10, 10, 2, 2, 8, 7); 455859137025721 / 283319837425200
bfloat(%), fpprec = 40; 1.608991241730729655954495520510088761201b0
/* Some larger example */ grid_resistor(20, 20, 1, 1, 20, 20); 129548954101732562831760781545158173626645023 / 33283688571680493510612137844679320717594861
bfloat(%), fpprec = 40; 3.89226554090400912102670691601064387507b0</lang>
Mathematica
<lang mathematica>gridresistor[p_, q_, ai_, aj_, bi_, bj_] :=
Block[{A, B, k, c, V}, A = ConstantArray[0, {p*q, p*q}]; Do[k = (i - 1) q + j; If[{i, j} == {ai, aj}, Ak, k = 1, c = 0; If[1 <= i + 1 <= p && 1 <= j <= q, c++; Ak, k + q = -1]; If[1 <= i - 1 <= p && 1 <= j <= q, c++; Ak, k - q = -1]; If[1 <= i <= p && 1 <= j + 1 <= q, c++; Ak, k + 1 = -1]; If[1 <= i <= p && 1 <= j - 1 <= q, c++; Ak, k - 1 = -1]; Ak, k = c], {i, p}, {j, q}]; B = SparseArray[(k = (bi - 1) q + bj) -> 1, p*q]; LinearSolve[A, B]k];
N[gridresistor[10, 10, 2, 2, 8, 7], 40]</lang>
- Output:
1.608991241730729655954495520510088761201
<lang mathematica>graphresistor[g_, a_, b_] :=
LinearSolve[ SparseArray[{{a, a} -> 1, {i_, i_} :> Length@AdjacencyList[g, i], Alternatives @@ Join[#, Reverse /@ #] &[ List @@@ EdgeList[VertexDelete[g, a]]] -> -1}, {VertexCount[ g], VertexCount[g]}], SparseArray[b -> 1, VertexCount[g]]]b;
N[graphresistor[GridGraph[{10, 10}], 12, 77], 40]</lang>
- Output:
1.608991241730729655954495520510088761201
Perl
<lang perl>use strict;
my ($w, $h) = (9, 9); my @v = map([ (0) x ($w + 1) ], 0 .. $h); # voltage my @f = map([ (0) x ($w + 1) ], 0 .. $h); # fixed condition my @d = map([ (0) x ($w + 1) ], 0 .. $h); # diff
my @n; # neighbors for my $i (0 .. $h) { push @{$n[$i][$_]}, [$i, $_ - 1] for 1 .. $w; push @{$n[$i][$_]}, [$i, $_ + 1] for 0 .. $w - 1; } for my $j (0 .. $w) { push @{$n[$_][$j]}, [$_ - 1, $j] for 1 .. $h; push @{$n[$_][$j]}, [$_ + 1, $j] for 0 .. $h - 1; }
sub set_boundary { $f[1][1] = 1; $f[6][7] = -1; $v[1][1] = 1; $v[6][7] = -1; }
sub calc_diff { my $total_diff; for my $i (0 .. $h) { for my $j (0 .. $w) { my ($p, $v) = $n[$i][$j]; $v += $v[$_->[0]][$_->[1]] for @$p; $d[$i][$j] = $v = $v[$i][$j] - $v / scalar(@$p); $total_diff += $v * $v unless $f[$i][$j]; } } $total_diff; }
sub iter { my $diff = 1; while ($diff > 1e-24) { # 1e-24 is overkill (12 digits of precision) set_boundary(); $diff = calc_diff(); print "error^2: $diff\r"; for my $i (0 .. $h) { for my $j (0 .. $w) { $v[$i][$j] -= $d[$i][$j]; } } } print "\n";
my @current = (0) x 3; for my $i (0 .. $h) { for my $j (0 .. $w) { $current[ $f[$i][$j] ] += $d[$i][$j] * scalar(@{$n[$i][$j]}); } } return ($current[1] - $current[-1]) / 2; }
print "R = @{[2 / iter()]}\n";</lang>
Perl 6
<lang perl6>my $S = 10;
my @fixed;
sub allocmesh ($w, $h) {
gather for ^$h {
take [0 xx $w];
}
}
sub force-fixed(@f) {
@f[1][1] = 1; @f[6][7] = -1;
}
sub force-v(@v) {
@v[1][1] = 1; @v[6][7] = -1;
}
sub calc_diff(@v, @d, Int $w, Int $h) {
my $total = 0; for ^$h X ^$w -> $i, $j { my @neighbors = grep *.defined, @v[$i-1][$j], @v[$i][$j-1], @v[$i+1][$j], @v[$i][$j+1]; my $v = [+] @neighbors; @d[$i][$j] = $v = @v[$i][$j] - $v / +@neighbors; $total += $v * $v unless @fixed[$i][$j]; } return $total;
}
sub iter(@v, Int $w, Int $h) {
my @d = allocmesh($w, $h); my $diff = 1e10; my @cur = 0, 0, 0;
while $diff > 1e-24 { force-v(@v); $diff = calc_diff(@v, @d, $w, $h); for ^$h X ^$w -> $i, $j { @v[$i][$j] -= @d[$i][$j]; } }
for ^$h X ^$w -> $i, $j { @cur[ @fixed[$i][$j] + 1 ] += @d[$i][$j] * (?$i + ?$j + ($i < $h - 1) + ($j < $w - 1)); }
return (@cur[2] - @cur[0]) / 2;
}
my @mesh = allocmesh($S, $S);
@fixed = allocmesh($S, $S); force-fixed(@fixed);
say 2 / iter(@mesh, $S, $S);</lang>
- Output:
1.60899124172989
Python
<lang python>DIFF_THRESHOLD = 1e-40
class Fixed:
FREE = 0 A = 1 B = 2
class Node:
__slots__ = ["voltage", "fixed"] def __init__(self, v=0.0, f=Fixed.FREE): self.voltage = v self.fixed = f
def set_boundary(m):
m[1][1] = Node( 1.0, Fixed.A) m[6][7] = Node(-1.0, Fixed.B)
def calc_difference(m, d):
h = len(m) w = len(m[0]) total = 0.0
for i in xrange(h): for j in xrange(w): v = 0.0 n = 0 if i != 0: v += m[i-1][j].voltage; n += 1 if j != 0: v += m[i][j-1].voltage; n += 1 if i < h-1: v += m[i+1][j].voltage; n += 1 if j < w-1: v += m[i][j+1].voltage; n += 1 v = m[i][j].voltage - v / n
d[i][j].voltage = v if m[i][j].fixed == Fixed.FREE: total += v ** 2 return total
def iter(m):
h = len(m) w = len(m[0]) difference = [[Node() for j in xrange(w)] for i in xrange(h)]
while True: set_boundary(m) # Enforce boundary conditions. if calc_difference(m, difference) < DIFF_THRESHOLD: break for i, di in enumerate(difference): for j, dij in enumerate(di): m[i][j].voltage -= dij.voltage
cur = [0.0] * 3 for i, di in enumerate(difference): for j, dij in enumerate(di): cur[m[i][j].fixed] += (dij.voltage * (bool(i) + bool(j) + (i < h-1) + (j < w-1)))
return (cur[Fixed.A] - cur[Fixed.B]) / 2.0
def main():
w = h = 10 mesh = [[Node() for j in xrange(w)] for i in xrange(h)] print "R = %.16f" % (2 / iter(mesh))
main()</lang>
- Output:
R = 1.6089912417307286
REXX
This version allows specification of the grid size and the location of the A & B points. <lang rexx>/*REXX pgm calculates resistance between any 2 points on a resister grid*/ numeric digits 20 /*use moderate digits (precision)*/ minVal=(1'e-'||(digits()*2)) / 1 /*calculate the threshold min val*/ if 1=='f1'x then ohms = 'ohms' /*EBCDIC machine? Use 'ohms'. */
else ohms = 'ea'x /* ASCII machine? Use Greek Ω.*/
parse arg wide high Arow Acol Brow Bcol . say 'minVal = ' format(minVal,,,,0) ; say say 'resistor mesh is of size: ' wide "wide, " high 'high.' ; say say 'point A is at (row,col): ' Arow","Acol say 'point B is at (row,col): ' Brow","Bcol @.=0; cell.=1
do until $ <= minVal; $=0; v = 0 @.Arow.Acol = +1 ; cell.Arow.Acol = 0 @.Brow.Bcol = -1 ; cell.Brow.Bcol = 0
do i =1 for high; im=i-1; ip=i+1 do j=1 for wide; jm=j-1; jp=j+1; n=0; v=0 if i\==1 then do; v=v+@.im.j; n=n+1; end if j\==1 then do; v=v+@.i.jm; n=n+1; end if i<high then do; v=v+@.ip.j; n=n+1; end if j<wide then do; v=v+@.i.jp; n=n+1; end v=@.i.j-v/n; #.i.j=v; if cell.i.j then $=$+v*v end /*j*/ end /*i*/ do r=1 for High do c=1 for Wide; @.r.c=@.r.c-#.r.c end /*c*/ end /*r*/ end /*until*/
say Acur = #.Arow.Acol * sides(Arow,Acol) Bcur = #.Brow.Bcol * sides(Brow,Bcol) say 'resistance between point A and point B is: ' 4/(Acur-Bcur) ohms exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────sides subroutine────────────────────*/ sides: parse arg i,j; !=0; if i\==1 & i\==high then !=!+2
else !=!+1 if j\==1 & j\==wide then !=!+2 else !=!+1
return !</lang> output when using the input of: 10 10 2 2 7 8
minVal = 1E-40 resistor mesh is of size: 10 wide, 10 high. point A is at (row,col): 2,2 point B is at (row,col): 7,8 resistance between point A and point B is: 1.6089912417307296559 Ω
Tcl
<lang tcl>package require Tcl 8.6; # Or 8.5 with the TclOO package
- This code is structured as a class with a little trivial DSL parser so it is
- easy to change what problem is being worked on.
oo::class create ResistorMesh {
variable forcePoints V fixed w h
constructor {boundaryConditions} {
foreach {op condition} $boundaryConditions { switch $op { size { lassign $condition w h set fixed [lrepeat $h [lrepeat $w 0]] set V [lrepeat $h [lrepeat $w 0.0]] } fixed { lassign $condition j i v lset fixed $i $j [incr ctr] lappend forcePoints $j $i $v } } }
}
method CalculateDifferences {*dV} {
upvar 1 ${*dV} dV set error 0.0 for {set i 0} {$i < $h} {incr i} { for {set j 0} {$j < $w} {incr j} { set v 0.0 set n 0 if {$i} { set v [expr {$v + [lindex $V [expr {$i-1}] $j]}] incr n } if {$j} { set v [expr {$v + [lindex $V $i [expr {$j-1}]]}] incr n } if {$i+1 < $h} { set v [expr {$v + [lindex $V [expr {$i+1}] $j]}] incr n } if {$j+1 < $w} { set v [expr {$v + [lindex $V $i [expr {$j+1}]]}] incr n } lset dV $i $j [set v [expr {[lindex $V $i $j] - $v/$n}]] if {![lindex $fixed $i $j]} { set error [expr {$error + $v**2}] } } } return $error
}
method FindCurrentFixpoint {epsilon} {
set dV [lrepeat $h [lrepeat $w 0.0]] set current {0.0 0.0 0.0} while true { # Enforce the boundary conditions foreach {j i v} $forcePoints { lset V $i $j $v } # Compute the differences and the error set error [my CalculateDifferences dV] # Apply the differences for {set i 0} {$i < $h} {incr i} { for {set j 0} {$j < $w} {incr j} { lset V $i $j [expr { [lindex $V $i $j] - [lindex $dV $i $j]}] } } # Done if the error is small enough if {$error < $epsilon} break } # Compute the currents from the error for {set i 0} {$i < $h} {incr i} { for {set j 0} {$j < $w} {incr j} { lset current [lindex $fixed $i $j] [expr { [lindex $current [lindex $fixed $i $j]] + [lindex $dV $i $j] * (!!$i+!!$j+($i<$h-1)+($j<$w-1))}] } } # Compute the actual current flowing between source and sink return [expr {([lindex $current 1] - [lindex $current 2]) / 2.0}]
}
# Public entry point method solveForResistance Template:Epsilon 1e-24 {
set voltageDifference [expr { [lindex $forcePoints 2] - [lindex $forcePoints 5]}] expr {$voltageDifference / [my FindCurrentFixpoint $epsilon]}
}
}</lang> Setting up and solving this particular problem: <lang tcl>ResistorMesh create mesh {
size {10 10} fixed {1 1 1.0} fixed {6 7 -1.0}
} puts [format "R = %.12g" [mesh solveForResistance]]</lang>
- Output:
R = 1.60899124173
XPL0
<lang XPL0>code real RlRes=46, RlOut=48; def S = 10;
proc SetBoundary(MV, MF); real MV; int MF; [MF(1,1):= 1; MV(1,1):= 1.0;
MF(6,7):= -1; MV(6,7):= -1.0;
];
func real CalcDiff(MV, MF, DV, W, H); real MV; int MF; real DV; int W, H; int I, J, N; real V, Total; [Total:= 0.0; for I:= 0 to H-1 do
for J:= 0 to W-1 do [V:= 0.0; N:= 0; if I then [V:= V + MV(I-1,J); N:= N+1]; if J then [V:= V + MV(I,J-1); N:= N+1]; if I+1 < H then [V:= V + MV(I+1,J); N:= N+1]; if J+1 < W then [V:= V + MV(I,J+1); N:= N+1]; V:= MV(I,J) - V/float(N); DV(I,J):= V; if MF(I,J) = 0 then Total:= Total + V*V; ];
return Total; ];
func real Iter(MV, MF, W, H); real MV; int MF, W, H; real DV, Diff, Cur; int I, J; [DV:= RlRes(W); for I:= 0 to W-1 do DV(I):= RlRes(H); Diff:= 1E10; Cur:= [0.0, 0.0, 0.0]; while Diff > 1E-24 do
[SetBoundary(MV, MF); Diff:= CalcDiff(MV, MF, DV, W, H); for I:= 0 to H-1 do for J:= 0 to W-1 do MV(I,J):= MV(I,J) - DV(I,J); ];
for I:= 0 to H-1 do
for J:= 0 to W-1 do Cur(MF(I,J)+1):= Cur(MF(I,J)+1) + DV(I,J) * float(-(I>0) - (J>0) - (I<H-1) - (J<W-1)); \middle=4; side=3; corner=2
return (Cur(2)-Cur(0))/2.0; ];
real MeshV(S,S); int MeshF(S,S); RlOut(0, 2.0 / Iter(MeshV, MeshF, S, S))</lang>
- Output:
1.60899