非精確Newton 方法中線性迭代收斂判據研究

馮選燕，燕振國，朱華君，馬燕凱，馮新龍

(1.新疆大學 數學與系統科學學院，烏魯木齊 830000；2.空天飛行空氣動力科學與技術全國重點實驗室，綿陽 621000)

0 引 言

非線性問題廣泛存在于各類學科和實際問題中，是科學研究中的難點之一[1]。其高效數值求解方法的研究不僅在計算數學基礎理論上具有重要意義，而且在科學與工程計算中處于重要地位，能夠促進計算流體力學（Computational Fluid Dynamics, CFD）等科學研究和工程應用的發展。

具體到CFD 領域，其控制方程（Navier-Stokes equations，即N-S 方程組）具有高度非線性特征，對于絕大部分實際流動問題，往往只能對控制方程進行離散求解，為了獲得全離散方程，需要同時進行空間離散和時間離散。通常的空間離散包括限差分法、有限體積法和間斷有限元（discontinuous Galerkin-finite element method, DG-FEM）等數值方法[2-4]，時間離散則通常被劃分為顯式和隱式兩類[5]。其中隱式時間推進方法是解決CFD 剛性流動模擬效率問題的重要手段[6]，然而隱式時間離散后不可避免地需要求解非線性離散系統。發展這類大型非線性離散系統的基礎理論和高效求解方法是提高CFD 計算效率的關鍵手段[7]。

Newton 法是求解非線性方程的常用方法之一，其二階的收斂效率在高效CFD 計算中十分具有吸引力[8]。然而傳統的Newton 法需要計算非線性殘差方程的Jacobian（雅可比）矩陣來構造線性迭代方程組，并且要求對線性方程組進行精確求解，上述兩點在大規模CFD 計算中都是較難實現的。為了避免上述困難，通常會在上述兩個環節引入近似。例如，CFD 中常用的Jacobian-free Newton-Krylov（JFNK）方法就采用有限差分來近似Jacobian 矩陣，引入了Jacobian 矩陣計算誤差[9-11]；而線化后的方程組系統通常采用迭代方法進行求解，迭代方法會引入線性迭代誤差[12-13]。所以CFD 等大型模擬中實際上采用的是非精確Newton 方法[14]。

對于非精確Newton 方法，國內外開展了近似方法的系列研究工作。Vanden 和Orkwis[15]討論了精確和數值Jacobian 矩陣的求解過程及其對計算的影響，數值測試結果表明：在這兩種情況下獲得了較一致的迭代收斂次數及相近的收斂殘差，并指出對于較為復雜的非線性系統，數值Jacobian 矩陣可能更具有優勢。Ezertas[16]等在此基礎上進一步研究了數值Jacobian 矩陣計算中誤差的來源，探究并給出了使得有限差分近似誤差最小化的有限差分擾動幅值。測試發現：在該條件下，采用數值Jacobian 矩陣表現出了與精確Jacobian 矩陣相當的收斂性能。基于此，數值Jacobian 矩陣具有實現簡單等優勢，被廣泛應用于CFD 大規模數值模擬中。然而，Yang 等[17]的數值測試結果表明：采用無矩陣（Jacobian-free, JF）方法的數值Jacobian 矩陣并不能獲得與精確Jacobian 矩陣相當的收斂特性，有限差分誤差的引入會對迭代的收斂過程產生明顯的影響。

另外，非精確Newton 方法在每一步的內迭代中通過線性迭代收斂判據來控制求解線性方程組時解的精確程度，在此過程中會產生迭代的收斂誤差。對于該問題，文獻[18]將其稱為Newton 迭代的強制項，其中強制項控制系數 ηk的取值大小對收斂性有很大的影響。一些學者[19]結合理論分析和數值測試給出了該參數的不同選取方式，并給出了強制項能夠保證Newton 方法二階收斂速度的理論分析。然而，在考慮Jacobian 矩陣誤差的情況下，線性迭代收斂判據的選取方式尚未得到深入研究，需要進一步考慮Jacobian矩陣誤差的影響并展開對比研究。

因此，本文重點對存在Jacobian 矩陣誤差情況下不同線性迭代收斂判據的表現進行了研究，分析了不同線性迭代判據所引起的過度求解問題。結合兩種類型的判據發展了一種新的線性迭代收斂判據，并利用數值測試對新提出線性迭代收斂判據的優勢進行了驗證。

對于本文所關心的CFD 問題，其求解的Euler（歐拉）和N-S 方程，可以描述為以下形式：

1 控制方程和數值方法

1.1 隱式時間推進方法

對于Euler 或N-S 方程，式（3）是非線性方程組。采用Newton 方法對該方程進行線性化后迭代求解，迭代式為：

其中，J(Q)是 殘 差 向量N(Q)的 Jacobian 矩 陣，ΔQn=Qn+1 ?Qn。

為了增強計算的穩定性，通常會在式（4）左端添加偽時間項，此時不考慮時間項的計算精度，時間離散通常采用實現簡單、存儲量小的Euler 隱式格式，此時迭代方程變為：

式（5）中時間步長 Δt→∞時，退化為式（4）的形式。

1.2 Newton-Krylov 方法

Newton-Krylov 方法[20]目前是CFD 領域常用的大型非線性方程組求解算法，基于Krylov 空間，其具體形式為：

該類方法在Krylov 空間內構造正交的基向量，并在該空間中尋找使得線性方程組誤差最小的數值解，因此有殘差單調下降等優勢。其收斂速度依賴于線性方程組系數矩陣的條件數（由矩陣的特征值分布決定），當矩陣的條件數較小時，Krylov 方法可以在幾次搜索內將殘差降到很小的量級。這里為了求解Euler 和N-S 方程產生的非對稱系統，采用廣義最小殘差法（generalized minimum residual，GMRES）[21-22]，具體如下：

其中：M代表最大內迭代次數， σ為GMRES 迭代的收斂判據。Jqi的計算可采用有限差分近似得到，具體形式在本文1.3 節中詳細介紹。N代表GMRES 的最大重啟次數，重啟可以有效避免存儲過多Krylov 空間向量，減少正交化操作，是大規模計算中常用的技術。

而采用Krylov 方法求解線性方程組時帶來的線性迭代誤差，主要依賴于線性迭代收斂判據：

其 中， ΔQk=Qk+1?Qk； σ可控制線性迭 代 的 求 解 精度，當其取值過大時，較大的線性迭代誤差可能對Newton 迭代的收斂造成較大影響，而取值過小又容易導致過度求解，造成計算浪費，因此 σ的選擇對于發展高效的數值方法來說是十分重要的。

1.3 線性化方法及其誤差

在每一步Newton 迭代的過程中，需要采取直接法或迭代法求解由Jacobian 矩陣J組成的線性方程組。在標準的Newton 方法中需要采用精確方法計算J，然而對于大規模CFD 計算而言，計算精確的雅可比矩陣所產生的計算量和存儲量都是較大的。因此，歷史上提出了多種近似雅可比矩陣方法，其中一種是對雅可比矩陣的表達式進行近似，例如在經典的LUSGS 方法中對流通量和黏性通量的Jacobain 矩陣分別采用一階精度近似和黏性譜半徑近似[5]；另外一種方法則是采用數值方法計算Jacobian 矩陣，例如CFD中常用的JF 方法就是采用有限差分方法近似Jacobian 矩陣；此外，還存在上述兩種近似的混合方法[11]。本文主要討論采用JF 方法的誤差及其影響。

在構造Krylov 子空間（式6）過程中，令q=Jj?1r(k)，則Jq的計算可采用有限差分近似計算Jacobian 矩陣與向量的乘積：

2 誤差影響研究

本文考慮的非精確Newton 迭代法在求解大型非線性方程組時相對精確Newton 方法存在兩類誤差：Jacobian 矩陣誤差和線性迭代誤差。這里首先給出存在兩種誤差情況下非精確Newton 方法的迭代式，其次給出了兩種誤差對Newton 迭代式的影響，最后分析了Jacobian 矩陣誤差來源驗證及其影響。

2.1 誤差及其分析

非精確Newton 迭代的過程中，Jacobian 矩陣誤差和線性迭代誤差都對整個非線性方程組求解過程中的收斂性造成一定影響。對于經典Newton 法，解的迭代式為：

作為示例，這里給出Euler 方程采用DG 方法離散后的具體離散形式和Jacobian 矩陣的形式，詳細矩陣系數和N-S 方程的離散形式參照文獻[11]。具體殘差形式為：

這里d為空間維數，M=BTΛ(wJ)B為 質量矩陣， Λ為對角矩陣，D j為第j個方向的微分矩陣，B為從求解點或者模態系數到通量點的變換矩陣、BΓ為到單元邊界上的通量點的變換矩陣，w和J分別代表數值積分權重和網格變換Jacobian，Mc為 單元邊界基函數 ?Γ和基函數 ?之間的映射矩陣，F^n為對流數值通量。

其中，空間殘差項N(Q) 關 于守恒變量Q求導的Jacobian 矩陣J(Qn)形式：

?L(Qn)/?Q具 有 塊 結 構，例 如 ?Le1(Qn)/?Qe2代 表 第e1行和第e2列的塊Jacobian 矩陣元素，其含義為第e1個單元的 L 對第e2個單元的Q的Jacobian 矩陣塊。其具體形式為：

2.2 Jacobian 矩陣誤差來源驗證及其影響

本小節首先簡單驗證了采用有限差分近似的Jacobian 矩陣的相對誤差變化規律，驗證了Jacobian矩陣誤差主要來源于差分格式帶來的截斷誤差和計算機的舍入誤差[16]，而且表明在特定的數據精度和計算格式下，都存在一個使得誤差最小的ε，即最優ε，如圖1 所示。圖中FD 表示一階前向差分，CD 表示二階中心差分。

圖1 數值Jacobian 矩陣的相對誤差變化規律Fig.1 Relative error variations of the numerical Jacobian matrices

為進一步測試Jacobian 矩陣誤差對收斂速度的影響，這里利用周期計算域內的一維Euler 系統所搭建的簡單模型問題進行研究，此處選取120 個網格點，通過控制網格拉伸比使得網格向計算域兩側匯聚，并通過控制計算網格的大小比值來控制問題的剛性，此處網格拉伸比為1.05。空間離散采用一階迎風格式，該模型問題Jacobian 矩陣的條件數約為1 ×108。圖2 中，橫坐標（Non-Iter）代表非線性迭代次數，結果驗證了當數值Jacobian 矩陣近似誤差較大時會引起效率降低甚至無法收斂的問題[16-17]，然而當Jacobian矩陣誤差較小時則可能避免上述問題。

圖2 不同量級Jacobian 矩陣誤差對非精確Newton 迭代的影響Fig.2 Influence of the Jacobian matrix error order on the inexact Newton iteration

3 線性迭代收斂判據研究

對于Newton 類非線性系統求解方法，線性迭代誤差雖然不影響計算最終收斂后的精度，但是對非線性系統能否收斂以及收斂時的計算效率有較大的影響。其中，相關研究[18]給出的線性系統迭代收斂判據的選取準則要求，在保證Newton 方法二階收斂效率的前提下盡量避免線性系統的過度求解。所謂過度求解就是當線性系統的求解精度達到一定程度以后，繼續提升線性系統的求解精度將無法顯著提升非線性系統的收斂效率，此時額外增加的線性系統迭代步數就造成了線性系統的過度求解。從提高計算效率角度，避免過度求解顯然是必須考慮的問題。

常用的迭代誤差收斂判據有兩種，一種是絕對量型收斂判據，即選取某個絕對的數值作為迭代收斂的判據，為了使得該策略中的參數更具有普適性，本文采用如下的形式：

另外一種是相對型收斂判據，即將線性系統的迭代收斂判據與非線性系統的殘差進行關聯，其形式如下：

此處 ηk為線性系統收斂判據的主要控制參數，其含義為線性系統迭代誤差要相對當前的非線性系統殘差小 ηk倍以上。文獻[19]對式（20）的線性迭代收斂判據進行了較系統的研究，并從能夠保持Newton 方法二階收斂效率的角度給出了相應的數學證明和數值測試；然而并未對式（19）的絕對量型收斂判據和存在Jacobian 矩陣誤差情況下式（20）的迭代收斂誤差的影響開展研究。

3.1 不同 ηk的對比研究

這里采用風雷高精度軟件平臺[24]開展測試研究。控制方程采用N-S 方程，算例為黏性NACA 0012翼型流動問題，來流Ma=0.1，雷諾數取100，攻角為0。，翼型表面采用無滑移邊界條件，計算區域采用20 倍翼型長度，在外邊界采用遠場邊界條件。采用2842個四邊形網格對計算域進行離散，如圖3 所示。在當前的定常計算中，添加了偽時間項來控制線性系統的剛性以增強計算的穩定性，計算結果如圖4 所示。本計算主要關注殘差的收斂過程。

圖3 NACA 0012 算例的計算網格分布Fig.3 Computational grid distribution for the case of NACA 0012

圖4 NACA 0012 算例計算結果密度分布Fig.4 Density distribution for the flow around NACA 0012

如圖5 所示，首先對比測試了 ηk分別取1×10?3、1×10?4和1×10?6情況下采用式（8）的數值Jacobian 矩陣（JF1）的迭代收斂情況，作為對比，本文同時給出了采用精確Jacobian 的結果（JF0）。從線性迭代步數（Iter）上來講，顯然采用更小的 ηk所需的線性迭代步數顯著增加，線性方程組的求解精度相應提升。然而除了 ηk=1×10?3以外，所有非線性系統殘差基本完全重疊；當 ηk=1×10?3時，收斂曲線在收斂初段和最終階段已經出現了較明顯的偏差。上述結果表明，當ηk≤1 ×10?4時線性迭代已經足夠精確，繼續降低 ηk則會造成過度求解問題，使計算效率降低。

圖5 NACA 0012 算例在不同Jacobian 矩陣下非線性迭代殘差和GMRES 迭代步數對比Fig.5 Comparison of the nonlinear iterative residual and GMRES iteration steps for different Jacobian matrices in the case of NACA 0012

另外需要強調的是，在迭代接近收斂時，基于ηk的迭代收斂判據都會造成過度求解問題。本文進一步通過增加GMRES 的迭代重啟次數來增加允許的最大線性系統迭代總步數，從而進一步研究造成上述問題的原因，結果如圖6 所示。在迭代的最終收斂階段，無論將允許的最大線性系統迭代步數增大到多大（在本文的測試范圍內），線性方程組均無法收斂到 ηk=1×10?6，相應的線性方程迭代的總步數達到允許的最大值。與之相反，采用精確Jacobian 矩陣的計算結果，雖然在收斂的最終階段迭代步數同樣增加，但是并不存在上述現象，說明上述問題是當所要求的線性系統迭代誤差太小時，Jacobian 矩陣誤差會造成線性系統無法收斂引起的。然而從非線性系統的殘差角度來看，各條曲線是基本重合的，進一步說明采用近似Jacobian 所增加的線性系統迭代是由于收斂判據過嚴引起的過度求解，需要發展更加合理的迭代收斂判據避免上述問題。

圖6 不同GMRES 允許的最大迭代步數下的線性迭代次數Fig.6 Linear iterative number under different allowable maximum GMRES iteration steps

3.2 不同 ζk的對比研究

類似3.1 節，這里對式（19）的收斂判據進行測試，結果如圖7 所示。本文分別給出了 ζk取1×10?6、1×10?10、1×10?14的測試結果，為了對比，圖中同時給出了 ηk=1×10?6情況下的收斂曲線。可以看到，采用ζk類 型的迭代收斂判據，只有在 ζk取很小的值時才能保證迭代收斂。但是，當 ζk取值很小時（如1×10?14），在迭代的初始階段會造成較嚴重的過度求解現象。而在迭代接近收斂的階段，采用 ζk類型的迭代收斂判據卻可以較好地避免過度求解的問題。

圖7 NACA 0012 算例中不同ζ 取值情況下非線性迭代殘差和GMRES 迭代步數對比Fig.7 Comparison of the nonlinear iterative residual and GMRES iteration steps for different ζ in the case of NACA 0012

3.3 一種新的線性迭代收斂判據

基于以上研究結果，本文發展了一種新的迭代收斂判據：

其中 max函數內兩項分別為式（19）和式（20）對應的絕對收斂判據和相對收斂判據，這里取兩項中的較大值。 在非線性迭代初期||N(Qk)||2較大，因此ηk||N(Qk)||2的量級會大于ζk||Q0||2，從而避免了后者在非線性迭代初期所引起的嚴重過度求解問題；隨著迭代逐漸接近收斂， ||N(Qk)||2逐漸減小，只要取合理的參數，可以使得 max函 數在接近收斂時取 ζk||Q0||2，此時可以有效避免圖5 中的過度求解問題。在此基礎上，本文還在 max函 數外添加了相應的 min函數，其作用是保證每個線性迭代至少要下降一定量級（本文中取ηmin=0.01），從而避免類似圖7 中迭代不能完全收斂的問題。本文中 ζk取1×10?13，主要從計算的舍入誤差角度進行選取，由于最終線性迭代誤差會加到Qk上， ζk=1×10?13的參數選擇可以避免線性方程組收斂到比Newton 迭代步本身舍入誤差更小的量級，從而避免類似圖5 中的過度求解問題。

圖8 同樣采用NACA 0012 算例對不同迭代收斂判據進行了對比，直至迭代第19 步采用式（21）的新迭代收斂判據（NewTol）得到了與 ηk=1×10?6完全相同的計算結果。此后采用 ηk=1×10?6的線性迭代步數迅速增加到允許的最大值，線性迭代無法收斂，而采用式（21）的線性迭代步數則開始下降，說明如設計的那樣，式（21）中 ζk||Q0||2部分開始發揮作用從而大幅減少了線性迭代的計算量。直至第23 步（殘差收斂至1×10?13）時采用新判據的非線性系統殘差依然與ηk=1×10?6重疊，隨后采用式（21）的線性迭代步數開始顯著減小，非線性殘差曲線不再重疊。表1 給出了殘差收斂到不同量級時所需要的線性迭代總步數，該指標能夠很好地反映出所消耗的總計算量。當殘差收斂到1×10?10以下時，JF0- ηk=1×10?6（精確Jacobian、ηk=1×10?6）、JF1- ηk=1×10?6（近似Jacobian、 ηk=1×10?6）和JF1-NewTol（近似Jacobian、新判據）所需線性迭代步數分別為13113、19670 和12692 步；當殘差收斂到1×10?13以下時，所需線性迭代步數分別為16964、29670和13870 步；當殘差收斂到1×10?16以下時，所需線性迭代步數分別為21473、44670 和15453 步。可以看到，新發展的線性迭代收斂判據能夠有效避免過度求解的問題，顯著降低所需要的線性系統迭代步數。線性迭代步數最多降低3 倍，從而提高計算效率。

表1 不同判據情況下迭代收斂所需線性迭代步數對比Table 1 Comparison of the linear iteration steps under different convergence criteria

圖8 新線性迭代收斂判據的殘差和迭代步數Fig.8 Nonlinear iterative residual and GMRES iteration steps for the new linear iteration convergence criterion

為了進一步驗證上述迭代收斂策略在不同計算狀態中的通用性，本文通過改變離散格式的計算精度來改變非線性系統的特性和規模。將DG 格式的計算精度提高到三階，此時總共求解的網格自由度變為17052 個。圖9 給出了三階精度的計算結果，可以看到，采用新的線性迭代收斂判據同樣能夠很好地避免接近收斂狀態下的過度求解問題，從而顯著提高計算效率。

圖9 三階DG 格式計算中新線性迭代收斂判據的殘差和迭代步數Fig.9 Nonlinear iterative residual and GMRES iteration steps for the new linear iteration convergence criterion with third-order DG scheme

4 結 論

針對CFD 中常見的基于Newton 類方法的非線性系統求解問題，本文通過數值測試，對存在Jacobian矩陣誤差情況下每個Newton 迭代步中線性迭代收斂判據進行了研究，最終結合兩種類型的迭代收斂判據，發展了一種新的判據。得到的主要結論如下：

1）存在Jacobian 矩陣誤差情況下，線性迭代收斂判據的影響會更加顯著；

2）基于絕對量的線性迭代收斂判據容易在非線性迭代初期引起過度求解問題；相對非線性殘差的線性迭代收斂判據容易在接近最終收斂時引起過度求解問題；

3）新型線性迭代收斂判據能夠有效地降低過度求解問題，顯著降低收斂所需要的線性迭代步數，從而提高計算效率。

下一步，我們將進一步對Jacobian 矩陣誤差對非線性系統迭代收斂效率的影響開展研究，并將發展的線性迭代收斂誤差應用到更加復雜的流動模擬中。