HackMD
克希荷夫定律與矩陣運算
作者：王一哲
日期：2024年7月8日
原理
當電路上只有電池與電阻器，但又不是單純的串聯、並聯電路時，可以利用克希荷夫定律 (Kirchhoff's circuit laws) 求所有的分支電流。克希荷夫定律有兩個定則
- 節點定則：對於電路上的某個節點而言，流入、流出的電流量值相等，實際上就是電荷守恆。
- 迴路定則：對於電路上的某個迴路而言，繞一圈回到出發點時電位相等，實際上就是能量守恆。
下圖是 WikiPedia 克希荷夫定律英文版條目的例子
Image Not Showing
Possible Reasons
- The image file may be corrupted
- The server hosting the image is unavailable
- The image path is incorrect
- The image format is not supported
由圖中右方的節點可得
由上方的迴路可得
由下方的迴路可得
整理以上3式可得
如果轉換成矩陣型式
假設代入的數值分別為
電流分別為
矩陣運算
假設要解的未知數有
以上面的電路圖為例，
高中數學課本上比較常見的反矩陣求法，是將矩陣
先將第2、3列除以100可得
將第1列乘以
將第2列乘以
將第2列乘以
將第3列乘以
將第2列加到第1列可得
接下來再將
高中數學課本上比較常見的寫法是將要解的矩陣
先將第2、3列除以100可得
將第1列乘以
將第2列乘以
將第3列乘以
將第3列加到第1列，第3列乘以
將第2列除以3可得
最後將第2列加到第1列可得
使用 Python 的 NumPy 套件求解
NumPy 有現成的函式 numpy.linalg.solve 可以求解，也可以先用 numpy.linalg.inv 求反矩陣
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.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 |
