Resistor mesh
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>
- 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>
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