FDM解常微分方程
問題描述
\[\frac{d^2\phi}{dx^2}=S_{\phi} \tag{1} \]
這是二階常微分方程（second-order Ordinary Differential Equation, ODE），考慮最簡單的情況即\(S=0\)，積分后可得\(\phi=c_1x+c_2\)，有兩個待定系數，因此要求解該方程必須提供兩個邊界條件（因為方程中不包含時間項，因此無初始條件），例如
\[\phi(x_L)=\phi_L \quad \phi(x_R)=\phi_R \]
場和邊界
我們把上述方程轉化為一個場（field）。試想存在一維（因為\(\phi\)僅與一個變量有關）直線，線上有若干等距分布的結點（node），每個結點都有唯一的\(\phi\)，那么\(\phi\)和\(x\)的關系滿足上述方程。
內部結點滿足方程，那么邊界上的\(\phi(x_L)\)和\(\phi(x_R)\)呢？事實上，邊界應是內部結點的延申，即邊界點也應該滿足上述方程。這也是為什么我們可以通過兩個邊界條件求解\(\phi=c_1x+c_2\)的原因，這隱含的假設就是邊界點滿足常微分方程。
結點上的值
用數值方法求解上述方程等價於求解場內每個結點上的\(\phi\)，結點\(i\)的上式表達為
\[\frac{d^2\phi}{dx^2}\bigg|_{i} \]
如何從結點值\(\phi_i\)得到導數值？很自然想到用Taylor展開。需要注意的是，Taylor展開隱含的假設是\(\phi\)無限可導。將\(i\)關於周圍兩點做Taylor展開，即
假設等距，則
\[\frac{d^2\phi}{dx^2}\bigg|_{i} =\frac{\phi_{i+1}+\phi_{i-1}-2\phi_i}{(\Delta x)^2}-\frac{(\Delta x)^2}{12} \frac{d^4\phi}{dx^4}\bigg|_{i} + ... \tag{2} \]
\[\varepsilon_i=-\frac{(\Delta x)^2}{12} \frac{d^4\phi}{dx^4}\bigg|_{i} + ... \]
其中\(\varepsilon_i\)為離散誤差或截斷誤差（discretization or truncation error），該誤差正比於距離的平方，因此我們稱上式對導數的逼近有二階精度。
邊界條件
上一部分是對場內部結點的推導。邊界上結點的值為邊界條件，有以下三種形式：
Dirichlet
\[\phi_{i=1}=\phi_L \]
Neumann
\[\frac{d\phi}{dx}\bigg|_{i=1}=J_L \tag{3} \]
將結點2關於結點1做Taylor展開，可以得到一階精度的梯度表達式
\[\frac{d\phi}{dx}\bigg|_{i=1}=\frac{\phi_2-\phi_1}{\Delta x}+\varepsilon(\Delta x) \tag{4} \]
將結點3關於結點1做Taylor展開，結合上式可以進一步得到二階精度
\[\frac{d\phi}{dx}\bigg|_{i=1}=\frac{4\phi_2-\phi_3-3\phi_1}{2\Delta x}+\varepsilon(\Delta x^2) \tag{5} \]
對於非Dirichlet的邊界條件，應盡量讓邊界條件也滿足內部結點的控制方程。
對\(\phi_2\)和\(\phi_3\)關於結點1展開可得
\[\phi_2=\phi_1+(\Delta x)\frac{d\phi}{dx}\bigg|_{i=1}+... \\ \phi_3=\phi_1+(2\Delta x)\frac{d\phi}{dx}\bigg|_{i=1}+... \]
將邊界條件公式(3)代入上式可得
\[\phi_2=\phi_1+(\Delta x)J_L+... \\ \phi_3=\phi_1+(2\Delta x)J_L+... \]
將上式整理可得結點1滿足的二階導數項：
\[\frac{d^2\phi}{dx^2}\bigg|_{i=1}=\frac{8\phi_2-\phi_3-7\phi_1-(6\Delta x)J_L}{2(\Delta x)^2}+\varepsilon(\Delta x^2) \tag{6} \]
上式為結點1滿足的控制方程，並用上了梯度邊界條件。需要說明的是，不用公式(6)而用公式(4,5)同樣可以求解問題，但是前者的精度更高。
Robin
\[\alpha \phi_1 + \beta \frac{d\phi}{dx}\bigg|_{i=1}=\gamma \]
類似地，令邊界結點1既滿足控制方程又滿足邊界條件。對上式移項可得
\[\frac{d\phi}{dx}\bigg|_{i=1}=\frac{\gamma-\alpha \phi_1}{\beta} \]
同樣代入\(\phi_2\)和\(\phi_3\)的Taylor展開，
\[\phi_2=\phi_1+(\Delta x) \left( \frac{\gamma-\alpha \phi_1}{\beta} \right) +... \\ \phi_3=\phi_1+(2\Delta x) \left( \frac{\gamma-\alpha \phi_1}{\beta} \right) +... \]
消去三階導數項並移項可得
\[\frac{d^2\phi}{dx^2}\bigg|_{i=1}=\frac{8\phi_2-\phi_3-7\phi_1-(6\Delta x)J_L}{2(\Delta x)^2}+\varepsilon(\Delta x^2) \tag{7} \]
其中
\[J_L=\frac{d\phi}{dx}\bigg|_{i=1}=\frac{\gamma-\alpha \phi_1}{\beta} \]
組建矩陣
有了內部結點的離散格式和邊界條件，我們便可對每個結點列方程。由於結點\(i\)離散格式中勢必會包含其他結點信息，例如對結點\(i\)
\[A_{i,1}\phi_1 + ... + A_{i,i-1}\phi_{i-1} + A_{i,i}\phi_i + A_{i,i+1}\phi_{i+1} + ... + A_{i,N}\phi_N =Q_i \]
將所有結點的方程組合起來，借助矩陣運算求解。
FDM求解泊松方程：二維問題
二維Poisson方程：
\[\nabla^2 \phi=\frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = S_{\phi} \tag{8} \]
如果\(S_{\phi}=0\)，則退化為Laplace方程。二階偏導數項同樣使用Taylor展開，只不過針對該問題應使用二維展開。如果使用正交網格結點，則相當於在每一維獨自做Taylor展開。
按上圖所示，網格結點可分為內部結點和邊界結點。
內部結點控制方程
由上文公式(2)可知，兩個二階偏導數項可分別寫為
\[\frac{\partial^2\phi}{\partial x^2}\bigg|_{i,j} =\frac{\phi_{i+1,j}+\phi_{i-1,j}-2\phi_{i,j}}{(\Delta x)^2} + \varepsilon(\Delta x^2) \tag{9} \]
\[\frac{\partial^2\phi}{\partial y^2}\bigg|_{i,j} =\frac{\phi_{i,j+1}+\phi_{i,j-1}-2\phi_{i,j}}{(\Delta y)^2} + \varepsilon(\Delta y^2) \tag{10} \]
去掉高階項，可得內部結點控制方程
\[\frac{\phi_{i+1,j}+\phi_{i-1,j}-2\phi_{i,j}}{(\Delta x)^2} + \frac{\phi_{i,j+1}+\phi_{i,j-1}-2\phi_{i,j}}{(\Delta y)^2} \approx S_{i,j} \tag{11} \]
下邊界
Robin邊界條件
\[\alpha \phi (x,0) + \beta \frac{\partial \phi}{\partial y}\bigg|_{x,0} = \gamma \]
對於下邊界的某個非邊角結點\((i,1)\)，其中\(i\neq 1, i\neq N\)，根據公式(7)來構建既滿足邊界條件又滿足內部結點控制方程的離散格式，即公式(7)
\[\frac{\partial^2\phi}{\partial y^2}\bigg|_{i,1} = \frac{8\phi_{i,2}-\phi_{i,3}-(7-6\Delta y \frac{\alpha}{\beta}) \phi_{i,1} - 6\Delta y (\gamma / \beta)}{2(\Delta y)^2} \tag{12} \]
而下邊界上\(x\)的二階偏導數項仍按照公式(9)離散。因此，下邊界結點的控制方程由公式(9)和(12)組成。
右邊界
Neumann邊界條件
\[\frac{\partial\phi}{\partial x}\bigg|_{L,y} = J_R \]
使用公式(6)得到\(x\)的二階偏導數項，即
\[\frac{\partial^2\phi}{\partial x^2}\bigg|_{N,j}=\frac{8\phi_{N-1,j}-\phi_{N-2,j}-7\phi_{N,j}-(6\Delta x)J_R}{2(\Delta x)^2} \tag{13} \]
而\(y\)的二階偏導數項仍用公式(10)。右下角的結點\((N,1)\)應同時滿足Robin和Neumann邊界條件，該結點的\(x\)偏導數項應使用公式(13)，\(y\)應使用公式(12)
上邊界和左邊界
Dirichlet邊界條件：
\[\phi(0,y)=\phi_L \\ \phi(x,H)=\phi_T \]
