<style>
.markdown-body table{
display: unset;
}
</style>
# 克希荷夫定律與矩陣運算
> 作者：王一哲
> 日期：2024年7月8日
<br />
## 原理
當電路上只有電池與電阻器，但又不是單純的串聯、並聯電路時，可以利用**克希荷夫定律 (Kirchhoff's circuit laws)** 求所有的分支電流。克希荷夫定律有兩個定則
1. **節點定則**：對於電路上的某個節點而言，流入、流出的電流量值相等，實際上就是**電荷守恆**。
2. **迴路定則**：對於電路上的某個迴路而言，繞一圈回到出發點時電位相等，實際上就是**能量守恆**。
下圖是 [WikiPedia 克希荷夫定律英文版條目](https://en.wikipedia.org/wiki/Kirchhoff%27s_circuit_laws)的例子
<img height="40%" width="40%" src="https://upload.wikimedia.org/wikipedia/commons/7/77/Kirshhoff-example.svg" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">圖片來源：<a href="https://upload.wikimedia.org/wikipedia/commons/7/77/Kirshhoff-example.svg" target="_blank">https://upload.wikimedia.org/wikipedia/commons/7/77/Kirshhoff-example.svg</a></div>
<br /><br />
由圖中右方的節點可得
$$
i_1 = i_2 + i_3
$$
由上方的迴路可得
$$
\varepsilon_1 - i_1 R_1 - i_2 R_2 = 0
$$
由下方的迴路可得
$$
-i_3 R_3 - \varepsilon_2 - \varepsilon_1 + i_2 R_2 = 0
$$
整理以上3式可得
$$
\begin{cases}
i_1 - i_2 - i_3 = 0 \\
i_1 R_1 - i_2 R_2 = \varepsilon_1 \\
i_2 R_2 - i_3 R_3 = \varepsilon_1 + \varepsilon_2
\end{cases} ~\Rightarrow~
\begin{cases}
1 \cdot i_1 + (-1) \cdot i_2 + (-1) \cdot i_3 = 0 \\
R_1 \cdot i_1 + (-R_2) \cdot i_2 + 0 \cdot i_3 = \varepsilon_1 \\
0 \cdot i_1 + R_2 \cdot i_2 + (-R_3) \cdot i_3 = \varepsilon_1 + \varepsilon_2
\end{cases}
$$
如果轉換成矩陣型式
$$
\begin{bmatrix}
1 & -1 & -1 \\
R_1 & R_2 & 0 \\
0 & R_2 & -R_3
\end{bmatrix}
\begin{bmatrix}
i_1 \\ i_2 \\ i_3
\end{bmatrix} = \begin{bmatrix}
0 \\ \varepsilon_1 \\ \varepsilon_1 + \varepsilon_2
\end{bmatrix}
$$
假設代入的數值分別為
$$
R_1 = 100 ~\mathrm{\Omega},~ R_2 = 200 ~\mathrm{\Omega},~ R_3 = 300 ~\mathrm{\Omega},~ \varepsilon_1 = 3 ~\mathrm{V},~ \varepsilon_2 = 4 ~\mathrm{V}
$$
電流分別為
$$
i_1 = \frac{1}{1100} ~\mathrm{A},~ i_2 = \frac{4}{275}~\mathrm{A},~ i_3 = -\frac{3}{220} ~\mathrm{A}
$$
<br /><br />
## 矩陣運算
假設要解的未知數有 $n$ 個，維度為 $n \times n$ 的矩陣 $A$ 是聯立方程式中每一個未知數的係數，維度為 $n \times 1$ 的矩陣 $x$ 是未知數，維度為 $n \times 1$ 的矩陣 $B$ 是每一列方程式對應的常數項。第一種計算方法是利用 $A$ 的反矩陣 $A^{-1}$，計算過程如下：
$$
Ax = B ~\Rightarrow~ A^{-1}Ax = A^{-1}B ~\Rightarrow~ x = A^{-1} B
$$
以上面的電路圖為例，
$$
A =
\begin{bmatrix}
1 & -1 & -1 \\
100 & 200 & 0 \\
0 & 200 & -300
\end{bmatrix} ~~~~~ A^{-1} =
\begin{bmatrix}
\frac{6}{11} & \frac{1}{220} & -\frac{1}{550} \\
-\frac{3}{11} & \frac{3}{1100} & \frac{1}{1100} \\
-\frac{2}{11} & \frac{1}{550} & -\frac{3}{1100}
\end{bmatrix}
$$
高中數學課本上比較常見的反矩陣求法，是將矩陣 $A$ 與對應的單位矩陣 $I$ 左、右並列，再用高斯消去法，將左側的矩陣化為單位矩陣，此時右側就是反矩陣 $A^{-1}$，計算過程如下
$$
\begin{bmatrix}
1 & -1 & -1 & | & 1 & 0 & 0 \\
100 & 200 & 0 & | & 0 & 1 & 0 \\
0 & 200 & -300 & | & 0 & 0 & 1
\end{bmatrix}
$$
先將第2、3列除以100可得
$$
\begin{bmatrix}
1 & -1 & -1 & | & 1 & 0 & 0 \\
1 & 2 & 0 & | & 0 & \frac{1}{100} & 0 \\
0 & 2 & -3 & | & 0 & 0 & \frac{1}{100}
\end{bmatrix}
$$
將第1列乘以 $-1$ 加到第2列可得
$$
\begin{bmatrix}
1 & -1 & -1 & | & 1 & 0 & 0 \\
0 & 3 & 1 & | & －1 & \frac{1}{100} & 0 \\
0 & 2 & -3 & | & 0 & 0 & \frac{1}{100}
\end{bmatrix}
$$
將第2列乘以 $-\frac{2}{3}$ 加到第3列可得
$$
\begin{bmatrix}
1 & -1 & -1 & | & 1 & 0 & 0 \\
0 & 3 & 1 & | & -1 & \frac{1}{100} & 0 \\
0 & 0 & -\frac{11}{3} & | & \frac{2}{3} & -\frac{1}{150} & \frac{1}{100}
\end{bmatrix}
$$
將第2列乘以 $\frac{1}{3}$、第3列乘以 $-\frac{3}{11}$ 可得
$$
\begin{bmatrix}
1 & -1 & -1 & | & 1 & 0 & 0 \\
0 & 1 & \frac{1}{3} & | & -\frac{1}{3} & \frac{1}{300} & 0 \\
0 & 0 & 1 & | & -\frac{2}{11} & \frac{1}{550} & -\frac{3}{1100}
\end{bmatrix}
$$
將第3列乘以 $-\frac{1}{3}$ 加到第2列、第3列加到第1列可得
$$
\begin{bmatrix}
1 & -1 & 0 & | & \frac{9}{11} & \frac{1}{550} & -\frac{3}{1100} \\
0 & 1 & 0 & | & -\frac{3}{11} & \frac{3}{1100} & \frac{1}{1100} \\
0 & 0 & 1 & | & -\frac{2}{11} & \frac{1}{550} & -\frac{3}{1100}
\end{bmatrix}
$$
將第2列加到第1列可得
$$
\begin{bmatrix}
1 & 0 & 0 & | & \frac{6}{11} & \frac{1}{220} & -\frac{1}{550} \\
0 & 1 & 0 & | & -\frac{3}{11} & \frac{3}{1100} & \frac{1}{1100} \\
0 & 0 & 1 & | & -\frac{2}{11} & \frac{1}{550} & -\frac{3}{1100}
\end{bmatrix}
$$
接下來再將 $A^{-1}$ 乘以 $B$ 即可求解。
$$A^{-1}B =
\begin{bmatrix}
\frac{6}{11} & \frac{1}{220} & -\frac{1}{550} \\
-\frac{3}{11} & \frac{3}{1100} & \frac{1}{1100} \\
-\frac{2}{11} & \frac{1}{550} & -\frac{3}{1100}
\end{bmatrix} \begin{bmatrix}
0 \\ 3 \\ 7
\end{bmatrix} = \begin{bmatrix}
\frac{1}{1100} \\ \frac{4}{275} \\ -\frac{3}{220}
\end{bmatrix}
$$
高中數學課本上比較常見的寫法是將要解的矩陣 $A$ 與每一列對應的常數項 $B$ 左右並列，再用高斯消去法，將左側的矩陣化為單位矩陣，此時右側的矩陣就是對應的解。計算過程如下：
$$
\begin{bmatrix}
1 & -1 & -1 & | & 0 \\
100 & 200 & 0 & | & 3 \\
0 & 200 & -300 & | & 7
\end{bmatrix}
$$
先將第2、3列除以100可得
$$
\begin{bmatrix}
1 & -1 & -1 & | & 0 \\
1 & 2 & 0 & | & \frac{3}{100} \\
0 & 2 & -3 & | & \frac{7}{100}
\end{bmatrix}
$$
將第1列乘以 $-1$ 加到第2列可得
$$
\begin{bmatrix}
1 & -1 & -1 & | & 0 \\
0 & 3 & 1 & | & \frac{3}{100} \\
0 & 2 & -3 & | & \frac{7}{100}
\end{bmatrix}
$$
將第2列乘以 $-\frac{2}{3}$ 加到第3列可得
$$
\begin{bmatrix}
1 & -1 & -1 & | & 0 \\
0 & 3 & 1 & | & \frac{3}{100} \\
0 & 0 & -\frac{11}{3} & | & \frac{1}{20}
\end{bmatrix}
$$
將第3列乘以 $-\frac{3}{11}$ 加到第3列可得
$$
\begin{bmatrix}
1 & -1 & -1 & | & 0 \\
0 & 3 & 1 & | & \frac{3}{100} \\
0 & 0 & 1 & | & -\frac{3}{220}
\end{bmatrix}
$$
將第3列加到第1列，第3列乘以 $-1$ 加到第2列可得
$$
\begin{bmatrix}
1 & -1 & 0 & | & -\frac{3}{220} \\
0 & 3 & 0 & | & \frac{12}{275} \\
0 & 0 & 1 & | & -\frac{3}{220}
\end{bmatrix}
$$
將第2列除以3可得
$$
\begin{bmatrix}
1 & -1 & 0 & | & -\frac{3}{220} \\
0 & 1 & 0 & | & \frac{4}{275} \\
0 & 0 & 1 & | & -\frac{3}{220}
\end{bmatrix}
$$
最後將第2列加到第1列可得
$$
\begin{bmatrix}
1 & 0 & 0 & | & \frac{1}{1100} \\
0 & 1 & 0 & | & \frac{4}{275} \\
0 & 0 & 1 & | & -\frac{3}{220}
\end{bmatrix}
$$
<br /><br />
## 使用 Python 的 NumPy 套件求解
NumPy 有現成的函式 numpy.linalg.solve 可以求解，也可以先用 numpy.linalg.inv 求反矩陣 $A^{-1}$，再用矩陣乘法計算 $A^{-1} B = x$。
```python=
import numpy as np
A = np.asarray([[1, -1, -1],
[100, 200, 0],
[0, 200, -300]])
B = np.asarray([[0], [3], [7]])
print(np.linalg.solve(A, B)) # 使用 numpy.linalg.solve 求解
AI = np.linalg.inv(A) # 使用 numpy.linalg.inv 求 A 的反矩陣 AI
print(AI @ B) # 再用矩陣乘法求解
print(AI.dot(B))
print(np.dot(AI, B))
print(np.matmul(AI, B))
print(1/1100, 4/275, -3/220) # 印出理論值
```
輸出為
```python
# 5種寫法的輸出皆為
[[ 0.00090909]
[ 0.01454545]
[-0.01363636]]
# 理論值
0.0009090909090909091 0.014545454545454545 -0.013636363636363636
```
<br /><br />
## 使用 Python 但不使用 NumPy 套件求解
以下的自訂函式可以參考前一篇文章〈[使用 Python 及 C++ 計算矩陣乘法、反矩陣](https://hackmd.io/@yizhewang/S1IUjnUDC)〉，可以得到一樣的答案。
```python=
import copy
def matmul(A, B):
m, n, p = len(A), len(A[0]), len(B[0])
C = [[0]*p for _ in range(m)]
for i in range(m):
for j in range(p):
for k in range(n):
C[i][j] += A[i][k]*B[k][j]
return C
def inv(A):
# 取 A 的列、欄數量，如果不是方陣，回傳 ValueError 並中止程式
n, m = len(A), len(A[0])
if n != m: raise ValueError
# 建立右側的單位矩陣 B
B = [[0]*n for _ in range(n)]
for i in range(n): B[i][i] = 1.0
# 使用 copy.deepcopy 複製 A
test = copy.deepcopy(A)
# 消除 test 左下角的數值
for i in range(n): # 依序處理第 0 ~ (n-1) 列
coef = test[i][i] # 以對角線上的元素作為係數
for j in range(i+1, n): # 處理第 (i+1) ~ (n-1) 列
mul = -test[j][i] / coef # 計算要將 test[i] 同乘的倍數 mul
for k in range(n): # 處理第 0 ~ (n-1) 欄
test[j][k] += mul*test[i][k] # 化簡 test 的左下角
B[j][k] += mul*B[i][k] # 同時處理 B
# 消除 test 右下角的數值
for i in range(n-1, -1, -1): # 依序處理第 (n-1) ~ 0 列
coef = test[i][i] # 以對角線上的元素作為係數
for j in range(i-1, -1, -1): # 處理第 (i-1) ~ (0) 列
mul = -test[j][i] / coef # 計算要將 test[i] 同乘的倍數 mul
for k in range(n-1, -1, -1): # 處理第 (n-1) ~ 0 欄
test[j][k] += mul*test[i][k] # 化簡 test 的右下角
B[j][k] += mul*B[i][k] # 同時處理 B
# 將 test 對角線上的元素都變為1
for i in range(n): # 依序處理第 0 ~ (n-1) 列
coef = test[i][i] # 以對角線上的元素作為係數
test[i][i] /= coef # 處理對角線上的元素
for j in range(n): # 依序處理 B 的第 0 ~ (n-1) 欄
B[i][j] /= coef
# 回傳反矩陣 B
return B
if __name__ == "__main__":
A = [[1, -1, -1],
[100, 200, 0],
[0, 200, -300]]
B = [[0], [3], [7]]
AI = inv(A)
x = matmul(AI, B)
print(x)
```
輸出為
```python
[[0.000909090909090908], [0.014545454545454545], [-0.013636363636363637]]
```
<br /><br />
## 使用 C++ 求解
以下的自訂函式可以參考前一篇文章〈[使用 Python 及 C++ 計算矩陣乘法、反矩陣](https://hackmd.io/@yizhewang/S1IUjnUDC)〉，可以得到一樣的答案。
```cpp=
#include <iostream>
#include <vector>
using namespace std;
void inv(vector<vector<float>> &A, vector<vector<float>> &AI) {
int n = (int)A.size(); // A 的列、欄數量
for(int i=0; i<n; i++) AI[i][i] = 1.0; // 將 AI 的對角線元素都設定為 1.0
vector<vector<float>> test (A.begin(), A.end()); // 複製 A 的資料到 test
/* 消除 test 左下角的數值 */
for(int i=0; i<n; i++) { // 依序處理第 0 ~ (n-1) 列
float coef = test[i][i]; // 以對角線上的元素作為係數
for(int j=i+1; j<n; j++) { // 處理第 (i+1) ~ (n-1) 列
float mul = -test[j][i] / coef; // 計算要將 test[i] 同乘的倍數 mul
for(int k=0; k<n; k++) { // 處理第 0 ~ (n-1) 欄
test[j][k] += mul*test[i][k]; // 化簡 test 的左下角
AI[j][k] += mul*AI[i][k]; // 同時處理 AI
}
}
}
/* 消除 test 右上角的數值 */
for(int i=n-1; i>=0; i--) { // 依序處理第 (n-1) ~ 0 列
float coef = test[i][i]; // 以對角線上的元素作為係數
for(int j=i-1; j>=0; j--) { // 處理第 (i-1) ~ (0) 列
float mul = -test[j][i] / coef; // 計算要將 A[i] 同乘的倍數 mul
for(int k=n-1; k>=0; k--) { // 處理第 (n-1) ~ 0 欄
test[j][k] += mul*test[i][k]; // 化簡 test 的右下角
AI[j][k] += mul*AI[i][k]; // 同時處理 AI
}
}
}
/* 將 test 對角線上的元素都變為1 */
for(int i=0; i<n; i++) { // 依序處理第 0 ~ (n-1) 列
float coef = test[i][i]; // 以對角線上的元素作為係數
test[i][i] /= coef; // 處理對角線上的元素
for(int j=0; j<n; j++) AI[i][j] /= coef; // 依序處理 AI 的第 0 ~ (n-1) 欄
}
}
void matmul(vector<vector<float>> &A, vector<vector<float>> &B, vector<vector<float>> &x) {
size_t m = A.size(), n = A[0].size(), p = B[0].size();
for(size_t i=0; i<m; i++)
for(size_t j=0; j<p; j++)
for(size_t k=0; k<n; k++)
x[i][j] += A[i][k]*B[k][j];
}
int main() {
vector<vector<float>> A = {{1.0, -1.0, -1.0},
{100.0, 200.0, 0.0},
{0.0, 200.0, -300.0}};
vector<vector<float>> B = {{0.0}, {3.0}, {7.0}};
vector<vector<float>> AI (A.size(), vector<float> (A[0].size(), 0.0));
vector<vector<float>> x (B.size(), vector<float> (B[0].size(), 0.0));
inv(A, AI);
matmul(AI, B, x);
for(size_t i=0; i<x.size(); i++)
for(size_t j=0; j<x[0].size(); j++)
cout << x[i][j] << " \n"[j == x[0].size()-1];
return 0;
}
```
輸出為
```cpp
0.000909091
0.0145455
-0.0136364
```
<br /><br />
### 執行時間
我在 Ubuntu 指令介面中使用 time 測試這幾個寫法的執行時間，很意外地發現用 Python list 的寫法執行速度居然有點快，C++ 的部分是測試編輯後的可執行檔，速度明顯快很多。
<div style="text-align:center">
<table>
| <center>程式碼</center> | real (s) | user (s) | sys (s) |
| ----- | -------- | -------- | ------- |
| np.linalg.solve(A, B) | 0.153 | 0.432 | 1.230 |
| 先用 np.linalg.inv 再用 AI @ B | 0.153 | 0.426 | 1.246 |
| 先用 np.linalg.inv 再用 AI.dot(B) | 0.128 | 0.422 | 1.282 |
| 先用 np.linalg.inv 再用 np.dot(AI, B) | 0.157 | 0.455 | 1.304 |
| 先用 np.linalg.inv 再用 np.matmul(AI, B) | 0.132 | 0.423 | 1.264 |
| Python list | 0.014 | 0.014 | 0.000 |
| C++ | 0.002 | 0.002 | 0.000 |
</table>
</div>
<br /><br />
## 參考資料
1. [WikiPedia Kirchhoff's circuit laws](https://en.wikipedia.org/wiki/Kirchhoff%27s_circuit_laws)
2. [numpy.linalg.solve](https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html)
3. [numpy.linalg.inv](https://numpy.org/doc/stable/reference/generated/numpy.linalg.inv.html)
---
###### tags:`Physics`、`Python`、`C++`、`Math`
Sign in with Wallet
Connect another wallet