这是上个学期没有全部完成的话题。借此来复习线性代数和科学计算。
建立模型
设在电阻网络中：
- 有着 $n$ 个节点,编号为 $1\sim n$ ;
- 在不同节点间,有着 $m$ 条接线,编号为 $1\sim m$ ;
- 设第$k$条接线是从节点 $from(k)$ 连接到节点 $to(k)$ 的,且 $from(k)<to(k)$ ,规定从节点 $from(k)$ 到节点 $to(k)$ 为接线 $k$ 的正方向,两个节点之间,至多存在一条接线;
- 设节点 $k$ 上的电压为 $u_k$ ,记 $U = [u_1,…,u_n]^T$ ;
- 设节点 $k$ 上的注入电流为 $i_k$ ,若为流出则 $i_k<0$ ,记 $I = [i_1,…,i_n]^T$ ;
- 设接线 $k$ 上的电阻为 $r_k$ ;
- 设接线 $k$ 上的电压差为 $u_{w_k}=u_{from(k)}-u_{to(k)}$ ,记 $u_w = [u_{w_1},…,u_{w_m}]^T$ ;
- 设接线 $k$ 上的电流为 $i_{w_k}$ ,以接线 $k$ 的正方向为电流正方向,记 $i_w = [i_{w_1},…,i_{w_m}]^T$ .
考虑连通电阻网络,令 $m \ge n -1$。
建立方程
首先考虑$u_w$和$u$的关系,由定义
$$\begin{cases}u_{w_{1}} = u_{from(1)} - u_{to(1)}\\…\\u_{w_{m}} = u_{from(m)} - u_{to(m)}\\\end{cases}\tag{1}$$
定义$m\times n$矩阵$A$，其中
$$\begin{cases}A_{ij} = 1, \;\;\quad\quad\text{若}j = from(i);\\ A_{ij} = -1, \;\;\;\quad\text{若}j = to(i);\\ A_{ij} = 0,\;\;\quad\quad\text{其他} .\end{cases}\tag{2}$$
那么有
$$u_w = A U\tag{3}$$
再考虑$u_w$和$i_w$的关系。欧姆定律得
$$\begin{cases}i_{w_{1}} = \frac{1}{r_1}u_{w_{1}}\\…\\i_{w_{m}} = \frac{1}{r_m}u_{w_{m}}\\\end{cases}\tag{4}$$
定义$m\times m$矩阵$C$，其中
$$\begin{cases}C_{ij} = \frac{1}{r_i}, \;\;\quad\quad\text{若}i = j;\\C_{ij} = 0,\;\;\;\quad\quad\text{若}i \ne j.\\\end{cases}\tag{5}$$
那么有
$$i_w = C u_w\tag{6}$$
最后考虑$i$和$i_w$的关系，由电流守恒方程
$$\begin{cases}i_1 = \underset{from(k) = 1} {\sum i_{w_k}} - \underset{to(k) = 1} {\sum i_{w_k}}\\…\\i_n = \underset{from(k) = n} {\sum i_{w_k}} - \underset{to(k) = n} {\sum i_{w_k}}\\\end{cases}\tag{7}$$
故有
$$I = A^T i_w\tag{8}$$
结合(3)(6)(8)式，得到
$$I = A^TC AU\tag{9}$$
这样就建立起各个节点注入电流与电压之间的关系。
拉普拉斯矩阵[数学准备]
$A^TCA$实际上为拉普拉斯矩阵。考虑其具体性质。
拉普拉斯矩阵的表示
由前，知$A_{m\times n}$:
$$\begin{cases}A_{ij} = 1, \;\;\quad\quad\text{若}j = from(i);\\A_{ij} = -1, \;\;\;\quad\text{若}j = to(i);\\A_{ij} = 0,\;\;\quad\quad\text{其他} .\end{cases}$$
由前，知$C_{m\times m}$:
$$\begin{cases}C_{ij} = \frac{1}{r_i}, \;\;\quad\quad\text{若}i = j;\\C_{ij} = 0,\;\;\;\quad\quad\text{若}i \ne j.\\\end{cases}$$
故记$A^TCA = L_{n\times n}$，有
$$L_{ij} = \displaystyle\sum_{k=1}^m C_{kk}A_{ki}A_{kj}=\begin{cases}\underset{\text{接线}k\text{的一端为}i}{\sum \frac{1}{r_k}} , \qquad\text{若}i=j;\\\underset{k\text{为}i,j\text{间的接线}}{\sum-\frac{1}{r_k}} ,\qquad \text{若}i\ne j.\end{cases}\tag{10}$$
显然，有
$$\sum_{j=1}^n L_{ij} = 0\tag{11}$$
拉普拉斯矩阵的秩
下面考虑$L$的秩。令$LU = I, I = \vec{0}$，解这个方程。
设$u_t=max\{u_1,…,u_n\}$。考虑 $i_t=(LU)_{t}$。有
$$\begin{align*}0 &= i_t = \sum_{j=1}^n L_{tj}u_j\\&=\underset{i\le j \le n,j\ne t}{\sum L_{tj}u_j} + L_{tt}u_t\\&\ge u_t\left(L_{tt} + \underset{i\le j \le n,j\ne t}{\sum L_{tj} }\right)\qquad [L_{tj} \le 0, t\ne j]\\&=u_t \sum_{j=1}^n L_{tj} = 0\tag{12}\end{align*}$$
故不等号取等，知
$$u_r = u_t,\text{若t,r间有接线}\tag{13}$$
而全图是连通的，因此扩散可得
$$u_1 = u_2 = … = u_n\tag{14}$$
故$LU = \vec{0}$的解集为
$$U = \theta [1,1,…,1]^T, \theta \in \mathbb{R}\tag{15}$$
这表明($C(L)$为$L$的列空间)
$$dim(C(L)) = rank(L) = n-1\tag{16}$$
这是符合物理解释的。(15)式表示电压可以同时加/减同一数值而不改变电流分布。(16)式表示电压的选取具有一个自由度，即(15)式所述。
$LU=I$的解
而对于线性方程组$LU = I$，若使得方程有解，必须有$ I \in C(L)$。由(16)式，知$C(L)$具有一个线性约束。而由(15)式，知
$$[1,1,…,1] (LU) = (L[1,1,…,1]^T)^TU = \vec{0}\tag{17}$$
因此，若使方程有解，需有
$$[1,1,…,1]^T I = 0\tag{18}$$
这是符合物理意义的。实际上，这反映了电流守恒。在这种情况下，方程的解
$$U = U^* + \theta [1,1,…,1]^T,\theta \in \mathbb{R},U^*\text{是}LU = I \text{的特解}\tag{19}$$
拉普拉斯矩阵的余子式
下面考虑拉普拉斯矩阵$L$的余子式。记$L$中去掉第$i$行，第$j$列的$(n-1)\times(n-1)$矩阵为$M_{ij}$。那么对于$M_{i_0i_0}(i_0 \in [n])$，结合(10)(11)式，有
$$\begin{align*}&\underset{1\le j \le n-1,j\ne i}{\sum |(M_{i_0i_0})_{ij}|}\qquad\qquad [1\le i \le n-1]\\ = &\underset{1\le j \le n,j\ne i^\prime,i_0}{\sum |L_{i^\prime j}|}\qquad\qquad\quad [\text{若}i< i_0,i^\prime = i; \text{否则} i^\prime = i + 1] \\=&\underset{1\le j \le n,j\ne i^\prime,i_0}{-\sum L_{i^\prime j}}\qquad\qquad\quad[L_{i^\prime j} \le 0, i^\prime \ne j]\\=&L_{i^\prime i^\prime} + L_{i^\prime i_0}\qquad\\\le &|L_{i^\prime i^\prime}| = |(M_{i_0i_0})_{ii}|\tag{20}\end{align*}$$
同时，取等要求$L_{i^\prime i_0} = 0$，这表明节点$i^\prime$和$i_0$之间无接线。对于所有$i^\prime$，若均取等，表明节点$i_0$和所有节点均无接线，这和全图连通矛盾。因此必然存在某个$i^\prime$使得不等式无法取等。因此
$$M_{i_0i_0}\text{是弱对角占优矩阵}\tag{21}$$
同时，由于全图是连通的，且$L_{ij}\ne 0(i\ne j)$当且仅当$i,j$间有接线时成立，因此易得$M_{i_0i_0}$对应的有向图是强连通的。因此可得
$$M_{i_0i_0}\text{是不可约矩阵}\tag{22}$$
由于弱对角占优不可约矩阵是非奇异的，故结合(21)(22)式得
$$det(M_{i_0i_0}) \ne 0\tag{23}$$
求解电阻网络的相关量
求两节点间的等效电阻
由(9)式知$LU = A^TCAU = I$。现在要计算节点$i_0,i_1$间的等效电阻(为了方便，约定$i_0<i_1$)。
由定义，应令
$$i_t=\begin{cases}i_{in},\qquad\;\;\; t= i_0;\\-i_{in},\qquad t=i_1;\\0,\qquad\quad\; t\ne i_0,i_1.\end{cases}\tag{24}$$
对此$I$，解出$U$。由(19)式，注意到解不是唯一的。事实上，关注的是电势差，因此不妨固定$u_{i_1}=0$，此时相当于寻求特解$U^*$。
当$u_{i_1}=0$时，$L$的第$i_1$列可以忽略(与$u_{i_1}=0$相乘)，同时忽略$L$的第$i_1$行，构建新方程。记$I$删去第$i_1$个分量后剩下的列向量为$\hat{I}$，记$U$删去第$i_1$个分量后剩下的列向量为$\hat{U}$，那么
$$M_{i_1i_1}\hat{U} = \hat{I}\tag{25}$$
由(23)式知$det(M_{i_1i_1})\ne 0$。因此(25)式有且仅有唯一解。求解等效电阻，目标是求出$u_{i_0}$。克莱姆法则，记$M_{i_1i_1}$中第$i_0$列置换成$\hat{I}$后的矩阵为$\hat{M}$，那么可得
$$u_{i_0} = \frac{det(\hat{M})}{det(M_{i_1i_1})}\tag{26}$$
将$det(\hat{M})$关于$\hat{M}$的第$i_0$列展开，由于$\hat{I}$仅有第$i_0$个分量不为0，因此记$L$去掉第$i_0,i_1$行和列后的矩阵为$M_{i_0i_0i_1i_1}$，有
$$det(\hat{M}) = (-1)^{i_0+i_0}\cdot i_{in}\cdot det(M_{i_0i_0i_1i_1})\tag{27}$$
因此，结合以上两式，得
$$u_{i_0} = i_{in}\frac{det(M_{i_0i_0i_1i_1})}{det(M_{i_1i_1})}\tag{28}$$
【事实上，注意到这里$i_0$和$i_1$是完全对称的，因此(28)式分母里的$det(M_{i_1i_1})$也可以换成$det(M_{i_0i_0})$。这说明了
$$det(M_{i_0i_0})=det(M_{i_1i_1})\tag{29}$$
由于$i_0,i_1$的选取是任意的，因此对任意的$i$，$det(M_{ii})$均相等。】
由(28)式，知$i_0,i_1$间等效电阻
$$R_{i_0i_1}=\frac{u_{i_0}-u_{i_1}}{i_{in}} =\frac{det(M_{i_0i_0i_1i_1})}{det(M_{i_1i_1})}\tag{30}$$
对于特定的两点间的等效电阻，利用高斯消元可以在线性时间内完成行列式的计算，复杂度是可以接受的。
给定节点电流求解节点电压分布
相当于求解方程组$LU=I$。由(18)式，应有 $[1,1,…,1]^T I = 0$，否则方程无解。在此条件下，需要先寻求特解。不妨$u_n=0$。
此时，与(25)式类似，可以忽略$L$的第$n$行和第$n$列。记$I$删去第$n$个分量后剩下的列向量为$\hat{I}$，记$U$删去第$n$个分量后剩下的列向量为$\hat{U}$，那么有
$$M_{nn}\hat{U}=\hat{I}\tag{25*}$$
由(23)式，$det(M_{nn})\ne 0$，方程有且仅有唯一解。这里若采用克莱姆法则求解，计算复杂度过大，无法接受。因此采用$Gauss-Seidel$迭代法求解。记
$$S_{(n-1)\times (n-1)}=\begin{bmatrix}L_{11}&0&0&…&0\\L_{21}&L_{22}&0&…&0\\L_{31}&L_{32}&L_{33}&…&0\\\vdots&\vdots&\vdots&\ddots&\vdots\\L_{n-1,1}&L_{n-1,2}&L_{n-1,3}&…&L_{n-1,n-1}\\\end{bmatrix}\tag{31.1}$$
$$R_{(n-1)\times (n-1)}=\begin{bmatrix}0&L_{12}&L_{13}&…&L_{1,n-1}\\0&0&L_{23}&…&L_{2,n-1}\\0&0&0&…&L_{3,n-1}\\\vdots&\vdots&\vdots&\ddots&\vdots\\0&0&0&…&0\\\end{bmatrix}\tag{31.2}$$
那么$Gauss-Seidel$迭代法，有
$$S\hat{U}^{(k+1)} = R\hat{U}^{(k)} + \hat{I}\tag{32}$$
每次迭代后，由于$S$是三角矩阵，因此可以快速解出$\hat{U}^{(k+1)}$。又因为(21)(22)证明了$M_{nn}$是弱对角占优不可约矩阵，因此
$$Gauss-Seidel\text{迭代法收敛}\tag{33}$$
据此，可以通过迭代较快地解出$\hat{U}$，得到特解$U^*$。据(19)式，可以给出通解(物理意义上，没有改变电势差)。