數值格式誤差以及收斂精度估計方法

目錄

• 數值格式誤差以及收斂精度估計方法
  • 1.簡介
  • 2.收斂精度 (階) 介紹
  • 3.收斂精度計算
  • 4.格式誤差計算
    • 4.1.范數誤差定義
    • 4.2.選取正確的范數誤差函數
    • 4.3.范數誤差計算表達式
  • Reference

1.簡介

隨着高精度格式越來越多應用到CFD中，如何判斷數值格式的收斂精度 (Converge Rate) 也逐漸成為一個重要問題。

2.收斂精度 (階) 介紹

當我們考慮采用數值方法計算一個精確解\(u\)時,數值解\(\tilde{u}_h\)與精確解近似程度一般和一個參數\(h\)相關，這個近似程度可以表示為

\[\begin{equation} \left| \tilde{u}_h - u \right| \le Ch^p \end{equation}\]

其中\(C\)是與\(h\)無關的常數。在式（1）中的冪次\(p\)就是我們常說的數值格式收斂精度（階數）。誤差\(\left| \tilde{u}_h - u \right|\)也可用\(h\)表示為

\[\begin{equation} \tilde{u}_h - u = Ch^p + O(h^{p+1}) \end{equation}\]

在判斷CFD數值模型結果精度時，需要比較的不再僅僅是兩個標量，而是數值解\(\tilde{u}_h(x)\)與精確解\(u(x)\)兩個函數之間插值，此時就需要引入泛函中范數的概念。如何計算誤差函數對應的范數將在第四節討論。一般而言，若數值格式具有\(p\)階空間與\(q\)階時間精度，那么殘差項應該滿足

\[\begin{equation} R(\Delta x, \Delta t) = E\left( \Delta x^p, \Delta t^q \right) + O(\Delta x^{p+1}, \Delta t^{q+1}) \end{equation}\]

其中\(E\)為線性方程，而\(p\)、\(q\)中最小值\({\mathrm{min}(p, q)}\)決定了格式收斂階大小。

3.收斂精度計算

通常而言，在計算NS方程時，首先會把方程在空間內進行離散，得到常微分方程（ODE）

\[\begin{equation} \frac{\partial u}{\partial t} = L(u) \end{equation}\]

而求解只包含時間的常微分方程時，通常可以采用具有更高的精度格式\(q>p\)來減小誤差（Runge-Kutta，AB）。因此，數值格式的精度主要受空間離散格式的限制。

為了求出空間離散格式收斂階\(p\)，可以根據方程（2）采用網格逐次加密方法。當我們采用不同尺寸的網格計算，可以得到誤差函數隨網格尺寸\(h\)的變化關系，如分別采用\((\Delta x, \Delta t), \frac{1}{2}(\Delta x, \Delta t), \frac{1}{4}(\Delta x, \Delta t)\)的步長進行計算，得到\(R(\Delta x), R(\frac{1}{2}\Delta x), R(\frac{1}{4} \Delta x)\)，由於

\[\begin{equation} R(\Delta x) = E(\Delta x) + O(\Delta x^{p+1}) \approx C \Delta x^p \end{equation}\]

那么

\[\begin{equation} \frac{R( \Delta x )}{R( \frac{1}{2}\Delta x )} = \frac{C\Delta x^p}{ C \left( \frac{\Delta x}{2} \right)^p } \end{equation}\]

因此收斂精度\(p\)可以采用下式進行估計

\[p = \frac{log \left( { R( \Delta x )}/{R( \frac{1}{2}\Delta x ) }\right) }{log2} \]

4.格式誤差計算

4.1.范數誤差定義

首先介紹范數概念。范數是表示是泛函空間內兩個元素距離的函數，在泛函空間內每個函數\(u(x)\)都是一個元素，而\(L_1\)、\(L_2\)和\(L_{\infty}\)等范數表示就是兩個元素之間距離。

\(L_1\)范數通常稱為最小絕對偏差LAD（least absolute deviations）或最小絕對誤差LAE（least absolute errors）。它統計目標函數\(y(x)\)與期望函數\(f(x)\)絕對誤差之和

\[\begin{equation} S = \int \left|y(x_i) - f(x_i) \right| \rm{dx} \end{equation}\]

\(L_2\)范數也被稱為最小二乘。它通常統計目標函數與期望函數平方誤差之和^[1]

\[\begin{equation} S = \int \left|y(x_i) - f(x_i) \right|^2 \rm{dx} \end{equation}\]

\(L_{\infty}\)范數用來統計目標函數與期望函數之間最大誤差

\[\begin{equation} S = \rm{max} \left|y(x_i) - f(x_i) \right| \end{equation}\]

4.2.選取正確的范數誤差函數

想要回答什么時候選取哪個范數誤差進行估計十分復雜，在CFD計算中，需要根據算例以及各個范數性質選取不同的誤差范數。

在CFD計算中，選取誤差范數需要考慮特性包括：^[2]

1. 魯棒性^[3]
   L1-范數誤差統計的數據包含異常值時，可以將數據中的異常值安全而有效地忽略。所以L1-范數適合評價具有只有極少特別大誤差出現的間斷解算例；
   而L2-范數會將誤差平方，總誤差會增加很多，所以L2誤差對於異常值更加敏感。因此L2-范數常用來評價幾乎不存在特別大誤差的連續解算例。
2. 穩定性
   L1-范數具有不穩定性，精確值一個較小的水平移動會導致回歸線產生特別大的改變。因此若算例沒有解析解，需要采用數值方法估計收斂精度時，采用L2-范數更合適。

Wolfgang Bangerth.
Professor of mathematics at Texas A&M University. Research focus is on numerical methods for partial differential equations, finite element software, the deal.II software library.

衡量偏微分方程求解誤差時候，自然地選擇是解所在的空間范數，因為范數作用就是衡量空間內兩元素之間距離。舉例來說，對於橢圓偏微分方程，方程解所在空間是\(H^1\)，所以使用\(H^1\)范數來估計誤差就是很合適的。這樣選取的意義在於，以上面例子來說，真實解不在\(W^{1,\infty}\)空間內，因此計算梯度的最大誤差就是沒有意義的，因為有可能存在真實解梯度不是有限值情況。換句話說，假如真實解在空間\(H^1\)中，那么用\(W^{1,\infty}\)范數來估計誤差就是沒有意義的。
另一方面，我們總是選擇解空間\(Y\)的子空間\(Z \supset Y\)進行誤差衡量，這里子空間指的就是L2。對於有些情況來說，這是因為其物理意義決定的，L2范數在有些情況下具有一定的物理意義：對電磁場的積分\(\int E(x)^2\)可以表示電磁場包含的能量；同樣的，對波動方程積分可以表示其儲存的勢能大小。其他情況使用L2范數進行誤差估計只是因為它方便。在有限元方法中，L2范數可以通過\(U^TMU\)方便的計算，其中\(M\)為質量矩陣。
例如在非恆定熱傳導方程里，計算L2范數誤差就是錯誤的，因為其沒有任何物理依據，總能量及總物質量所在空間都是L1空間，在這個例子里，計算L2范數除了方便以外並沒有其他含義。^[4]

4.3.范數誤差計算表達式

以\(L_2\)空間范數為例，其原始表達式為

\[\begin{equation} \|f(x) \|_2= \frac{1}{A_{\Omega}} \sqrt{ \int_{\Omega}|f(x)|^2 \rm{dx} } \end{equation}\]

當我們使用FEM或DGM方法來計算\(L_2\)計算到的結果\(u_h\)，是原始函數\(u(x)\)的近似表達形式

\[\begin{equation} u_h = \sum_i u_i \varphi_i(x) \end{equation}\]

其中\(\varphi_i(x)\)是解空間內的基函數。因此，計算結果\(L_2\)誤差表達式為

\[\begin{equation} \| \Delta u \|_2 = \frac{1}{A_{\Omega}} \sqrt{ \int_{\Omega} \left|u(x) - u_h \right|^2 \rm{dx} } = \frac{1}{A_{\Omega}} \sqrt{ \int_{\Omega} \left|u(x) - \sum_i u_i \varphi_i(x) \right|^2 \rm{dx} } \end{equation}\]

實際上，當基函數為Lagrange函數時，更簡單的方法是令精確解也采用基函數進行表示，然后代入方程中進行計算

\[\begin{equation} \| \Delta u \|_2 = \frac{1}{A_{\Omega}} \sqrt{ \int_{\Omega} \left|u(x) - u_h \right|^2 \rm{dx} } = \frac{1}{A_{\Omega}} \sqrt{ \int_{\Omega} \left( \sum_i \left|u(x_i) - u_i \right|^2 \varphi_i(x) \right) \rm{dx} } \end{equation}\]

Reference

---

1. Correct way of computing norm L2 for a finite difference scheme ↩︎

2. Differences between the L1-norm and the L2-norm (Least Absolute Deviations and Least Squares) ↩︎

3. L1 norm and L2 norm ↩︎

4. What norm to choose when? ↩︎