Resistor mesh: Difference between revisions

m
→‎{{header|Wren}}: Changed to Wren S/H
(Rename Perl 6 -> Raku, alphabetize, minor clean-up)
m (→‎{{header|Wren}}: Changed to Wren S/H)
 
(22 intermediate revisions by 12 users not shown)
Line 8:
 
;See also:
*   (humor, nerd sniping)   [http://xkcd.com/356/ xkcd.com cartoon] (you can solve that for extra credits)
*   [https://www.paulinternet.nl/?page=resistors An article on how to calculate this and an implementation in Mathematica]
<br><br>
 
 
=={{header|11l}}==
{{trans|Python}}
 
<syntaxhighlight lang="11l">-V DIFF_THRESHOLD = 1e-40
 
T.enum Fixed
FREE
A
B
 
T Node
Float voltage
Fixed fixed
F (v = 0.0, f = Fixed.FREE)
.voltage = v
.fixed = f
 
F set_boundary(&m)
m[1][1] = Node( 1.0, Fixed.A)
m[6][7] = Node(-1.0, Fixed.B)
 
F calc_difference(m, &d)
V h = m.len
V w = m[0].len
V total = 0.0
 
L(i) 0 .< h
L(j) 0 .< w
V v = 0.0
V n = 0
I i != 0 {v += m[i - 1][j].voltage; n++}
I j != 0 {v += m[i][j - 1].voltage; n++}
I i < h-1 {v += m[i + 1][j].voltage; n++}
I j < w-1 {v += m[i][j + 1].voltage; n++}
v = m[i][j].voltage - v / n
d[i][j].voltage = v
I m[i][j].fixed == FREE
total += v ^ 2
R total
 
F iter(&m)
V h = m.len
V w = m[0].len
V difference = [[Node()] * w] * h
 
L
set_boundary(&m)
I calc_difference(m, &difference) < :DIFF_THRESHOLD
L.break
L(di) difference
V i = L.index
L(dij) di
V j = L.index
m[i][j].voltage -= dij.voltage
 
V cur = [0.0] * 3
L(di) difference
V i = L.index
L(dij) di
V j = L.index
cur[Int(m[i][j].fixed)] += (dij.voltage *
(Int(i != 0) + Int(j != 0) + (i < h - 1) + (j < w - 1)))
 
R (cur[Int(Fixed.A)] - cur[Int(Fixed.B)]) / 2.0
 
V w = 10
V h = 10
V mesh = [[Node()] * w] * h
print(‘R = #.16’.format(2 / iter(&mesh)))</syntaxhighlight>
 
{{out}}
<pre>R = 1.6089912417307285</pre>
 
=={{header|Ada}}==
{{trans|C}}
<langsyntaxhighlight Adalang="ada">with Ada.Text_IO; use Ada.Text_IO;
procedure ResistMesh is
H, W : constant Positive := 10;
Line 81 ⟶ 156:
diff := 4.0 / (curA - curB);
IIO.Put (diff, Exp => 0); New_Line;
end ResistMesh;</langsyntaxhighlight>
{{out}}<pre> 1.60899124173073</pre>
 
Line 87 ⟶ 162:
{{trans|Maxima}}
{{works with|BBC BASIC for Windows}}
<langsyntaxhighlight lang="bbcbasic"> INSTALL @lib$+"ARRAYLIB"
*FLOAT 64
@% = &F0F
Line 117 ⟶ 192:
PROC_invert(A())
B() = A().B()
= B(k%, 0)</langsyntaxhighlight>
{{out}}
<pre>Resistance = 1.60899124173071 ohms</pre>
 
=={{header|C}}==
<langsyntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
Line 186 ⟶ 261:
printf("R = %g\n", 2 / iter(mesh, S, S));
return 0;
}</langsyntaxhighlight>
 
=={{header|C sharp|C#}}==
{{trans|Java}}
<langsyntaxhighlight lang="csharp">using System;
using System.Collections.Generic;
 
Line 297 ⟶ 372:
}
}
}</langsyntaxhighlight>
{{out}}
<pre>R = 1.608991241729890</pre>
Line 303 ⟶ 378:
=={{header|C++}}==
{{trans|C#}}
<langsyntaxhighlight lang="cpp">#include <iomanip>
#include <iostream>
#include <vector>
Line 427 ⟶ 502:
 
return 0;
}</langsyntaxhighlight>
{{out}}
<pre>R = 1.60899124172989</pre>
Line 433 ⟶ 508:
=={{header|D}}==
{{trans|C}}
<langsyntaxhighlight lang="d">import std.stdio, std.traits;
 
enum Node.FP differenceThreshold = 1e-40;
Line 510 ⟶ 585:
 
writefln("R = %.19f", 2 / mesh.iter);
}</langsyntaxhighlight>
{{out}}
<pre>R = 1.6089912417307296556</pre>
Line 516 ⟶ 591:
=={{header|ERRE}}==
We'll solve the linear system. We'll write [[wp:Kirchhoff's circuit laws|Kirchhoff's circuit laws]] at each node and search for a voltage distribution that creates a 1A current coming from A exiting in B. The difference of voltage between B and A is then the resistance.
<syntaxhighlight lang="erre">
<lang ERRE>
PROGRAM RESISTENCE_MESH
Line 593 ⟶ 668:
PRINT("Resistence=";ABS(A[A,NN+1]-A[B,NN+1]))
END PROGRAM
</syntaxhighlight>
</lang>
{{out}}
<pre>Nodes 12 68
Line 602 ⟶ 677:
The functions for this have been implemented in Euler already. Thus the following commands solve this problem.
 
<syntaxhighlight lang="euler math toolbox">
<lang Euler Math Toolbox>
>load incidence;
>{u,r}=solvePotentialX(makeRectangleX(10,10),12,68); r,
1.60899124173
</syntaxhighlight>
</lang>
 
The necessary functions in the file incidence.e are as follows. There are versions with full matrices. But the functions listed here use compressed matrices and the conjugate gradient method.
 
<syntaxhighlight lang="text">
function makeRectangleX (n:index,m:index)
## Make the incidence matrix of a rectangle grid in compact form.
Line 655 ⟶ 730:
return {u,2/f}
endfunction
</syntaxhighlight>
</lang>
 
Here is the code for the conjugate gradient method for compressed, sparse matrices from cpx.e.
 
<syntaxhighlight lang="text">
function cgX (H:cpx, b:real column, x0:real column=none, f:index=10)
## Conjugate gradient method to solve Hx=b for compressed H.
Line 699 ⟶ 774:
return x;
endfunction
</syntaxhighlight>
</lang>
 
=={{header|FreeBASIC}}==
<langsyntaxhighlight lang="freebasic">' version 01-07-2018
' compile with: fbc -s console
 
Line 782 ⟶ 857:
Print : Print "hit any key to end program"
Sleep
End</langsyntaxhighlight>
{{out}}
<pre>Nodes a: 12 b: 68
Line 790 ⟶ 865:
=={{header|Go}}==
{{trans|C}}
<langsyntaxhighlight lang="go">package main
 
import "fmt"
Line 892 ⟶ 967:
fmt.Printf("R = %g\n", 2/iter(mesh, S, S))
}
</syntaxhighlight>
</lang>
 
=={{header|Haskell}}==
{{trans|Octave}} All mutations are expressed as monoidal operations.
<langsyntaxhighlight lang="haskell">{-# LANGUAGE ParallelListComp #-}
import Numeric.LinearAlgebra (linearSolve, toDense, (!), flatten)
import Data.Monoid ((<>), Sum(..))
Line 926 ⟶ 1,001:
x `when` p = if p then x else mempty
 
current = toDense [ ((a, 0), -1) , ((b, 0), 1) , ((n^2-1, 0), 0) ]</langsyntaxhighlight>
 
{{Out}}
Line 936 ⟶ 1,011:
We represent the mesh as a [[wp:Ybus matrix|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. (It has to be the first element because A was the first node. And we can "ignore" B because we made everything be relative to B.) Most of the work is defining <code>Y</code> which represents the Ybus.
 
<langsyntaxhighlight Jlang="j">nodes=: 10 10 #: i. 100
nodeA=: 1 1
nodeB=: 6 7
Line 947 ⟶ 1,022:
Yii=: (* =@i.@#) #/.~ {."1 wiring NB. diagonal of Y represents connections to B
Yij=: -1:`(<"1@[)`]}&(+/~ 0*i.1+#ref) wiring NB. off diagonal of Y represents wiring
Y=: _1 _1 }. Yii+Yij</langsyntaxhighlight>
 
Here, the result of <code>nodes conn offset</code> represents all pairs of nodes where we can connect the argument nodes to neighboring nodes at the specified offset, and <code>wiring</code> 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.
Line 953 ⟶ 1,028:
So:
 
<langsyntaxhighlight Jlang="j"> {.{. %. Y
1.60899</langsyntaxhighlight>
 
Or, if we want an exact answer (this is slow), we can assume our resistors are perfect:
 
<langsyntaxhighlight Jlang="j"> {.{.%. x:Y
455859137025721r283319837425200</langsyntaxhighlight>
 
(here, the letter 'r' separates the numerator from the denominator)
Line 965 ⟶ 1,040:
To get a better feel for what the <code>conn</code> operation is doing, here is a small illustration:
 
<langsyntaxhighlight Jlang="j"> 3 3 #: i.9
0 0
0 1
Line 992 ⟶ 1,067:
 
2 1
2 2</langsyntaxhighlight>
 
In other words, each coordinate pair is matched up with the coordinate pair that you would get by adding the offset to the first of the pair. In actual use, we use this four times, with four offsets (two horizontal and two vertical) to get our complete mesh.
Line 998 ⟶ 1,073:
=={{header|Java}}==
{{trans|Kotlin}}
<langsyntaxhighlight Javalang="java">import java.util.ArrayList;
import java.util.List;
 
Line 1,104 ⟶ 1,179:
System.out.printf("R = %.15f", r);
}
}</langsyntaxhighlight>
{{out}}
<pre>R = 1.608991241729889</pre>
 
 
=={{header|JavaScript}}==
Kirchhoff's circuit laws on the resistor mesh are represented as a linear
equation system for the electric potential at each node of the grid.
The linear equation is then solved using the [https://en.wikipedia.org/wiki/Conjugate_gradient_method conjugate gradient method]:
<syntaxhighlight lang=JavaScript>
// Vector addition, scalar multiplication & dot product:
const add = (u, v) => {let i = u.length; while(i--) u[i] += v[i]; return u;};
const sub = (u, v) => {let i = u.length; while(i--) u[i] -= v[i]; return u;};
const mul = (a, u) => {let i = u.length; while(i--) u[i] *= a; return u;};
const dot = (u, v) => {let s = 0, i = u.length; while(i--) s += u[i]*v[i]; return s;};
 
const W = 10, H = 10, A = 11, B = 67;
 
function getAdjacent(node){ // Adjacency lists for square grid
let list = [], x = node % W, y = Math.floor(node / W);
if (x > 0) list.push(node - 1);
if (y > 0) list.push(node - W);
if (x < W - 1) list.push(node + 1);
if (y < H - 1) list.push(node + W);
return list;
}
 
function linOp(u){ // LHS of the linear equation
let v = new Float64Array(W * H);
for(let i = 0; i < v.length; i++){
if ( i === A || i === B ) {
v[i] = u[i];
continue;
}
// For each node other then A, B calculate the net current flow:
for(let j of getAdjacent(i)){
v[i] += (j === A || j === B) ? u[i] : u[i] - u[j];
}
}
return v;
}
 
function getRHS(phiA = 1, phiB = 0){ // RHS of the linear equation
let b = new Float64Array(W * H);
// Setting boundary conditions (electric potential at A and B):
b[A] = phiA;
b[B] = phiB;
for(let j of getAdjacent(A)) b[j] = phiA;
for(let j of getAdjacent(B)) b[j] = phiB;
return b;
}
 
function init(phiA = 1, phiB = 0){ // initialize unknown vector
let u = new Float64Array(W * H);
u[A] = phiA;
u[B] = phiB;
return u;
}
 
function solveLinearSystem(err = 1e-20){ // conjugate gradient solver
 
let b = getRHS();
let u = init();
let r = sub(linOp(u), b);
let p = r;
let e = dot(r,r);
 
while(true){
let Ap = linOp(p);
let alpha = e / dot(p, Ap);
u = sub(u, mul(alpha, p.slice()));
r = sub(linOp(u), b);
let e_new = dot(r,r);
let beta = e_new / e;
 
if(e_new < err) return u;
 
e = e_new;
p = add(r, mul(beta, p));
}
}
 
function getResistance(u){
let curr = 0;
for(let j of getAdjacent(A)) curr += u[A] - u[j];
return 1 / curr;
}
 
let phi = solveLinearSystem();
let res = getResistance(phi);
console.log(`R = ${res} Ohm`);
</syntaxhighlight>
{{out}}
<pre>R = 1.608991241730955 Ohm</pre>
 
=={{header|jq}}==
'''Works with jq and gojq, that is, the C and Go implementations of jq.'''
 
'''Adapted from [[#Wren|Wren]]'''
<syntaxhighlight lang=jq>
# Create a $rows * $columns matrix initialized with the input value
def matrix($rows; $columns):
. as $in
| [range(0;$columns)|$in] as $row
| [range(0;$rows)|$row];
 
def Node($v; $fixed):
{$v, $fixed};
 
# input: a suitable matrix of Nodes
def setBoundary:
.[1][1].v = 1
| .[1][1].fixed = 1
| .[6][7].v = -1
| .[6][7].fixed = -1 ;
 
# input: {d, m} where
# .d and .m are matrices (as produced by matrix(h; w)) of Nodes
# output: {d, m, diff} with d updated
def calcDiff($w; $h):
def adjust($cond; action): if $cond then action | .n += 1 else . end;
 
reduce range(0; $h) as $i (.diff = 0;
reduce range(0; $w) as $j (.;
.v = 0
| .n = 0
| adjust($i > 0; .v += .m[$i-1][$j].v)
| adjust($j > 0; .v += .m[$i][$j-1].v)
| adjust($i + 1 < $h; .v += .m[$i+1][$j].v)
| adjust($j + 1 < $w; .v += .m[$i][$j+1].v)
| .v = .m[$i][$j].v - .v/.n
| .d[$i][$j].v = .v
| if (.m[$i][$j].fixed == 0) then .diff += .v * .v else . end ) ) ;
 
# input: a mesh of width w and height h, i.e. a matrix as prodcued by matrix(h;w)
def iter:
length as $h
| (.[0]|length) as $w
| { m : .,
d : (Node(0;0) | matrix($h; $w)),
cur: [0,0,0],
diff: 1e10 }
| until (.diff <= 1e-24;
.m |= setBoundary
| calcDiff($w; $h)
| reduce range(0;$h) as $i (.;
reduce range(0;$w) as $j (.;
.m[$i][$j].v += (- .d[$i][$j].v) )) )
| reduce range(0; $h) as $i (.;
reduce range(0; $w) as $j (.;
.k = 0
| if ($i != 0) then .k += 1 else . end
| if ($j != 0) then .k += 1 else . end
| if ($i < $h - 1) then .k += 1 else . end
| if ($j < $w - 1) then .k += 1 else . end
| .cur[.m[$i][$j].fixed + 1] += .d[$i][$j].v * .k ))
| (.cur[2] - .cur[0]) / 2 ;
 
def task($S):
def mesh: Node(0; 0) | matrix($S; $S);
(2 / (mesh | iter)) as $r
| "R = \($r) ohms";
 
task(10)
</syntaxhighlight>
{{output}}
<pre>
R = 1.608991241729889 ohms
</pre>
 
 
=={{header|Julia}}==
Line 1,113 ⟶ 1,355:
Because the graph is a rectangular grid, we can in turn write the incidence matrix D in terms of Kronecker products ⊗ (<code>kron</code> in Julia) of "one-dimensional" D<sub>1</sub> matrices (the incidence matrix of a 1d resistor network).
We use Julia's built-in sparse-matrix solvers (based on SuiteSparse) to solve the resulting sparse linear system efficiently
<langsyntaxhighlight lang="julia">N = 10
D1 = speye(N-1,N) - spdiagm(ones(N-1),1,N-1,N)
D = [ kron(D1, speye(N)); kron(speye(N), D1) ]
Line 1,119 ⟶ 1,361:
b = zeros(N^2); b[i], b[j] = 1, -1
v = (D' * D) \ b
v[i] - v[j]</langsyntaxhighlight>
{{out}}
<pre>
Line 1,129 ⟶ 1,371:
=={{header|Kotlin}}==
{{trans|C}}
<langsyntaxhighlight lang="scala">// version 1.1.4-3
 
typealias List2D<T> = List<List<T>>
Line 1,190 ⟶ 1,432:
val r = 2.0 / iter(mesh, S, S)
println("R = $r")
}</langsyntaxhighlight>
 
{{out}}
Line 1,197 ⟶ 1,439:
</pre>
 
=={{header|Mathematica}}/{{header|Wolfram Language}}==
 
{{trans|Maxima}}
{{works with|Mathematica|13.0}}
<lang mathematica>gridresistor[p_, q_, ai_, aj_, bi_, bj_] :=
Use <b>KirchhoffMatrix</b> and <b>DrazinInverse</b> to compute the effective resistance matrix of a graph:
Block[{A, B, k, c, V}, A = ConstantArray[0, {p*q, p*q}];
 
Do[k = (i - 1) q + j;
<syntaxhighlight lang="mathematica">ResistanceMatrix[g_Graph] := With[{n = VertexCount[g], km = KirchhoffMatrix[g]},
If[{i, j} == {ai, aj}, A[[k, k]] = 1, c = 0;
Table[ ReplacePart[ Diagonal[ DrazinInverse[ ReplacePart[km, k -> UnitVector[n, k]]]], k -> 0],
If[1 <= i + 1 <= p && 1 <= j <= q, c++; A[[k, k + q]] = -1];
{k, n}]
If[1 <= i - 1 <= p && 1 <= j <= q, c++; A[[k, k - q]] = -1];
]
If[1 <= i <= p && 1 <= j + 1 <= q, c++; A[[k, k + 1]] = -1];
 
If[1 <= i <= p && 1 <= j - 1 <= q, c++; A[[k, k - 1]] = -1];
rm = ResistanceMatrix[GridGraph[{10, 10}]];
A[[k, k]] = c], {i, p}, {j, q}];
 
B = SparseArray[(k = (bi - 1) q + bj) -> 1, p*q];
N[rm[[12, 68]], 40]</syntaxhighlight>
LinearSolve[A, B][[k]]];
N[gridresistor[10, 10, 2, 2, 8, 7], 40]</lang>
{{Out}}
<pre>1.608991241730729655954495520510088761201</pre>
 
{{works with|Mathematica|8.0}}
Use <b>KirchhoffMatrix</b> and <b>PseudoInverse</b> to compute the effective resistance matrix of a graph to the desired precision:
<syntaxhighlight lang="mathematica">ResistanceMatrix[g_, prec_:$MachinePrecision]:= With[{m = PseudoInverse[N[KirchhoffMatrix[g], prec]]},
Outer[Plus, Diagonal[m], Diagonal[m]] - m - Transpose[m]
]
 
rm = ResistanceMatrix[GridGraph[{10, 10}], 40];
 
rm[[12, 68]]</syntaxhighlight>
{{works with|Mathematica|9.0}}
<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>
{{Out}}
<pre>1.608991241730729655954495520510088761201</pre>
 
=={{header|Maxima}}==
<langsyntaxhighlight 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.
Line 1,288 ⟶ 1,529:
 
bfloat(%), fpprec = 40;
3.89226554090400912102670691601064387507b0</langsyntaxhighlight>
 
=={{header|Modula-2}}==
<langsyntaxhighlight lang="modula2">MODULE ResistorMesh;
FROM RConversions IMPORT RealToStringFixed;
FROM Terminal IMPORT WriteString,WriteLn,ReadChar;
Line 1,401 ⟶ 1,642:
 
ReadChar;
END ResistorMesh.</langsyntaxhighlight>
 
=={{header|Nim}}==
{{trans|Kotlin}}
<syntaxhighlight lang="nim">const S = 10
 
type
 
NodeKind = enum nodeFree, nodeA, nodeB
 
Node = object
v: float
fixed: NodeKind
 
Mesh[H, W: static int] = array[H, array[W, Node]]
 
 
func setBoundary(m: var Mesh) =
m[1][1].v = 1.0
m[1][1].fixed = nodeA
m[6][7].v = -1.0
m[6][7].fixed = nodeB
 
 
func calcDiff[H, W: static int](m,: Mesh[H, W]; d: var Mesh[H, W]): float =
for i in 0..<H:
for j in 0..<W:
var v = 0.0
var n = 0
if i > 0:
v += m[i - 1][j].v
inc n
if j > 0:
v += m[i][j - 1].v
inc n
if i + 1 < m.H:
v += m[i + 1][j].v
inc n
if j + 1 < m.W:
v += m[i][j + 1].v
inc n
v = m[i][j].v - v / n.toFloat
d[i][j].v = v
if m[i][j].fixed == nodeFree:
result += v * v
 
 
func iter[H, W: static int](m: var Mesh[H, W]): float =
var
d: Mesh[H, W]
cur: array[NodeKind, float]
diff = 1e10
 
while diff > 1e-24:
m.setBoundary()
diff = calcDiff(m, d)
for i in 0..<H:
for j in 0..<W:
m[i][j].v -= d[i][j].v
 
for i in 0..<H:
for j in 0..<W:
var k = 0
if i != 0: inc k
if j != 0: inc k
if i < m.H - 1: inc k
if j < m.W - 1: inc k
cur[m[i][j].fixed] += d[i][j].v * k.toFloat
 
result = (cur[nodeA] - cur[nodeB]) / 2
 
 
when isMainModule:
 
var mesh: Mesh[S, S]
let r = 2 / mesh.iter()
echo "R = ", r</syntaxhighlight>
 
{{out}}
<pre>R = 1.608991241729889</pre>
 
=={{header|Octave}}==
We'll solve the linear system. We'll write [[wp:Kirchhoff's circuit laws|Kirchhoff's circuit laws]] at each node and search for a voltage distribution that creates a 1A current coming from A exiting in B. The difference of voltage between B and A is then the resistance.
 
<langsyntaxhighlight lang="octave">N = 10;
NN = N*N;
G = sparse(NN, NN);
Line 1,445 ⟶ 1,765:
VB = voltage( B );
 
full( abs( VA - VB ) )</langsyntaxhighlight>
{{out}}
<pre>ans = 1.6090</pre>
Line 1,451 ⟶ 1,771:
=={{header|Perl}}==
{{trans|C}}
<langsyntaxhighlight lang="perl">use strict;
use warnings;
 
Line 1,510 ⟶ 1,830:
}
 
printf "R = %.6f\n", 2 / iter();</langsyntaxhighlight>
{{out}}
<pre>R = 1.608991</pre>
Line 1,518 ⟶ 1,838:
uses inverse() from [[Gauss-Jordan_matrix_inversion#Phix]]
and matrix_mul() from [[Matrix_multiplication#Phix]]
<!--<syntaxhighlight lang="phix">-->
<lang Phix>function resistormesh(integer ni, nj, ai, aj, bi, bj)
<span style="color: #008080;">function</span> <span style="color: #000000;">resistormesh</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">ni</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">nj</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">ai</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">aj</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">bi</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">bj</span><span style="color: #0000FF;">)</span>
integer n = ni*nj, k, c
<span style="color: #004080;">integer</span> <span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">ni</span><span style="color: #0000FF;">*</span><span style="color: #000000;">nj</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">k</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">c</span>
sequence A = repeat(repeat(0,n),n),
<span style="color: #004080;">sequence</span> <span style="color: #000000;">A</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span><span style="color: #000000;">n</span><span style="color: #0000FF;">),</span>
B = repeat({0},n)
<span style="color: #000000;">B</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">({</span><span style="color: #000000;">0</span><span style="color: #0000FF;">},</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
for i=1 to ni do
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">ni</span> <span style="color: #008080;">do</span>
for j=1 to nj do
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">nj</span> <span style="color: #008080;">do</span>
k = (i-1)*nj + j--1
<span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">nj</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">j</span><span style="color: #000080;font-style:italic;">--1</span>
if i=ai and j=aj then
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">ai</span> <span style="color: #008080;">and</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">aj</span> <span style="color: #008080;">then</span>
A[k,k] = 1
<span style="color: #000000;">A</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
else
<span c style="color: 0#008080;">else</span>
if i<nispan thenstyle="color: #000000;">c</span> +<span style="color: 1#0000FF;">=</span> A[k,k+nj]<span style="color: -1 end if#000000;">0</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;"><</span><span style="color: #000000;">ni</span> <span style="color: #008080;">then</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">;</span> <span style="color: #000000;">A</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">nj</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
if i>1 then c += 1; A[k,k-nj] = -1 end if
<span style="color: #008080;">if</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">;</span> <span style="color: #000000;">A</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">nj</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
if j<nj then c += 1; A[k,k+1] = -1 end if
<span style="color: #008080;">if</span> <span style="color: #000000;">j</span><span style="color: #0000FF;"><</span><span style="color: #000000;">nj</span> <span style="color: #008080;">then</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">;</span> <span style="color: #000000;">A</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
if j>1 then c += 1; A[k,k-1] = -1 end if
<span style="color: #008080;">if</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">></span><span style="color: #000000;">1</span> <span style="color: #008080;">then</span> <span style="color: #000000;">c</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span><span style="color: #0000FF;">;</span> <span style="color: #000000;">A</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">-</span><span style="color: #000000;">1</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
A[k,k] = c
<span style="color: #000000;">A</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">k</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">c</span>
end if
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
end for
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
k = (bi-1)*nj +bj
<span style="color: #000000;">k</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">bi</span><span style="color: #0000FF;">-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">nj</span> <span style="color: #0000FF;">+</span><span style="color: #000000;">bj</span>
B[k,1] = 1
<span style="color: #000000;">B</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
A = inverse(A)
<span style="color: #000000;">A</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">inverse</span><span style="color: #0000FF;">(</span><span style="color: #000000;">A</span><span style="color: #0000FF;">)</span>
B = matrix_mul(A,B)
<span style="color: #000000;">B</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">matrix_mul</span><span style="color: #0000FF;">(</span><span style="color: #000000;">A</span><span style="color: #0000FF;">,</span><span style="color: #000000;">B</span><span style="color: #0000FF;">)</span>
return B[k,1]
<span style="color: #008080;">return</span> <span style="color: #000000;">B</span><span style="color: #0000FF;">[</span><span style="color: #000000;">k</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span>
end function
<span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
 
printf(1,"Resistance = %.13f ohms\n",resistormesh(10, 10, 2, 2, 8, 7))</lang>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Resistance = %.13f ohms\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">resistormesh</span><span style="color: #0000FF;">(</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">10</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">8</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">7</span><span style="color: #0000FF;">))</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Line 1,552 ⟶ 1,874:
=={{header|Python}}==
{{trans|D}}
<langsyntaxhighlight lang="python">DIFF_THRESHOLD = 1e-40
 
class Fixed:
Line 1,615 ⟶ 1,937:
print "R = %.16f" % (2 / iter(mesh))
 
main()</langsyntaxhighlight>
{{out}}
<pre>R = 1.6089912417307286</pre>
Line 1,621 ⟶ 1,943:
{{trans|Maxima}}
 
<langsyntaxhighlight lang="python">import sys, copy
from fractions import Fraction
 
Line 1,691 ⟶ 2,013:
# python grid.py 10 10 1 1 7 6
# 455859137025721/283319837425200
# 1.6089912417307297</langsyntaxhighlight>
 
=={{header|Racket}}==
Line 1,698 ⟶ 2,020:
This version avoids mutation... possibly a little more costly than C, but more functional.
 
<langsyntaxhighlight lang="racket">#lang racket
(require racket/flonum)
 
Line 1,767 ⟶ 2,089:
 
(module+ main
(printf "R = ~a~%" (mesh-R '((1 1) (6 7)) 10 10)))</langsyntaxhighlight>
 
{{out}}
Line 1,774 ⟶ 2,096:
=={{header|Raku}}==
(formerly Perl 6)
{{trans|cC}}
<syntaxhighlight lang="raku" perl6line>my $S*TOLERANCE = 101e-12;
 
sub set-boundary(@mesh,@p1,@p2) {
my @fixed;
@mesh[ @p1[0] ; @p1[1] ] = 1;
 
@mesh[ @p2[0] ; @p2[1] ] = -1;
sub allocmesh ($w, $h) {
gather for ^$h {
take [0 xx $w];
}
}
 
sub force-fixedsolve(@fp1, @p2, Int \w, Int \h) {
my @f[1][1]d = [0 xx w] xx 1h;
my @fV = [6][70 xx w] =xx -1h;
my @fixed = [0 xx w] xx h;
}
set-boundary(@fixed,@p1,@p2);
 
loop {
sub force-v(@v) {
set-boundary(@V,@p1,@p2);
@v[1][1] = 1;
@v[6][7] my $diff = -10;
for (flat ^h X ^w) -> \i, \j {
}
my @neighbors = (@V[i-1;j], @V[i;j-1], @V[i+1;j], @V[i;j+1]).grep: *.defined;
@d[i;j] = my \v = @V[i;j] - @neighbors.sum / @neighbors;
sub calc_diff(@v, @d, Int $w, Int $h) {
my $totaldiff += 0v × v unless @fixed[i;j];
for (flat ^$h X ^$w) -> $i, $j {}
last if $diff =~= 0;
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 for (flat ^h X ^w) -24> \i, \j {
force @V[i;j] -v(= @v)d[i;j];
$diff = calc_diff(@v, @d, $w, $h);
for (flat ^$h X ^$w) -> $i, $j {
@v[$i][$j] -= @d[$i][$j];
}
}
 
my @current;
for (flat ^$h X ^$w) -> $i, $j {
for (flat ^h X @cur[^w) @fixed[$-> \i][$j], + 1\j ]{
@current[ @fixed[i;j]+1 ] += @d[$i][$;j] *× (?$i + ?$j + ($i < $h - 1) + ($j < $w - 1) );
}
(@current[2] - @current[0]) / 2
 
return (@cur[2] - @cur[0]) / 2;
}
my @mesh = allocmesh($S, $S);
 
@fixed = allocmesh($S, $S);
force-fixed(@fixed);
 
say 2 / itersolve (@mesh1, $S1), $S(6,7), 10, 10;</lang>
</syntaxhighlight>
{{out}}
<pre>1.60899124172989</pre>
Line 1,841 ⟶ 2,142:
 
Dropping the decimal digits precision &nbsp; ('''numeric digits''') &nbsp; to &nbsp; '''10''' &nbsp; makes the execution &nbsp; '''3''' &nbsp; times faster.
<langsyntaxhighlight lang="rexx">/*REXX program calculates the resistance between any two points on a resistor grid.*/
if 2=='f2'x then ohms = "ohms" /*EBCDIC machine? Then use 'ohms'. */
else ohms = "Ω" /* ASCII " " " Greek Ω.*/
Line 1,885 ⟶ 2,186:
sides: parse arg row,col; z=0; if row\==1 & row\==high then z= z+2; else z= z+1
if col\==1 & col\==wide then z= z+2; else z= z+1
return z</langsyntaxhighlight>
{{out|output|text=&nbsp; when using the default inputs:}}
<pre>
Line 1,900 ⟶ 2,201:
=={{header|Sidef}}==
{{trans|Perl}}
<langsyntaxhighlight lang="ruby">var (w, h) = (10, 10)
 
var v = h.of { w.of(0) } # voltage
Line 1,949 ⟶ 2,250:
}
 
say "R = #{2 / iter()}"</langsyntaxhighlight>
{{out}}
<pre>
Line 1,957 ⟶ 2,258:
=={{header|Tcl}}==
{{trans|C}}
<langsyntaxhighlight 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
Line 2,051 ⟶ 2,352:
expr {$voltageDifference / [my FindCurrentFixpoint $epsilon]}
}
}</langsyntaxhighlight>
Setting up and solving this particular problem:
<langsyntaxhighlight lang="tcl">ResistorMesh create mesh {
size {10 10}
fixed {1 1 1.0}
fixed {6 7 -1.0}
}
puts [format "R = %.12g" [mesh solveForResistance]]</langsyntaxhighlight>
{{out}}
<pre>
R = 1.60899124173
</pre>
 
=={{header|Wren}}==
{{trans|Kotlin}}
<syntaxhighlight lang="wren">class Node {
construct new(v, fixed) {
_v = v
_fixed = fixed
}
 
v { _v }
v=(value) { _v = value }
 
fixed { _fixed }
fixed=(value) { _fixed = value }
}
 
var setBoundary = Fn.new { |m|
m[1][1].v = 1
m[1][1].fixed = 1
m[6][7].v = -1
m[6][7].fixed = -1
}
 
var calcDiff = Fn.new { |m, d, w, h|
var total = 0
for (i in 0...h) {
for (j in 0...w) {
var v = 0
var n = 0
if (i > 0) {
v = v + m[i-1][j].v
n = n + 1
}
if (j > 0) {
v = v + m[i][j-1].v
n = n + 1
}
if (i + 1 < h) {
v = v + m[i+1][j].v
n = n + 1
}
if (j + 1 < w) {
v = v + m[i][j+1].v
n = n + 1
}
v = m[i][j].v - v/n
d[i][j].v = v
if (m[i][j].fixed == 0) total = total + v*v
}
}
return total
}
 
var iter = Fn.new { |m, w, h|
var d = List.filled(h, null)
for (i in 0...h) {
d[i] = List.filled(w, null)
for (j in 0...w) d[i][j] = Node.new(0, 0)
}
var cur = [0] * 3
var diff = 1e10
while (diff > 1e-24) {
setBoundary.call(m)
diff = calcDiff.call(m, d, w, h)
for (i in 0...h) {
for (j in 0...w) m[i][j].v = m[i][j].v - d[i][j].v
}
}
for (i in 0...h) {
for (j in 0...w) {
var k = 0
if (i != 0) k = k + 1
if (j != 0) k = k + 1
if (i < h - 1) k = k + 1
if (j < w - 1) k = k + 1
cur[m[i][j].fixed + 1] = cur[m[i][j].fixed + 1] + d[i][j].v * k
}
}
return (cur[2] - cur[0]) / 2
}
 
var S = 10
var mesh = List.filled(S, null)
for (i in 0...S) {
mesh[i] = List.filled(S, null)
for (j in 0...S) mesh[i][j] = Node.new(0, 0)
}
var r = 2 / iter.call(mesh, S, S)
System.print("R = %(r)")</syntaxhighlight>
 
{{out}}
<pre>
R = 1.6089912417299
</pre>
 
=={{header|XPL0}}==
{{trans|C}}
<langsyntaxhighlight XPL0lang="xpl0">code real RlRes=46, RlOut=48;
def S = 10;
 
Line 2,114 ⟶ 2,509:
real MeshV(S,S); int MeshF(S,S);
RlOut(0, 2.0 / Iter(MeshV, MeshF, S, S))</langsyntaxhighlight>
 
{{out}}
Line 2,123 ⟶ 2,518:
=={{header|Yabasic}}==
{{trans|ERRE}}
<langsyntaxhighlight Yabasiclang="yabasic"> N=10
NN=N*N
DIM A(NN,NN+1)
Line 2,188 ⟶ 2,583:
NEXT
NEXT
PRINT "Resistence = "; : PRINT ABS(A(A,NN+1)-A(B,NN+1)) USING "%1.13f"</langsyntaxhighlight>
 
=={{header|zkl}}==
{{trans|Maxima}}
Uses the GNU Scientific Library.
<langsyntaxhighlight lang="zkl">var [const] GSL=Import("zklGSL"); // libGSL (GNU Scientific Library)
fcn onGrid(i,j,p,q){ ((0<=i<p) and (0<=j<q)) }
fcn gridResistor(p,q, ai,aj, bi,bj){
Line 2,212 ⟶ 2,607:
b[k:=bi*q + bj]=1;
A.AxEQb(b)[k];
}</langsyntaxhighlight>
<langsyntaxhighlight lang="zkl">gridResistor(10,10, 1,1, 7,6).println();</langsyntaxhighlight>
{{out}}
<pre>
9,476

edits