Resistor mesh
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 Fixed : size_t { free = 0, A = 1, B = 2 }
FP voltage = 0.0; private Fixed fixed__ = Fixed.free;
this(in FP v, in Fixed f) pure nothrow {
this.voltage = v;
this.fixed__ = f;
}
@property Fixed fixed() const pure nothrow { return fixed__; }
}
Node.FP iter(size_t w, size_t h)(ref Node[w][h] m) pure nothrow {
static void enforceBoundaryConditions(size_t w, size_t h)
(ref Node[w][h] m)
pure nothrow {
m[1][1].voltage = 1.0;
m[6][7].voltage = -1.0;
}
static Node.FP calcDifference(size_t w, size_t h)
(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;
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.Fixed.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.Fixed).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.Fixed.A] - cur[Node.Fixed.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.Fixed.A); mesh[6][7] = Node(-1.0, Node.Fixed.B);
writefln("R = %.19f", 2 / iter(mesh));
}</lang>
- Output:
R = 1.6089912417307296556
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 ; say @.=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 /*r*/
end /*c*/
end /*until $<=minVal*/
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