引入

我们早就学过了计算串并联电阻的等效电阻的公式:

针对混合的电阻网络,有些时候我们依然可以反复运用这两个公式来计算等效电阻。在有些电阻网络中,这样是无法得到答案的:

如果你熟悉星三角变换,你应该能很快算出答案。(或者你能看出这就是一个平衡桥,那就能口算了)

星三角变换固然是一种手工计算等效电阻的很好方法,但是对于计算机而言,利用这个方法计算电阻会很麻烦。下面,我将使用高斯消元法来计算两点之间的等效电阻。

建模

我们可以把两点之间的纯电阻网络看作是一张有向图,节点表示导线,边表示电阻,边权表示电阻的大小。要求源点和汇点之间的等效电阻大小。

为什么要建成有向边呢,因为虽然电阻是无向的,但是电流是有向的,我们一开始并不知道电流的方向,但是可以任意地假定一个方向,然后给每一条边赋上一个电流IiI_i,当IiI_i为正,表示电流方向与箭头同向,当IiI_i为负,表示电流方向与箭头反向。然后我们再假定汇点的电势为00,其余每个节点赋上一个电势φi\varphi_i

若要求源点到汇点的等效电阻,只需要知道源点的电势φs\varphi_s和汇点的电势φt\varphi_t,以及汇入汇点的总电流II,就可以计算出等效电阻R=φsφtIR = \frac{\varphi_s - \varphi_t}{I}。为了计算这些量,我们要把整个网络的电流和电势都算出来。

设这张图共有nn节点mm条边,则我们一共有mmIiI_i未知,n2n-2φi\varphi_i未知。一共m+n2m+n-2个未知量,需要m+n2m+n-2个方程才能解出来。我们可以利用所学的物理知识来找到这些方程。

欧姆定律

最简单的就是安培定律了,对于编号为ii的边(u,v)(u, v),有

φuφv=IiRi\begin{align*} \varphi_u - \varphi_v = I_i R_i\\ \end{align*}

这一共有mm个方程。

基尔霍夫电流定律(KCL)

除了源点和汇点外,每个结点流入的电流都等于流出的电流。

对于编号为uu的节点,有

iAdj(u)Ii=0\sum_{i \in Adj(u)} I_i = 0

n2n-2 个方程,再加上前面的 mm 个,正好有 m+n2m+n-2 个方程,我们可以用高斯消元法来解出所有的电流。

高斯消元法

现在不妨令 φs=1\varphi_s = 1φt=0\varphi_t = 0,则有

{φuφv=IiRi(u,v)EiAdj(u)±Ii=0uV{s,t}φs=1φt=0\begin{cases} \varphi_u - \varphi_v = I_i R_i && \forall (u, v) \in E \\ \sum_{i \in Adj(u)} \pm I_i = 0 && \forall u \in V - \{s, t\} \\ \varphi_s = 1 \\ \varphi_t = 0 \\ \end{cases}

有了这 m+nm+n 个方程,我们就可以用高斯消元法解出所有的电势和电流,进而求出等效电阻了。

高斯消元法是一种求解线性方程组的算法,在计算机中的实现也很简单。此处复杂度O((m+n)3)O((m+n)^3)

此处的优化空间:从上面可以看出,这些方程一般是很稀疏的,如果改用计算稀疏线性方程组的算法可以达到更高的效率。如LU分解、bfgs、lbfgs

代码实现

使用纯C实现

Github地址:https://github.com/cyrus28214/resistance

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

double **array2d(int rows, int cols);
void gauss(double **a, int n);

int main() {
int n, //节点数量
m, //边的数量
s, //源点
t; //汇点
//图不能有自环,可以有重边,需要联通

scanf("%d%d%d%d", &n, &m, &s, &t);

//为了代码方便,增加一个变量I_总
int k = n + m + 1; //方程的个数

double **a = array2d(k, k+1);//多出的一列是等号右边的常数
//每一列代表的变量分别为:
// I_0 ... I_{m-1} | phi_0 ... phi_{n-1} | I_总
//I_i = a[][i]
//phi_i = a[][m+i]

//各个方程为
//0 ~ m-1 欧姆定律: phi_u - phi_v - I_i R_i = 0
//m ~ m+n-1 KCL: sum I_i = 0 (对非源点汇点)
// phi_s = 1 (对源点,加1V电压)
// phi_t = 0 (对汇点,接地)
//m+n: 计算总电流I_总=汇入t电流

//m行,每行三个整数u, v, r表示u到v有一个阻值为r的电阻
for (int i = 0; i < m; i++) {
int u, v;
double r;
scanf("%d%d%lf", &u, &v, &r);

//欧姆定律
a[i][m+u] = 1;
a[i][m+v] = -1;
a[i][i] = -r;
//第i个方程 phi_u - phim+n- I_i R_i = 0

//KCL 注意电流方向,这里按流出为负,流入为正
if (u != s && u != t) a[m+u][i] = -1;
else if (u == t) a[m+n][i] = -1;
if (v != s && v != t) a[m+v][i] = 1;
else if (v == t) a[m+n][i] = 1;
}

a[m+s][m+s] = 1, a[m+s][k] = 1; //源点电电势为1V
a[m+t][m+t] = 1, a[m+t][k] = 0; //汇点电势为0
a[m+n][m+n] = -1; //计算总电流I_总

gauss(a, k); //求解方程组

double I_total = a[k-1][k];
printf("%6e\n", 1 / I_total);//由于我们设定了整体电压为1,最后的等效电阻就是1/I_总
}

//生成一个rows行cols列的二维数组,全部初始化成0
double **array2d(int rows, int cols) {
double **arr = malloc(rows * sizeof(double *));
for (int i = 0; i < rows; i++) {
arr[i] = calloc(cols, sizeof(double));
}
return arr;
}

void gauss(double **a, int n) {
//高斯-约旦消元法求解线性方程组
for (int i = 0; i < n; i++) { //当前行
//将当前行替换为列主元最大的行
for (int j = i + 1; j < n; j++) {
if (fabs(a[j][i]) > fabs(a[i][i])) {
double *temp = a[i];
a[i] = a[j];
a[j] = temp;
}
}

double c = a[i][i];
//通常的算法此处需要特判c是否为0,但根据物理意义,本问题一定有唯一解,可以不用。
for (int k = 0; k <= n; k++)
a[i][k] /= c;

//消元
for (int j = 0; j < n; j++) {
if (i == j) continue;
c = a[j][i];
for (int k = 0; k <= n; k++)
a[j][k] -= c * a[i][k];
}
}
}

验证

对于开头的这个例子,我们可以验证一下:

输入

1
2
3
4
5
6
7
4 5 0 1
0 2 20
0 3 20
1 2 60
1 3 60
3 2 50
4.000000e+01

输出4.000000e+01,和手算结果一致

再来一个例子,经典高中物理题,每个电阻都是1Ω1\Omega,求对角线电阻:

用我们的程序算出:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
8 12 0 6
0 1 1
1 2 1
2 3 1
3 0 1
4 5 1
5 6 1
6 7 1
7 4 1
0 4 1
1 5 1
2 6 1
3 7 1
8.333333e-01

而这一题的答案是5/6Ω5 / 6 \,\Omega,答案正确。

一些私货

qzqz机房的hxd们可以在cofun上看到我出的这一题,快去做(cofun服务器好像又爆炸了)

如果你不知道cofun是什么,请忽略这段私货。