Resistor mesh

Resistor mesh is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Given 10 × 10 grid nodes interconnected by 1Ω resistors as shown, find the resistance between point A and B.

See also [[1]]

C

<lang c>#include <stdio.h>

1. include <stdlib.h>

1. define S 10

typedef struct { double v; int fixed; } node;

1. 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

Translation of: C

<lang d>import std.traits: EnumMembers;

enum Node.FP differenceThreshold = 1e-40;

bool isRectangular(T)(in T[][] mat) pure nothrow {

   foreach (const(T[]) row; mat)
       if (row.length != mat[0].length)
           return false;
   return true;

}

struct Node {

   alias real FP;
   FP voltage = 0.0;
   enum Fixed : size_t { free = 0, A = 1, B = 2 }
   /*const*/ Fixed fixed = Fixed.free;

}

void setBoundary(Node[][] m) pure nothrow {

   m[1][1] = Node( 1.0, Node.Fixed.A);
   m[6][7] = Node(-1.0, Node.Fixed.B);

}

Node.FP calcDifference(in Node[][] m, Node[][] d) pure nothrow in {

   assert(isRectangular(m) && isRectangular(d));
   // Are equal sized?
   assert(m.length == d.length && m[0].length == d[0].length);

} body {

   immutable h = m.length;
   immutable w = m[0].length;
   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.Fixed.free)
               total += v ^^ 2;
       }

   return total;

}

Node.FP iter(Node[][] m) pure nothrow in {

   assert(isRectangular(m));

} body {

   immutable h = m.length;
   immutable w = m[0].length;
   auto difference = new Node[][](h, w);

   while (true) {
       setBoundary(m); // Enforce boundary conditions.
       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[Node.Fixed] cur = 0.0;
   alias EnumMembers!(Node.Fixed) FixedFields;
   enum size_t nFixed = FixedFields.length;
   Node.FP[nFixed] 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() {

   import std.stdio;
   enum size_t w = 10,
               h = 10;
   auto mesh = new Node[][](h, w);
   writefln("R = %.19f", 2 / iter(mesh));

}</lang> Output:

R = 1.6089912417307296556

Perl

Translation of: C

<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

Translation of: c

<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

Translation of: D

<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