Resistor mesh: Difference between revisions

;Task:

Given &nbsp; <big> 10&times;10 </big> &nbsp; grid nodes &nbsp; (as shown in the image) &nbsp; interconnected by &nbsp; <big> 1Ω </big> &nbsp; resistors as shown,

 

 

;See also:

<br><br>

 

=={{header|Ada}}==

{{trans|C}}

procedure ResistMesh is

H, W : constant Positive := 10;

diff := 4.0 / (curA - curB);

IIO.Put (diff, Exp => 0); New_Line;

{{out}}<pre> 1.60899124173073</pre>

 

{{trans|Maxima}}

{{works with|BBC BASIC for Windows}}

*FLOAT 64

@% = &F0F

PROC_invert(A())

B() = A().B()

{{out}}

<pre>Resistance = 1.60899124173071 ohms</pre>

 

=={{header|C}}==

#include <stdlib.h>

printf("R = %g\n", 2 / iter(mesh, S, S));

return 0;

 

=={{header|D}}==

{{trans|C}}

 

enum Node.FP differenceThreshold = 1e-40;

 

writefln("R = %.19f", 2 / mesh.iter);

{{out}}

<pre>R = 1.6089912417307296556</pre>

 

=={{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.

PROGRAM RESISTENCE_MESH

PRINT("Resistence=";ABS(A[A,NN+1]-A[B,NN+1]))

END PROGRAM

{{out}}

<pre>Nodes 12 68

The functions for this have been implemented in Euler already. Thus the following commands solve this problem.

 

>load incidence;

>{u,r}=solvePotentialX(makeRectangleX(10,10),12,68); r,

1.60899124173

 

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.

 

function makeRectangleX (n:index,m:index)

## Make the incidence matrix of a rectangle grid in compact form.

return {u,2/f}

endfunction

 

Here is the code for the conjugate gradient method for compressed, sparse matrices from cpx.e.

 

function cgX (H:cpx, b:real column, x0:real column=none, f:index=10)

## Conjugate gradient method to solve Hx=b for compressed H.

return x;

endfunction

 

=={{header|Go}}==

{{trans|C}}

 

import "fmt"

fmt.Printf("R = %g\n", 2/iter(mesh, S, S))

}

 

=={{header|Haskell}}==

{{trans|Octave}} All mutations are expressed as monoidal operations.

import Numeric.LinearAlgebra (linearSolve, toDense, (!), flatten)

import Data.Monoid ((<>), Sum(..))

x `when` p = if p then x else mempty

 

 

{{Out}}

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.

 

nodeA=: 1 1

nodeB=: 6 7

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

 

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.

So:

 

 

Or, if we want an exact answer (this is slow), we can assume our resistors are perfect:

 

 

(here, the letter 'r' separates the numerator from the denominator)

To get a better feel for what the <code>conn</code> operation is doing, here is a small illustration:

 

0 0

0 1

 

2 1

 

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.

 

=={{header|Julia}}==

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₁ 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

D1 = speye(N-1,N) - spdiagm(ones(N-1),1,N-1,N)

D = [ kron(D1, speye(N)); kron(speye(N), D1) ]

b = zeros(N^2); b[i], b[j] = 1, -1

v = (D' * D) \ b

{{out}}

<pre>

=={{header|Kotlin}}==

{{trans|C}}

 

typealias List2D<T> = List<List<T>>

val r = 2.0 / iter(mesh, S, S)

println("R = $r")

 

{{out}}

R = 1.608991241729889

</pre>

 

=={{header|Maxima}}==

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.

is known up to a constant. For A we will simply write the equation V(A) = 0, to

keep the program simple.

 

bfloat(%), fpprec = 40;

 

 

 

 

=={{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.

 

NN = N*N;

G = sparse(NN, NN);

VB = voltage( B );

 

{{out}}

<pre>ans = 1.6090</pre>

=={{header|Perl}}==

{{trans|C}}

 

my ($w, $h) = (9, 9);

sub iter {

my $diff = 1;

set_boundary();

$diff = calc_diff();

for my $i (0 .. $h) {

for my $j (0 .. $w) {

}

}

 

my @current = (0) x 3;

}

 

 

{{out}}

 

=={{header|Python}}==

{{trans|D}}

 

class Fixed:

print "R = %.16f" % (2 / iter(mesh))

 

{{out}}

<pre>R = 1.6089912417307286</pre>

{{trans|Maxima}}

 

from fractions import Fraction

 

# python grid.py 10 10 1 1 7 6

# 455859137025721/283319837425200

 

=={{header|Racket}}==

This version avoids mutation... possibly a little more costly than C, but more functional.

 

(require racket/flonum)

 

 

(module+ main

 

{{out}}

<pre>R = 1.6089912417301238</pre>

 

=={{header|REXX}}==

This version allows specification of the grid size, &nbsp; the locations of the &nbsp; '''A''' &nbsp; and &nbsp; '''B''' &nbsp; points, &nbsp; and the number of decimal digits precision.

 

if 2=='f2'x then ohms = "ohms" /*EBCDIC machine? Then use 'ohms'. */

else ohms = "Ω" /* ASCII " " " Greek Ω.*/

parse arg high wide Arow Acol Brow Bcol digs . /*obtain optional arguments from the CL*/

numeric digits digs /*use moderate decimal digs (precision)*/

say ' minimum value is ' format(minVal,,,,0) " using " digs ' decimal digits'; say

say ' point A is at (row,col): ' Arow"," Acol

say ' point B is at (row,col): ' Brow"," Bcol

say

say ' resistance between point A and point B is: ' 4 / (Acur - Bcur) ohms

exit /*stick a fork in it, we're all done. */

/*──────────────────────────────────────────────────────────────────────────────────────*/

{{out|output|text=&nbsp; when using the default inputs:}}

<pre>

=={{header|Sidef}}==

{{trans|Perl}}

 

var v = h.of { w.of(0) } # voltage

}

 

{{out}}

<pre>

=={{header|Tcl}}==

{{trans|C}}

 

# This code is structured as a class with a little trivial DSL parser

expr {$voltageDifference / [my FindCurrentFixpoint $epsilon]}

}

Setting up and solving this particular problem:

size {10 10}

fixed {1 1 1.0}

fixed {6 7 -1.0}

}

{{out}}

<pre>

R = 1.60899124173

</pre>

 

=={{header|XPL0}}==

{{trans|C}}

def S = 10;

 

real MeshV(S,S); int MeshF(S,S);

 

{{out}}

1.60899

</pre>

 

=={{header|zkl}}==

{{trans|Maxima}}

Uses the GNU Scientific Library.

fcn onGrid(i,j,p,q){ ((0<=i<p) and (0<=j<q)) }

fcn gridResistor(p,q, ai,aj, bi,bj){

b[k:=bi*q + bj]=1;

A.AxEQb(b)[k];

{{out}}

<pre>