Блог

Главная - Вопросы по программированию - Матрица Лапласа, решающая проблему сетки резисторов

Матрица Лапласа, решающая проблему сетки резисторов

#r #igraph #physics #laplacian

Вопрос:

Я пытаюсь решить для потоков в сети, используя матрицу лапласа. Я начал с тестирования этой проблемы здесь: https://rosettacode.org/wiki/Resistor_mesh#Python

И это идеально подходит для решения, когда все веса равны 1, R = 1,6089 независимо от того, что. Я хочу иметь возможность решать для резисторов, которые не равны одному, хотя! затем можно будет получить токи на каждом резисторе, поэтому я попытался создать набор случайных значений для проводимости (веса) и проверить, работает ли он, чтобы увидеть, что сумма токов, поступающих из узла, в который вводится ток, равна току, вводимому, как и должно быть. К сожалению, это было не так. Я осмотрел лапласиан, и он выглядит примерно так, как я ожидал, но помимо этого я совершенно потерян, может ли кто-нибудь пролить свет на то, есть ли у меня здесь правильный лапласиан? или если я упускаю что-то действительно очевидное, например, неверную индексацию?

Код приведен ниже (извините за очень любительский способ, которым я создал сетку, любые советы и рекомендации по этому вопросу тоже приветствуются!):

 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

Ответ №1:

Ответ: функция 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)

давая токи, которые суммируются в один:

 > 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
Метки: ,