日期：2024年7月8日
原理
當電路上只有電池與電阻器，但又不是單純的串聯、並聯電路時，可以利用克希荷夫定律 (Kirchhoff's circuit laws) 求所有的分支電流。克希荷夫定律有兩個定則
- 節點定則：對於電路上的某個節點而言，流入、流出的電流量值相等，實際上就是電荷守恆。
- 迴路定則：對於電路上的某個迴路而言，繞一圈回到出發點時電位相等，實際上就是能量守恆。
由圖中右方的節點可得 $$ 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 |
參考資料
HackMD 版本連結：https://hackmd.io/@yizhewang/HkhhK8HDR
沒有留言:
張貼留言