數值常微分方程(BVP)

Relaxation Method

在使用 Finite Difference Method 解 BVP 時會扯到解矩陣的問題,如果想要避開解矩陣而使用遞迴的方式,Relaxation Method 就是這樣的一套思考邏輯。延續 Finite Difference Method 的三對角結構,假設我們要解的方程組是

\(a_{i-1}y_{i-1}+d_{i}y_{i}+c_{i+1}y_{i+1} = b_{i}\)

我們可以重寫這個式子為

\(y_i = \frac{b_i-a_{i-1}y_{i-1}-c_{i+1}y_{i+1}}{d_i}\)
若我們把式子左邊當作未來的時間,式子右邊當作過去的時間,我們就有一個遞迴公式可以期待新產生的 \(y_i\) 值,可以收斂到我們期待的解。這裡有兩個重要形式的迭代公式
- Jacobi iteration:使用第 \(j\) 次的資料去更新 \(j+1\) 次的資料

\(y_i^{(j+1)} = \frac{b_i-a_{i-1}y_{i-1}^{(j)}-c_{i+1}y_{i+1}^{(j)}}{d_i}\)

- Gauss-Seidel iteration:使用第 \(j\) 次以及計算過已經可以使用的 \(j+1\) 次的資料來更新第 \(j+1\) 次的資料

\(y_i^{(j+1)} = \frac{b_{i}-a_{i-1}y_{i-1}^{(j+1)}-c_{i+1}y_{i+1}^{(j)}}{d_i}\)


請注意當我們使用遞迴公式時,我們在實際運作上並沒有把矩陣給展開。因為每一個 \(y_i\) 只會和前一個時間點的 \(y_{i-1}\), \(y_i\) 以及\(y_{i+1}\) 有關。所以,我們只會存向量型的資料,並且用這資料來更新我們的 \(y_i\)。這裡展現的效率,是記憶體管理上的效率。我們應該都有一個經驗,當我們跑一個程式把 Ram 佔滿時會發生什麼情況?沒錯,你或聽到電腦的硬碟嘎嘎作響。因為這時一般地作業系統,會基於保護資料的原則,開啟電腦硬碟來模擬實體記憶體。倘若計算到需要使用虛擬記憶體時,整個計算時間就會嚴重拖慢。理由是硬碟的存取速度是遠遠小於 Ram 的存取速度,而 Ram 的存取速度,又小於 CPU 裡面的 caches 的速度。這就是為什麼疊代法迷人的地方。

其實關於遞迴式的寫法可以有很多種變形,重點是我們會把它看做是

\(y^{j+1} = A y^j + b = A^j y^0 + (A^{j-1}+A^{j-2}+\dots+A)b\)

的遞迴形式。這樣的形式如何保證會收斂呢?答案當然是對應的 \(A\) 矩陣每一項 eigenvalue 都小於 \(1\) 。上述的兩種方法, Gauss-Seidel 的方法會比 Jacobi 的方法來的有效率,然而若我們加入一個overcorrection 參數,我們可以得到更有效率的收斂公式

\(y_i^{(j+1)} = (1-w)y_i^j + w\frac{b_i-a_{i-1}y_{i-1}^{(j+1)}-c_{i+1}y_{i+1}^{(j)}}{d_i}\)

這個公式當 \(0<w<2\) 時會收斂。為甚麼稱為 Relaxation Method 就是基於解一個固定點(fixed point)問題,我們可以在遞迴式的兩邊加上\(\lambda y_i\),其中等號左邊的 \(\lambda y_i\) 設成未來時間\(\lambda y_i^{j+1}\),等號右邊的 \(\lambda y_i\) 設成過去時間\(\lambda y_i^j\)。這感覺好像只是玩一個恆等式遊戲,但是卻可以把這個 overcorrection 參數給引進來。多了這一個參數,我們就可以調整 \(w\) 使得對應的矩陣特徵值都小於 \(1\) 來保證收斂,甚至是影響收斂速度。

假設我們有一個猜想的BVP 解 \(\{y_i\}\),這個解帶入finite difference 的差分式時會產生一個誤差

(RL式一):\(e_i = \frac{y_{i+1}-2y_i+y_{i-1}}{h^2} + a_{1,i}\frac{y_{i+1}-y_{i-1}}{2h}+a_{0,i}y_i-u_i\)


我們想透過遞迴的方式修正這個 \(e_i\) 使得每一個 \(e_i\) 都是 \(0\)。這樣的 \(e_i\) 要如何修正呢?我們把每一個 \(e_i\) 看作是 \(\{ y_{i-1}, y_i, y_{i+1} \}\) 的函數。若對每一個 \(y_i\) 進行一個修正 \(\Delta y_i\),我們得到

\(e_i(y_{i-1}+ \Delta y_{i-1}, y_i+ \Delta y_i, y_{i+1}+ \Delta y_{i+1}) \simeq e_i + \frac{\partial e_i}{\partial y_{i-1}} \Delta y_{i-1} + \frac{\partial e_i}{\partial  y_i} \Delta y_i + \frac{\partial e_i}{\partial y_{i+1}} \Delta y_{i+1}\\ = e_i + \alpha_i \Delta y_{i-1} + \beta_i \Delta y_i + \gamma_i \Delta y_{i+1}\)

這裡的 \(\alpha_i\), \(\beta_i\) 和 \(\gamma_i\) 很容易用(RL式一)計算出來。