MATLAB自动生成电阻网格网络：关联矩阵的计算方法

原代码的核心问题 #

1. 边数计算错误：n×n的网格包含水平、垂直两类边，总边数应为 2*n*(n-1)，原代码仅计算了一半。
2. 关联矩阵索引逻辑混乱：原代码的行列映射错误，无法正确对应每条边与节点的连接关系。
3. 网络方程构建错误：原代码的方程形式不符合基尔霍夫定律，导致无法正确求解节点电压。

正确构建关联矩阵的步骤 #

关联矩阵是边数×节点数的矩阵，每条边对应一行：边连接的两个节点中，一个记+1，另一个记-1（方向不影响最终结果）。我们分两类边构建：

1. 垂直边（上下节点连接）

共n*(n-1)条，每条边连接第i行第j列节点与第i-1行第j列节点。

• 行索引范围：1 ~ n*(n-1)
• 对于第k条垂直边：
  • 上方节点：(floor((k-1)/n))*n + mod(k-1,n)+1
  • 下方节点：上方节点 + n
  • 关联矩阵对应行：上方节点列填-1，下方节点列填+1

2. 水平边（左右节点连接）

共n*(n-1)条，每条边连接第i行第j列节点与第i行第j+1列节点。

• 行索引范围：n*(n-1)+1 ~ 2*n*(n-1)
• 对于第t条水平边（t从1开始）：
  • 左侧节点：(floor((t-1)/(n-1)))*n + mod(t-1,n-1)+1
  • 右侧节点：左侧节点 + 1
  • 关联矩阵对应行：左侧节点列填-1，右侧节点列填+1

修正后的完整代码 #

% Network Size
n = 49; % Grid side length (n×n nodes)

% Resistance Value
R = 100; % Resistance per edge

% Choose injection side (1=top, 2=right, 3=bottom, 4=left)
sourceSide = 3;

% Calculate total nodes and edges
numNodes = n^2;
numEdges = 2*n*(n-1); % Total edges: vertical + horizontal

% Initialize incidence matrix (edges × nodes)
incidenceMatrix = zeros(numEdges, numNodes);

% ----------------------
% 1. Build vertical edges (top-bottom connections)
% ----------------------
for k = 1:n*(n-1)
    upperNode = floor((k-1)/n)*n + mod(k-1,n) + 1;
    lowerNode = upperNode + n;
    incidenceMatrix(k, upperNode) = -1;
    incidenceMatrix(k, lowerNode) = 1;
end

% ----------------------
% 2. Build horizontal edges (left-right connections)
% ----------------------
offset = n*(n-1); % Start index of horizontal edges
for t = 1:n*(n-1)
    row = floor((t-1)/(n-1)) + 1;
    col = mod(t-1, n-1) + 1;
    leftNode = (row-1)*n + col;
    rightNode = leftNode + 1;
    incidenceMatrix(offset + t, leftNode) = -1;
    incidenceMatrix(offset + t, rightNode) = 1;
end

% ----------------------
% Set up network equations (Kirchhoff's laws)
% ----------------------
% Define node current vector: 1A injected at source node, -1A at reference (last node)
nodeCurrents = zeros(numNodes, 1);
% Select source node based on side
switch sourceSide
    case 1 % Top side center
        sourceNode = ceil(n/2);
    case 2 % Right side center
        sourceNode = n*ceil(n/2);
    case 3 % Bottom side center
        sourceNode = numNodes - ceil(n/2) + 1;
    case 4 % Left side center
        sourceNode = n*floor(n/2) + 1;
end
nodeCurrents(sourceNode) = 1;
nodeCurrents(end) = -1; % Set last node as reference (0V implicitly)

% Build system matrix: (incidenceMatrix' * incidenceMatrix)/R
systemMatrix = (incidenceMatrix' * incidenceMatrix)/R;
% Remove last row/column to avoid singular matrix (reference node voltage=0)
systemMatrix(end,:) = [];
systemMatrix(:,end) = [];
nodeCurrents(end) = [];

% Solve for node potentials
nodePotentials = systemMatrix \ nodeCurrents;
% Add back reference node voltage (0)
nodePotentials = [nodePotentials; 0];

% ----------------------
% Plot node potentials
% ----------------------
figure;
potGrid = reshape(nodePotentials, n, n);
imshow(potGrid, []);
title('Node Potential Distribution');
colormap(jet);
colorbar;
axis equal tight;

代码关键说明 #

1. 参考节点处理：将最后一个节点设为参考（电压0），移除对应方程避免矩阵奇异。
2. 方程推导：基于基尔霍夫电流定律，结合欧姆定律推导得到线性方程组，确保解的正确性。
3. 索引逻辑：通过行列映射清晰对应每条边与节点的连接，避免原代码的索引混乱问题。

内容的提问来源于stack exchange，提问作者Jan