使用Mathematica，在任意有限的电阻网络中计算两点之间等效电阻

代码来自 https://www.paulinternet.nl/downloads/resistors.pdf

voltagesOneAmp[graph_, vertexIn_, vertexOut_] := Module[
{g = graph, indexIn, indexOut, m, i}
,
indexIn = VertexIndex[g, vertexIn];
indexOut = VertexIndex[g, vertexOut];
If[WeightedGraphQ[g],
AnnotationValue[g, EdgeWeight] = 1 / AnnotationValue[g, EdgeWeight]
];
m = WeightedAdjacencyMatrix[g];
m = DiagonalMatrix[Total[m]] - m;
m = Drop[m, {indexOut}, {indexOut}];
i = Drop[SparseArray[{indexIn -> 1}, VertexCount[g]], {indexOut}];
Insert[LinearSolve[m, i], 0, indexOut]
];
equivalentResistance[graph_, vertex1_, vertex2_] :=
voltagesOneAmp[graph, vertex1, vertex2][[VertexIndex[graph, vertex1]]];

showExample[edges_, weights_, from_, to_] := Module[{e, g, r},
e = UndirectedEdge @@ # & /@ edges;
g = EdgeTaggedGraph[e, EdgeWeight -> weights, EdgeLabels -> "EdgeWeight",
VertexLabels -> Automatic, ImageSize -> Small];
r = equivalentResistance[g, from, to] // Simplify;
Print[g];
Print[r];
];

showExample[{a.b, b.c}, {r1, r2}, a, c];
showExample[{a.b, a.c, b.c, b.d, c.d}, {rab, rac, rbc, rbd, rcd}, a, d];

a = 19;
b = 74;
g = GridGraph[{10, 10}];
v = voltagesOneAmp[g, a, b];
r = v[[a]];
colors = ColorData["TemperatureMap"];
Graph[g,
VertexLabels -> {a -> Placed["A", Center], b -> "B"},
VertexSize -> MapThread[#1 -> #2/r &, {VertexList[g], v}],
VertexStyle -> MapThread[#1 -> colors[#2/r] &, {VertexList[g], v}]
]
TildeTilde[r, N[r]]

Clear["Global`*"];
resistanceMatrix[graph_] := 
  Module[{g = graph, m}, 
   If[WeightedGraphQ[g], 
    AnnotationValue[g, EdgeWeight] = 1/AnnotationValue[g, EdgeWeight]];
   m = WeightedAdjacencyMatrix[g];
   m = DiagonalMatrix[Total[m]] - m;
   Table[Insert[Drop[m, {i}, {i}] // Inverse // Diagonal, 0, i], {i, 
     VertexCount[g]}]];

edges = {a . b, a . c, b . c, b . d, c . d};
e = UndirectedEdge @@ # & /@ edges;
weights = {rab, rac, rbc, rbd, rcd};
g = EdgeTaggedGraph[e, EdgeWeight -> weights, 
  EdgeLabels -> "EdgeWeight", VertexLabels -> Automatic, 
  ImageSize -> Small]
resistanceMatrix[g] // MatrixForm

解释resistanceMatrix这个函数：

Let's break it down:

Module[{g = graph, m}, ...]: This initializes a local scope where g is initialized to the input graph and m is a variable to be defined later.

If[WeightedGraphQ[g], AnnotationValue[g, EdgeWeight] = 1/AnnotationValue[g, EdgeWeight]];: This checks if the graph g is a weighted graph. If it is, the weights of the edges are inverted (i.e., each weight is replaced with its reciprocal).

m = WeightedAdjacencyMatrix[g];: This creates a weighted adjacency matrix for the graph g and assigns it to m. In a weighted adjacency matrix, the element at row i and column j is the weight of the edge connecting node i and node j.

m = DiagonalMatrix[Total[m]] - m;: Here, Total[m] sums up the elements of each row of m, giving a list of the sums. DiagonalMatrix[Total[m]] creates a diagonal matrix with the sums on the diagonal. Subtracting m from this diagonal matrix gives a matrix that is the Laplacian matrix of the graph. The Laplacian matrix is defined as the degree matrix (a diagonal matrix containing the degree of each vertex) minus the adjacency matrix.

Table[Insert[Drop[m, {i}, {i}] // Inverse // Diagonal, 0, i], {i, VertexCount[g]}]: This is the main calculation of the resistance matrix. For each vertex i (from 1 to the total number of vertices), it does the following:

Drop[m, {i}, {i}]: Drops the ith row and column from the Laplacian matrix.
Inverse: Calculates the inverse of the resulting matrix.
Diagonal: Extracts the diagonal elements of the inverted matrix.
Insert[..., 0, i]: Inserts a zero at the ith position in the diagonal elements.
This is repeated for each vertex i, creating a list of lists (a matrix) where the ith list (row) is the ith row of the resistance matrix.

In the context of electrical networks, the resistance matrix is a matrix where the i, j entry is the effective resistance between node i and node j when they are viewed as nodes in an electrical network with resistors along the edges of the graph.