I'm trying to solve for flows on a network using a laplacian matrix. I started out by testing on this problem here: https://rosettacode.org/wiki/Resistor_mesh#Python
And it comes out with the solution perfectly, when all weights are 1, R = 1.6089whatever. I want to be able to solve for resistors that don't equal one though! then be able to get the currents on each resistor, so I tried producing a set of random values for the conductances (weights) and to check it was working see that the sum of the currents coming from the node where the current is injected is equal to the current injected as it should be. Unfortunately it was not. I've inspected the Laplacian and it looks about how I was expecting, but beyond that I'm totally lost, can anyone shed any light on whether I've even got the right laplacian here? or if I'm missing something really obvious such as some incorrect indexing?
Code is below (excuse the very amateur way I've generated the grid, any tips and tricks on that appreciated too!):
library(igraph)
library(SparseM)
# setup edgelist + grid to look at
sample_grid = function(N_x, N_y, xlim = c(-1,1), ylim = c(-1,1)) {
# bounding box
min_x <- xlim[1]
max_x <- xlim[2]
min_y <- ylim[1]
max_y <- ylim[2]
x_locs <- seq(0, N_x)
y_locs <- seq(0, N_y)
N_grid <- N_x * N_y
x <- rep(0, length(N_grid))
y <- rep(0, length(N_grid))
for(i in 1:N_x){
for(j in 1:N_y) {
x[N_y * (i-1) + j] <- x_locs[i]
y[N_y * (i-1) + j] <- y_locs[j]
}
}
locations <- cbind(x,y)
return(locations)
}
N_x <- 10
N_y <- 10
# node locations
Vert <- sample_grid(N_x,N_y)
# edges...
# horizontally...
fromH <- c()
for(i in 1:(N_x-1)){
fromH <- c(fromH, seq(0,(N_y*N_x-N_x), length.out = (N_y)) + i)
}
toH <- fromH + 1
# vertically...
fromV <- 1:(N_x*N_y-N_x)
toV <- fromV + N_x
# --------------------------------------------------------------------------------------------
# crux
Edges <- data.frame(from = c(fromH, fromV), to = c(toH, toV)) # , weights = rep(1,(2*N_x*N_y - (N_x+N_y))))
# change the weights up
set.seed(1)
weight <- rlnorm(n = nrow(Edges), meanlog = 0, sdlog = 1)
Edges$weight <- weight
the_graph <- graph_from_data_frame(Edges, directed = FALSE)
lo <- layout.norm(as.matrix(Vert))
plot(the_graph, layout = lo, directed = FALSE, edge.arrow.size=0)
# solving
L <- laplacian_matrix(the_graph, weights = weight)
# boundary conditions on 68, 12:
# draw 1 amp @ 12
# inject 1 amp @ 68
q <- matrix(rep(0, nrow(Vert)), ncol = 1)
q[68,] <- +1
q[12,] <- -1
# solve
p <- solve(L,q)
R <- p[68,] - p[12,]
# investigating why the weights aren't working! (wrong first attempt)
# neighbours <- c(78,58,67,69) # neighbours of node 68
# neighbours_weights <- weight[neighbours]
# neighbours_potential_diffs <- p[68,] - p[neighbours,]
# neighbours_currents <- neighbours_potential_diffs * neighbours_weights
neighbours <- which(Edges$from == 68 | Edges$to == 68)
potential_diffs <- p[Edges$from,] - p[Edges$to,]
currents <- potential_diffs * Edges$weight
what <- cbind(Edges,p[Edges$from,], p[Edges$to,], currents)
what[neighbours,]
# exact solution for effective resistance between 68 and 12 with 10x10 and all 1ohm
exact <- 455859137025721/283319837425200
The problem being that the calculated currents do not sum to one as expected:
> neighbours_currents
[1] 0.05035087 0.09044874 0.03309549 0.21782845
Answer: the function laplacian_matrix() reorders the nodes in descending degree, which meant that all the nodes were jumbled up.
The fix was to specify the sequence of node labels so
the_graph <- graph_from_data_frame(Edges, directed = FALSE)
became
Verts <- data.frame(label = 1:(N_x*N_y))
the_graph <- graph_from_data_frame(Edges, directed = FALSE, vertices = Verts)
giving currents that sum to one:
> what[neighbours,]
from to weight p[Edges$from, ] p[Edges$to, ] currents
67 67 68 0.1644813 0.2582772 0.8279230 -0.09369607
77 68 69 0.6419198 0.8279230 0.4943076 0.21415437
148 58 68 1.0175478 0.3902397 0.8279230 -0.44536367
158 68 78 0.5372635 0.8279230 0.3685843 0.24678589