问题描述
我正在尝试使用拉普拉斯矩阵求解网络上的流。我首先在这里测试这个问题：https://rosettacode.org/wiki/Resistor_mesh#Python
它完美地给出了解决方案，当所有权重都为 1 时，R = 1.6089 无论如何。我希望能够解决不等于 1 的电阻器！然后能够获得每个电阻上的电流，所以我尝试为电导（权重）生成一组随机值并检查它是否正常工作，看到来自注入电流的节点的电流总和相等到注入的电流。不幸的是，事实并非如此。我已经检查了拉普拉斯算子，它看起来和我预期的一样，但除此之外我完全迷失了，有人能解释一下我是否在这里找到了正确的拉普拉斯算子吗？或者如果我遗漏了一些非常明显的东西，例如一些不正确的索引？
代码如下（请原谅我生成网格的非常业余的方式，任何提示和技巧也值得赞赏！）：
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
问题在于计算出的电流和预期的不一样：
> neighbours_currents
[1] 0.05035087 0.09044874 0.03309549 0.21782845
解决方法
答案：函数laplacian_matrix()按降序对节点重新排序，这意味着所有节点都被打乱了。
解决方法是指定节点标签的顺序
the_graph <- graph_from_data_frame(Edges,directed = FALSE)
变成了
Verts <- data.frame(label = 1:(N_x*N_y))
the_graph <- graph_from_data_frame(Edges,directed = FALSE,vertices = Verts)
求和为 1 的电流：
> 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