克希荷夫定律與矩陣運算

作者：王一哲
日期：2024年7月8日

原理

當電路上只有電池與電阻器，但又不是單純的串聯、並聯電路時，可以利用克希荷夫定律 (Kirchhoff's circuit laws) 求所有的分支電流。克希荷夫定律有兩個定則

1. 節點定則：對於電路上的某個節點而言，流入、流出的電流量值相等，實際上就是電荷守恆。
2. 迴路定則：對於電路上的某個迴路而言，繞一圈回到出發點時電位相等，實際上就是能量守恆。

下圖是 WikiPedia 克希荷夫定律英文版條目的例子

圖片來源：https://upload.wikimedia.org/wikipedia/commons/7/77/Kirshhoff-example.svg

由圖中右方的節點可得 $$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}$$

矩陣運算

假設要解的未知數有 $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}$$

使用 Python 的 NumPy 套件求解

NumPy 有現成的函式 numpy.linalg.solve 可以求解，也可以先用 numpy.linalg.inv 求反矩陣 $A^{-1}$，再用矩陣乘法計算 $A^{-1} B = x$。

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)  # 印出理論值

輸出為

# 5種寫法的輸出皆為
[[ 0.00090909]
 [ 0.01454545]
 [-0.01363636]]
# 理論值
0.0009090909090909091 0.014545454545454545 -0.013636363636363636

使用 Python 但不使用 NumPy 套件求解

以下的自訂函式可以參考前一篇文章〈使用 Python 及 C++ 計算矩陣乘法、反矩陣〉，可以得到一樣的答案。

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)

輸出為

[[0.000909090909090908], [0.014545454545454545], [-0.013636363636363637]]

使用 C++ 求解

以下的自訂函式可以參考前一篇文章〈使用 Python 及 C++ 計算矩陣乘法、反矩陣〉，可以得到一樣的答案。

#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;
}

輸出為

0.000909091
0.0145455
-0.0136364

執行時間

我在 Ubuntu 指令介面中使用 time 測試這幾個寫法的執行時間，很意外地發現用 Python list 的寫法執行速度居然有點快，C++ 的部分是測試編輯後的可執行檔，速度明顯快很多。

程式碼	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.422	1.282
先用 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

參考資料

1. WikiPedia Kirchhoff's circuit laws
2. numpy.linalg.solve
3. numpy.linalg.inv

---

HackMD 版本連結：https://hackmd.io/@yizhewang/HkhhK8HDR