/* 取得增廣矩陣 * 輸入係數矩陣 M 及列向量 b 組合 * 回傳增廣矩陣 A */ let augmented = function(M, b) { let A = Array.from(M, function(v, k) { // 建立增廣矩陣 A return v.concat(b[k]); // 將 M 的每一列都加上 b 的列 }); return A; // 回傳增廣矩陣 A };
為了某些時候需要， 也可以撰寫將增廣矩陣 A 分解成係數矩陣 M 與列向量 b 的函數。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
/* 分解增廣矩陣 * 輸入增廣矩陣 A 分解 * 回傳為包含係數矩陣 M 及列向量 b 的物件 */ let unaugmented = function(A) { return { // 回傳包含 M 及 b 的物件 M: Array.from(A, function(r) { // 從增廣矩陣 A 建立 return r.filter(function(v, i) { // 過濾掉 A 每一列的最後一個元素則為係數矩陣 M return i != (r.length - 1); }); }), b: Array.from(A, function(r) { // 從增廣矩陣 A 建立 return r.filter(function(v, i) { // 過濾掉非 A 每一列的最後一個元素則為列向量 b return i == (r.length - 1); }); }) }; };
/* 兩列交換 * 輸入矩陣 A 交換 a 及 b 兩列 */ let exchange = function(A, a, b) { let T = Array.from(A[a]); // 複製第 a 列 A[a] = A[b]; // 令第 a 列為第 b 列 A[b] = T; // 令第 b 列為複製的 a 列 };
兩列加法
將一列乘上一個倍數加到另一列上。
1 2 3 4 5 6 7 8
/* 兩列加法 * 輸入矩陣 A 將第 a 列乘 scalar 加到第 b 列 */ let addition = function(A, a, b, scalar) { A[b] = A[b].map(function(v, k) { // 使第 b 列變化 return v + A[a][k] * scalar; // 將第 a 列乘 scalar 加上 }); };
單列縮放
將一列乘上一個常數。
1 2 3 4 5 6 7 8
/* 單列縮放 * 輸入矩陣 A 將第 a 列乘上 scalar */ let scalar = function(A, a, scalar) { A[a] = A[a].map(function(v, k) { // 將第 a 列變化 return v * scalar; // 將元素乘上 scalar }); };
消去法
通過剛剛寫的副程式，可以使消去法變得更容易完成， 但考慮到一些情形，必須要在撰寫一些函數。
取得軸
取得矩陣 A 的從第 k 列開始的軸。 考慮到如果當前軸那個位置的元素為零，這個函數可以幫我們找到可以交換的列。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
/* 取得軸 * 輸入矩陣 A 及第 k 列 * 建立一個 row(A) - k 長度的陣列 * 陣列元素值 e 為第 k 列開始第 e 行元素不為零 */ let findPivots = function(A, k) { returnArray.from({ length: (A.length - k) // 建立長度為 row(A) - k 的陣列 }, function(e, i) { for (let j = 0; j < A[k + i].length; j++) { // 遍歷 A 第 k + i 列的每個元素 if (A[k + i][j] != 0) { // 如果 A 第 k + i 列的第 j 個元素不為零 return j; // 返回 j } } return A[0].length; // 第 k + i 列的所有元素都為零返回行的數量 }); };
/* 高斯消去法 * 輸入增廣矩陣 A 使之變成上三角矩陣 */ letGauss = function(A) { let m = A.length; // A 有 m 個列 let n = A[0].length; // A 有 n 個行 let r = Math.min(m, n); // A 最大會有 r 個軸 for (let i = 0; i < m; i++) { // 對第 i 個列操作 let pivots = findPivots(A, i); // 計算從第 i 列開始的軸位置 let p = { // 建立目前軸 p 的位置 i: 0, j: 0 }; p.i = minIndex(pivots); // 取得最小軸位置的列 p.j = pivots[p.i]; // 取得最小軸位置 if (p.j == n) { // p.j 為 n 即整列是零 return; // 已經不能再做下去了 } exchange(A, i, i + p.i); // 交換最小軸位置的列到本列 scalar(A, i, 1 / A[i][p.j]); // 縮放本列的大小使軸為 1 for (let j = i + 1; j < m; j++) { // 遍歷軸下面的列 addition(A, i, j, -A[j][p.j] / A[i][p.j]); // 將其變為 0 } } };
/* 代入法 * 重複代入使軸上方元素為零 */ letJordan = function(A) { let m = A.length; // A 有 m 個列 let n = A[0].length; // A 有 n 個行 for (let i = m - 1; i >= 0; i--) { // 對第 i 個列操作 let pivots = findPivots(A, i); // 計算從第 i 列開始的軸位置 let p = { // 建立目前軸 p 的位置 i: null, j: pivots[0] }; if(p.j == n) { // p.j 為 n 即整列是零 continue; // 跳過這個列 } for (let j = i - 1; j >= 0; j--) { // 遍歷軸上面的列 addition(A, i, j, -A[j][p.j] / A[i][p.j]); // 將其變為 0 } } };