有限差分法

有限差分法，係一種攞嚟解微分方程嘅數值方法。

起源

自從微積分現世後，好多科學上嘅定理都可以用微分方程寫出嚟。為咗明白一條微分方程所表達嘅現象，解方程就係好重要嘅一環。但係，隨之而嚟嘅問題就係，呢個世界上大部分嘅微分方程都冇初等函數解，好多甚至連解析解都冇。所以，以前嘅人通常會簡化呢啲微分方程，例如用尺度分析鏟走啲冇咁重要嘅項，或者將啲方程線性化之類，然後就可以用一般嘅方法嚟搵方程嘅解析解。不過依然都有啲問題，就係處理太多初始條件同埋邊界條件等等大量嘅數據嘅時候，用人手計數實在太麻煩。

隨著電腦嘅發明同埋普及，解微分方程呢啲複雜嘅嘢就當然慢慢交俾運算量比人手大好多嘅電腦做，可以慳返好多時間。但係，電腦只可以處理離散嘅數據，單係入初始條件同埋邊界條件就要用喺唔同嘅格點輸入數據，而唔係將一個連續嘅函數沊落去就算。而且，電腦唔係人腦，唔識計微積分，只識計加減乘除。所以，我哋有必要將呢啲微分方程用基本代數嘅運算法則嚟估算佢，然後將呢啲只有四則運算嘅公式變成一條程式碼，放落電腦度 run 個解出嚟。呢啲公式就係有限差分方程。

函數嘅有限差分表示格式

喺有限差分裏面，函數唔係連續嘅，而係喺每一個格點上面對應一個數值。所以，我哋要用第二種方式嚟表示離散化咗嘅函數。

首先，我哋要用一啲符號或者字母代表三維空間同埋一維時間四個維度。習慣上，$n$ 代表時間維度上嘅格點指數，$i$ 代表沿著 $x$ 軸嘅格點指數，$j$ 代表沿著 $y$ 軸嘅格點指數，$k$ 代表沿著 $z$ 軸嘅格點指數。假設一個 $(x,y,z,t)$ 喺正一個格點，而呢個格點嘅 $x$ 方向指數、$y$ 方向指數，$z$ 方向指數同埋時間方向指數分別係 $i$、$j$、$k$ 同埋 $n$，噉一個定義域喺三維空間同埋一維時間嘅函數 $u$ 嘅有限差分表示法就係：

$u_{i,j,k}^{n}=u(x,y,z,t)$

留意時間指數永遠寫喺右上角，而空間指數永遠寫喺右下角。而且，呢度假設四個維度嘅格點分別按照 $\Delta x$、$\Delta y$、$\Delta z$ 同埋 $\Delta t$ 嘅間隔平均分佈。

如果個函數嘅域唔覆蓋哂全部四個維度，噉對應嘅指數可以省略唔寫。例如如果一個函數只對應一維空間同埋一維時間（即係 $u(x,t)$），佢嘅有限差分表示法就係 $u_{i}^{n}$；如果連時間維度都冇埋，就係 $u(x)=u_{i}$；如果剩係有一維時間，就係 $u(t)=u^{n}$。我哋要小心千祈唔好將呢兩個寫法分別同數列同埋次方炒埋一碟！

導數喺有限差分嘅表示法

一堆有限差分法裏面一階同二階導數唔同嘅表示式，除咗準確到唔同階之外，仲有分向前差分、向後差分同埋中央差分三種格式。

利用導數嘅基本原理，我哋可以得到有限差分法裏面最簡單嘅公式。因為一階導數嘅定義係：

$f'(x)=\displaystyle \lim _{\Delta x\rightarrow 0}{\frac {f(x+\Delta x)-f(x)}{\Delta x}}$

對於一階導數其中一條有限差分公式就係：

$f'(x)\approx {\frac {f(x+\Delta x)-f(x)}{\Delta x}}$

我哋再將 $f(x+\Delta x)$ 變做主項，就有：

$f(x+\Delta x)\approx f(x)+f'(x)\Delta x$

或者，用格點表示法就係：

${\frac {df}{dx}}\approx {\frac {f_{i+1}-f_{i}}{\Delta x}}$

所以，當我哋知道喺某一個格點 $f(x)$ 同埋 $f'(x)$ 嘅數值，就可以用自己定義嘅格點距離搵下一個格點上面 $f$ 嘅數值。歐拉法就係用緊呢個概念。

同樣地，如果將 $\Delta x$ 變做 $-\Delta x$，$\Delta x\rightarrow 0$ 都代表 $-\Delta x\rightarrow 0$，所以導數嘅基本原理亦都可以寫成：

$f'(x)=\displaystyle \lim _{\Delta x\rightarrow 0}{\frac {f(x)-f(x-\Delta x)}{\Delta x}}$

由呢度可以推出另一條一階導數嘅差分公式，就係：

$f'(x)\approx {\frac {f(x)-f(x-\Delta x)}{\Delta x}}$

我哋再將 $f(x)$ 變做主項，就有：

$f(x)\approx f(x-\Delta x)+f'(x)\Delta x$

或者，用格點表示法就係：

${\frac {df}{dx}}\approx {\frac {f_{i}-f_{i-1}}{\Delta x}}$

呢個概念同歐拉法都幾似，不過係用緊相反方向嘅導數，所以叫做反向歐拉法。

呢兩種有限差分方法都係準確至一階。

我哋仲可以用到上面嘅兩條差分公式推到另一條差分公式出嚟。通過整合上面得到嘅兩條公式：

$f'(x)\approx {\frac {f(x+\Delta x)-f(x)}{\Delta x}}\cdots (1)$

$f'(x)\approx {\frac {f(x)-f(x-\Delta x)}{\Delta x}}\cdots (2)$

我哋可以推出：

$(1)+(2):$

$2f'(x)\approx {\frac {f(x+\Delta x)-f(x)}{\Delta x}}+{\frac {f(x)-f(x-\Delta x)}{\Delta x}}$

$2f'(x)\approx {\frac {f(x+\Delta x)-f(x)+f(x)-f(x-\Delta x)}{\Delta x}}$

$2f'(x)\approx {\frac {f(x+\Delta x)-f(x-\Delta x)}{\Delta x}}$

$f'(x)\approx {\frac {f(x+\Delta x)-f(x-\Delta x)}{2\Delta x}}$

或者，用格點表示法就係：

${\frac {df}{dx}}\approx {\frac {f_{i+1}-f_{i-1}}{2\Delta x}}$

呢條新嘅公式喺中央差分法裏面好有用，而且比起先前兩條公式誤差更加細，因為佢係準確到二階。

差分法嘅公式仲可以由泰勒級數推出嚟。首先將佢應用到 $f(x+\Delta x)$，我哋有：

$f(x+\Delta x)=\displaystyle \sum _{n=0}^{\infty }{\frac {f^{(n)}(x)}{n!}}\Delta x^{n}=f(x)+f'(x)\Delta x+{\frac {f''(x)}{2!}}\Delta x^{2}+{\frac {f'''(x)}{3!}}\Delta x^{3}+{\frac {f''''(x)}{4!}}\Delta x^{4}+\cdots \,\,\,\,\,\cdots (3)$

再將佢應用到 $f(x-\Delta x)$，就有：

$f(x-\Delta x)=\displaystyle \sum _{n=0}^{\infty }{\frac {(-1)^{n}f^{(n)}(x)}{n!}}\Delta x^{n}=f(x)-f'(x)\Delta x+{\frac {f''(x)}{2!}}\Delta x^{2}-{\frac {f'''(x)}{3!}}\Delta x^{3}+{\frac {f''''(x)}{4!}}\Delta x^{4}-\cdots \,\,\,\,\,\cdots (4)$

將 $(3)$ 準確到一階，我哋可以重寫佢，變成：

$f(x+\Delta x)=f(x)+f'(x)\Delta x+O(\Delta x^{2})$

其中用大O符號嚟代表誤差。我哋抹走 $O(\Delta x^{2})$，就會彈返歐拉法條式出嚟。因為 $O(\Delta x^{2})$ 代表二階或以上嘅誤差，最大嘅誤差都係去到二階，所以可以睇得出呢條式喺一階嘅領域裏面絕對準確。

用同樣嘅方法，我哋重寫 $(4)$ 令到佢準確到一階，就會有：

$f(x-\Delta x)=f(x)-f'(x)\Delta x+O(\Delta x^{2})$

用返上面嘅睇法，我哋可以睇到呢條式都係喺一階嘅領域裏面絕對準確。

如果要準確到二階，我哋用 $(3)-(4)$，就會有：

$f(x+\Delta x)-f(x-\Delta x)=2f'(x)\Delta x+{\frac {2f'''(x)}{3!}}\Delta x^{3}+\cdots =2f'(x)\Delta x+O(\Delta x^{3})$

我哋抹走 $O(\Delta x^{3})$，就會彈返中央差分法條式出嚟。因為 $O(\Delta x^{3})$ 代表三階或以上嘅誤差，最大嘅誤差都係去到三階，所以可以睇得出呢條式喺二階嘅領域裏面絕對準確。

至於二階導數，我哋可以考慮佢係一階導數嘅導數嚟推出差分公式。首先，我哋用中央差分版嘅導數基本原理嚟表示二階導數：

$f''(x)=\displaystyle \lim _{\Delta x\rightarrow 0}{\frac {f'(x+\Delta x)-f'(x-\Delta x)}{2\Delta x}}$

拆走個極限，再將個格點間隔由 $2\Delta x$ 變返做 $\Delta x$（因為 $\Delta x\rightarrow 0$ 代表 ${\frac {\Delta x}{2}}\rightarrow 0$），我哋有：

$f''(x)\approx {\frac {f'(x+{\frac {\Delta x}{2}})-f'(x-{\frac {\Delta x}{2}})}{\Delta x}}$

我哋再對 $f'(x+{\frac {\Delta x}{2}})$ 同埋 $f'(x-{\frac {\Delta x}{2}})$ 用格點間隔係 $\Delta x$ 嘅中央差分法，就會有：

$f'(x+{\frac {\Delta x}{2}})\approx {\frac {f(x+\Delta x)-f(x)}{\Delta x}}$ （條式同歐拉法一樣）

$f'(x-{\frac {\Delta x}{2}})\approx {\frac {f(x)-f(x-\Delta x)}{\Delta x}}$ （條式同反向歐拉法一樣）

將佢哋代返入去，就有：

$f''(x)\approx {\frac {{\frac {f(x+\Delta x)-f(x)}{\Delta x}}-{\frac {f(x)-f(x-\Delta x)}{\Delta x}}}{\Delta x}}$

$f''(x)\approx {\frac {\frac {f(x+\Delta x)-f(x)-[f(x)-f(x-\Delta x)]}{\Delta x}}{\Delta x}}$

$f''(x)\approx {\frac {f(x+\Delta x)-2f(x)+f(x-\Delta x)}{\Delta x^{2}}}$

或者，用格點表示法就係：

${\frac {d^{2}f}{dx^{2}}}\approx {\frac {f_{i+1}-2f_{i}+f_{i-1}}{\Delta x^{2}}}$

呢條就係二階導數嘅差分公式，準確到二階。

用泰勒級數都可以得到呢個結果。我哋考慮 $(3)+(4)$：

$f(x+\Delta x)+f(x-\Delta x)=2f(x)+2\left[{\frac {f''(x)}{2!}}\Delta x^{2}\right]+2\left[{\frac {f''''(x)}{4!}}\Delta x^{4}\right]+\cdots$

$f(x+\Delta x)+f(x-\Delta x)=2f(x)+2\left[{\frac {f''(x)}{2!}}\Delta x^{2}\right]+O(\Delta x^{4})$

當 $O(\Delta x^{4})\rightarrow 0$，二階導數嘅差分公式就會出返嚟。

呢部分最頂嗰幅圖有一堆唔同導數嘅有限差分法公式，佢哋嘅證明方式同上面大致一樣。

有限差分法唔同嘅格式

上面講到，有限差分法大致上可以分做向前差分、向後差分同埋中央差分三種格式。但係，通常啲微分方程會同時覆蓋時間同埋空間嘅維度，我哋對時間做差分同埋對空間做差分其實可以用唔同嘅格式。所以，喺時間同埋空間上面都可以獨立細分三種格式，而總共就會衍生出九種格式。呢九種格式分別係時間向前差分空間向前差分（FTFS）格式、時間向前差分空間中央差分（FTCS）格式、時間向前差分空間向後差分（FTBS）格式、時間中央差分空間向前差分（CTFS）格式、時間中央差分空間中央差分（CTCS）格式、時間中央差分空間向後差分（CTBS）格式、時間向後差分空間向前差分（BTFS）格式、時間向後差分空間中央差分（BTCS）格式同埋時間向後差分空間向後差分（BTBS）格式。

有限差分法嘅應用

唔同嘅有限差分法格式可以用嚟搵唔同類型嘅微分方程嘅數值解，佢哋大致上嘅做法都係將條方程裏面唔同階嘅導數全部塞哂啲有限差分法嘅樣落去，然後再用加減乘除逐格逐格同埋逐個時間點搵哂所有格點嘅數值出嚟。例如解一條微分方程 $y{\frac {d^{2}y}{dx^{2}}}+{\frac {dy}{dx}}-{\frac {1}{y}}=1$，驟眼望落去呢條方程好難搵解析解。我哋用中間差分法將佢重寫做：

$y_{i}\left({\frac {y_{i+1}-2y_{i}+y_{i-1}}{\Delta x^{2}}}\right)+{\frac {y_{i+1}-y_{i-1}}{2\Delta x}}-{\frac {1}{y_{i}}}=1$

為咗搵前面嗰格嘅數值，我哋要將 $y_{i+1}$ 變做主項，條式就會慢慢變成：

${\frac {y_{i+1}-2y_{i}+y_{i-1}}{\Delta x^{2}}}={\frac {1}{y_{i}}}\left(1+{\frac {1}{y_{i}}}-{\frac {y_{i+1}-y_{i-1}}{2\Delta x}}\right)$

${\frac {y_{i+1}-2y_{i}+y_{i-1}}{\Delta x^{2}}}={\frac {1}{y_{i}}}+{\frac {1}{y_{i}^{2}}}-{\frac {y_{i+1}-y_{i-1}}{2y_{i}\Delta x}}$

$y_{i+1}\left({\frac {1}{\Delta x^{2}}}+{\frac {1}{2y_{i}\Delta x}}\right)={\frac {2y_{i}-y_{i-1}}{\Delta x^{2}}}+{\frac {1}{y_{i}}}+{\frac {1}{y_{i}^{2}}}+{\frac {y_{i-1}}{2y_{i}\Delta x}}$

${\frac {y_{i+1}}{\Delta x}}\left({\frac {1}{\Delta x}}+{\frac {1}{2y_{i}}}\right)={\frac {2y_{i}-y_{i-1}}{\Delta x^{2}}}+{\frac {1}{y_{i}}}+{\frac {1}{y_{i}^{2}}}+{\frac {y_{i-1}}{2y_{i}\Delta x}}$

$y_{i+1}={\frac {2y_{i}\Delta x^{2}}{2y_{i}+\Delta x}}\left({\frac {2y_{i}-y_{i-1}}{\Delta x^{2}}}+{\frac {1}{y_{i}}}+{\frac {1}{y_{i}^{2}}}+{\frac {y_{i-1}}{2y_{i}\Delta x}}\right)$

$y_{i+1}={\frac {2y_{i}\Delta x^{2}}{2y_{i}+\Delta x}}\left[{\frac {2y_{i}-y_{i-1}}{\Delta x^{2}}}+{\frac {1}{y_{i}}}\left(1+{\frac {1}{y_{i}}}+{\frac {y_{i-1}}{2\Delta x}}\right)\right]$

將呢條式用喺定義域裏面嘅每個格點，然後逐格逐格推，就可以搵到喺成個定義域裏面呢條微分方程嘅解。

喺大部分嘅時候，有限差分法通常係用嚟解偏微分方程，例如移流方程、熱方程、波動方程、拉普拉斯方程，甚至納維-斯托克斯方程之類，因為佢哋喺好多科學現象，尤其係物理現象裏面會出現。至於要用唔同嘅格式，係因為唔同種類嘅偏微分方程喺唔同嘅格式之下會有唔同嘅數值穩定性。

就用一維空間嘅移流方程做例子，其中一個可以用嘅格式係 CTCS 格式。條方程原本嘅樣係：

${\frac {\partial u}{\partial t}}+c{\frac {\partial u}{\partial x}}=0$

用 CTCS 格式個樣就會變成：

${\frac {u_{i}^{n+1}-u_{i}^{n-1}}{2\Delta t}}+c{\frac {u_{i+1}^{n}-u_{i-1}^{n}}{2\Delta x}}=0$

再將 $u_{i}^{n+1}$ 變成主項就會有：

$u_{i}^{n+1}=u_{i}^{n-1}-{\frac {c\Delta t}{\Delta x}}(u_{i+1}^{n}-u_{i-1}^{n})$

呢條式就可以直接入落電腦程式度解咗佢。