《Fluid Engine Development》 學習筆記1-求解線性方程組

我我的對基於物理的動畫很感興趣，最近在嘗試閱讀《Fluid Engine Development》，因爲內容涉及太多的數學問題，而單純學習數學又過於枯燥，難以堅持學習（我中途放棄好屢次了），打算嘗試經過編寫博客總結知識的學習方法來學習。

在計算數值問題時，咱們常常遇到線性方程，好比基於網格的流體模擬在求解擴散和壓強，須要求解線性方程組。

線性方程組

線性方程組 \(\left\{ \begin{matrix} 2 * x - y =3 \\ -x + 2 * y = 6 \end{matrix}\right\}\)​

能夠轉換成形如 \(A \times x=b\) 的矩陣矢量形式，轉換結果以下

$\left[
\begin{matrix}
2 & -1 \
-1 & 2
\end{matrix}
\right] \left[
\begin{matrix}
x \
y
\end{matrix}
\right] = \left[
\begin{matrix}
3 \
6
\end{matrix}
\right] $

其中A是系統矩陣，x 爲所求解，b是一個矢量，也就是線性方程的常數項

至於如何求解x，一般咱們是經過 是讓等式兩邊乘以系統矩陣的逆矩陣

高斯消元法

求解線性方程組咱們一般使用高斯消元法來求解逆矩陣，這種方法雖然足夠直接，然而把高斯消元法做爲一種算法來看待，這種算法的時間複雜度達到了驚人的 \(O(n^3)\),n是矩陣的尺寸，由此咱們能夠明確高斯消元法並不適合有着許多數值問題的較大系統。

若是不選用高斯消元法直接求解逆矩陣，若是不經過高斯消元法，咱們又該經過什麼方法計算逆矩陣呢？《Fluid Engine Development》給出了就4種不那麼直接的方法。經過對解的猜想和屢次迭代獲得近似解

雅可比方法

雅克比方法（Jacobi方法）是用於肯定對角佔優的解的迭代算法。求解每一個對角元素，並插入近似值。而後迭代該過程直到它收斂

使用雅克比方法，須要將\(A \times x=b\) 轉換成 \(（D + R） \times x=b\) （矩陣A被拆成對角矩陣D和矩陣R）

因此易得解\({x}^{(k+1)} = D^{-1} (\mathbf{b} - R \mathbf{x}^{(k)})\) k是迭代的次數

一樣做爲矢量的解的每個元素能夠經過$ x^{(k+1)}i = \frac{1}{a{ii}} \left(b_i -\sum_{j\ne i}a_{ij}x^{(k)}_j\right),\quad i=1,2,\ldots,n.$計算獲得

Gauss-Seidel法

Gauss-Seidel方法和jacobi方法有些像

將\(A \times x=b\) 轉換成 \(（L + U） \times x=b\) （矩陣A被拆成包含對角部分的下三角形和又矩陣A上三角形部分構成的矩陣U）

具體以下

$\left[
\begin{matrix}
2 & -1 \
-1 & 2
\end{matrix}
\right]= \left[
\begin{matrix}
2 & 0 \
-1 & 2
\end{matrix}
\right]+ \left[
\begin{matrix}0 & -1 \
0 & 0
\end{matrix}
\right] $

經過這種形式的轉換咱們能夠輕易的發現解的第一個元素能輕易的被計算出來

$ x^{(k+1)}1 = \frac{1}{a{11}} \left(b_1 -\sum_{j>1}a_{1j}x^{(k)}_j\right)$

將上式代入第二行，咱們能夠獲得$ x^{(k+1)}2 = \frac{1}{a{22}} \left(b_i -a_{21}x_{1} ^{k+1}-\sum_{j>1}a_{1j}x^{(k)}_j\right)$

而後能夠獲得一樣的迭代解 \(x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j=1}^{i-1}a_{ij}x^{(k+1)}_j - \sum_{j=i+1}^{n}a_{ij}x^{(k)}_j \right),\quad i=1,2,\dots,n.\)

梯度降低法

梯度降低法將線性方程組求解問題轉換成求解最小值問題 \(F(x) = |Ax - b|^{2}\)

若是x有解，則F(x) = 0,若是無解，則能夠一直迭代直到x有解，

若是從x1開始，沿着梯度降低，逐步迭代來靠近零點

共軛梯度法

相對梯度降低法，共軛梯度法不使用梯度方向迭代，而是使用共軛方向

何爲共軛，若是 \(a *( A b) = 0\) 則咱們稱向量a，b關於矩陣A 共軛

預處理共軛梯度法

在此基礎上，咱們能夠進一步加速： 這個方法稱爲預處理共軛梯度法。大致的思想是，引入一個preconditioning filter到系統中:

\(M^{-1}Ax = M^{-1}b\)

具體算法實現以下

d就是新的共軛方向，alpha的推導參考連接