2.1 高斯消元
數(shù)學(xué)原理
??雷米茲算法的基礎(chǔ)是N元一次方程的求解。此外,離散傅里葉變換也用到了N元一次方程的求解。當(dāng)然,實(shí)際應(yīng)用中求解多元一次方程組的例子非常多,我就不一一列舉。求解的方法主要有兩種,一種是高斯消元,一種是遞歸法。
??遞歸法不推薦,已經(jīng)被我刪除。而高斯消元不僅僅是用來(lái)解方程,還可以用來(lái)求行列式、逆矩陣。
??高斯消元是利用矩陣表示N元一次方程組。然后將矩陣變成上三角矩陣或下三角矩陣。
??所謂上三角矩陣就是左下角全部為0的矩陣,下三角反之。
??如下,就是轉(zhuǎn)換過(guò)程:
[111110211111?121110?321214]\left[ \begin{matrix} 1 & 1 & 1 & 1 & 10\\ 2 & 1 & 1 & 1 & 11\\ -1 & 2 & 1 & 1 & 10\\ -3 & 2 & 1 & 2 & 14\\ \end{matrix} \right]\\ ?12?1?3?1122?1111?1112?10111014??
??將第一列變成0:
[111110011190?3?2?2?200?5?4?5?44]\left[ \begin{matrix} 1&1&1&1&10\\ 0&1&1&1&9\\ 0&-3&-2&-2&-20\\ 0&-5&-4&-5&-44\\ \end{matrix} \right]\\ ?1000?11?3?5?11?2?4?11?2?5?109?20?44??
??將第二列變成0:
[1111100111900?1?1?700?10?1]\left[ \begin{matrix} 1&1&1&1&10\\ 0&1&1&1&9\\ 0&0&-1&-1&-7\\ 0&0&-1&0&-1\\ \end{matrix} \right]\\ ?1000?1100?11?1?1?11?10?109?7?1??
??將第三列變成0:
[1111100111900?1?1?700016]\left[ \begin{matrix} 1&1&1&1&10\\ 0&1&1&1&9\\ 0&0&-1&-1&-7\\ 0&0&0&1&6\\ \end{matrix} \right] ?1000?1100?11?10?11?11?109?76??
??這個(gè)轉(zhuǎn)換過(guò)程,就是先用第一行和后續(xù)的所有行進(jìn)行計(jì)算,讓后續(xù)所有行以0開始。
??然后將右下角視為一個(gè)子矩陣,繼續(xù)這個(gè)運(yùn)算。這個(gè)算法可以用遞歸實(shí)現(xiàn),也可以用循環(huán)實(shí)現(xiàn)。但是最好實(shí)用循環(huán)實(shí)現(xiàn),也就是高斯消元方法。
??然后使用代入法就可以解開,但是這時(shí)候千萬(wàn)不要去再轉(zhuǎn)換為下三角矩陣,并且進(jìn)一步轉(zhuǎn)為單位矩陣。因?yàn)槟菢訒?huì)有很多不必要的運(yùn)算。
轉(zhuǎn)換為三角矩陣
??這個(gè)核心的行變換方法其實(shí)就是行倍加。假設(shè)第一行的系數(shù)是aaa,第二行的系數(shù)是bbb,那么就是把第一行乘以?ba-\fraca?ab?加到第二行上。
??行例如第一行是x+2y+3z=6x+2y+3z =6x+2y+3z=6,第二行是3x?3y+9z=93x-3y+9z=93x?3y+9z=9,那么第一行就乘以?31=?3-\frac31=-3?13?=?3,然后加到第二行。
??代碼如下:
??而python代碼則如此:
def to_upper_triangle(self):# 將矩陣變成上三角矩陣lines = copy.deepcopy(self.__lines)for i in range(len(lines) - 1):for j in range(i + 1, len(lines)):Matrix.row_operate(lines[i], lines[j], lines[j][i], lines[i][i], i, len(lines[i]))return Matrix(lines)行變換代碼
@staticmethod def row_operate(line1, line2, coefficient1, coefficient2, start_column, end_column):for i in range(start_column, end_column):line2[i] = line1[i] * coefficient1 - line2[i] * coefficient2代入法
??而獲取了上三角矩陣之后。
??有兩種方法可以解開方程,一種辦法是代入法,另一種方法是從下到上使用行倍加變換。代入法性能高一點(diǎn),性能高的原因是避免了不必要的乘以0的運(yùn)算。
??就是套進(jìn)去解開了。解開的方法代碼如下:
??其具體的思路是用最后一行代進(jìn)去得到一個(gè)數(shù)字,比如如下的上三角矩陣:
[1111100111900?1?1?700016]\left[ \begin{matrix} 1&1&1&1&10\\ 0&1&1&1&9\\ 0&0&-1&-1&-7\\ 0&0&0&1&6 \end{matrix} \right] ?1000?1100?11?10?11?11?109?76??
??從最后一行開始,得到x4=6x_4=6x4?=6,然后x4=6x_4=6x4?=6代入上面的每一行,矩陣變成了這樣
[11104011030010100016]\left[ \begin{matrix} 1&1&1&0&4\\ 0&1&1&0&3\\ 0&0&1&0&1\\ 0&0&0&1&6\\ \end{matrix} \right] ?1000?1100?1110?0001?4316??
??然后把x3x_3x3?求出來(lái),代入上面的每一行,最終,整個(gè)矩陣變成
[10001010020010100016]\left[ \begin{matrix} 1&0&0&0&1\\ 0&1&0&0&2\\ 0&0&1&0&1\\ 0&0&0&1&6 \\ \end{matrix} \right] ?1000?0100?0010?0001?1216??
??所以這個(gè)方程的解就是:
x1=1x2=2x3=1x4=6x_1=1\\ x_2=2\\ x_3=1\\ x_4=6\\ x1?=1x2?=2x3?=1x4?=6
最終代碼
??代入法計(jì)算性能更高,如果再轉(zhuǎn)為下三角會(huì)有很多非必要的預(yù)算,再轉(zhuǎn)為單位矩陣,又是一堆沒(méi)必要的運(yùn)算。高斯消元總體代碼如下:
public T[] gaussianReduction() {final Matrix<T> tMatrix = this.toUpperTriangular();// 進(jìn)行代入法變換。T[] solution = newArray(this.array.length);final T[][] ts = tMatrix.copyArray();final int n = ts.length;System.out.println();System.out.println(tMatrix.toBeautifulString());// 從最后一行出來(lái)-for (int i = n - 1; i >= 0; i--) {// 首先計(jì)算出自己solution[i] = divide(ts[i][n], ts[i][i]);// 遍歷之上的所有行for (int j = 0; j < i; j++) {ts[j][n] = subtract(ts[j][n], multiply(ts[j][i], solution[i]));}}return solution; }??python代碼如下:
def gaussian_reduction(self):matrix = self.to_upper_triangle()solution = [0 for _ in matrix.__lines]n = len(matrix.__lines)for i in range(n - 1, -1, -1):solution[i] = matrix.__lines[i][n] / matrix.__lines[i][i]for j in range(0, i):matrix.__lines[j][n] = matrix.__lines[j][n] - matrix.__lines[j][i] * solution[i]return solution總結(jié)
- 上一篇: 国内9大免费CDN汇总,除了加速乐,你还
- 下一篇: 数据分析师要掌握SQL到什么程度?