時空多項式配點法求解三維Burgers方程*

曹艷華，張姊同，李 楠

（華東交通大學 理學院，南昌 330013）

引 言

Burgers 方程可以用來模擬沖擊波的傳播和反射，是一個被廣泛應用于流體力學、非線性聲學、氣體動力學等領域的非線性偏微分方程[1-3].由于Burgers 方程是從Navier-Stokes 方程簡化而來，因此對于該方程的研究眾多，其中不乏非常有意義的研究成果.Hopf[2]和Cole[4]提出了著名的Hopf-Cole 變換來研究Burgers 方程的基本性質.他們發現Reynolds 數對于平衡方程中的非線性對流項和黏性擴散項起著舉足輕重的地位.當Reynolds 數變大的時候非線性對流項占主導地位，由此導致了Burgers 方程的強非線性性，使得其求解異常困難.盡管通過Hopf-Cole 變換可以求得一些特殊類型Burgers 方程的精確解，但對于絕大多數Burgers 方程而言，其精確解是難以得到的.學生在學習Burgers 方程的強非線性理論時，不易理解其含義[5-7].

Gulsu 使用有限差分法求解具有小Reynolds 數的Burgers 方程[8]，對于具有大Reynolds 數的Burgers 方程，有限差分法的求解精度很差[9].Turgut 等采用半步長方法[10]，即時間方向采用半隱式有限差分格式、空間方向采用漸進展開格式求解具有大Reynolds 數的Burgers 方程.除了采用傳統的網格類方法，近些年來，無網格類方法也被廣泛用來求解Burgers 方程[11-17].

所謂時空類方法，即在求解發展方程時將時間變量t視為普通的空間變量x，其中著名的時空類方法有時空有限元法[18]、間斷Galerkin 有限元法[19]以及Mohammed 等提出的局部時空徑向基函數等方法[20].

2017年，Dangal 等提出了一種新的多項式特解法用于求解偏微分方程[21].2020年，Cao 等在多項式特解方法的基礎上，提出了一種時空多項式特解法求解一維和二維的發展方程[22-23].受其啟發，本文采用時空多項式特解法求解三維具有大Reynolds 數的Burgers 方程.相比于其他的無網格類方法，多項式特解配點法只需知道配點的位置信息，別的任何未知參數都不需要確定，其簡便性遠遠超過現有的其他時空類無網格方法，例如，使用時空徑向基函數方法時，不僅需要確定徑向基函數的形狀參數，還需要確定配點的中心位置、個數以及插值半徑，這些參數的確定沒有一般規律可循.數值例子表明，隨著Reynolds 數的增加，時空多項式特解法不僅可以保持數值解的高精確性，其穩定性也沒有受到任何影響.

注1所有以多項式作為基函數的方法中，由于系數矩陣隨著多項式階的增加，其條件數迅速增大，很快數值溢出，因此本文的所有算例均采用了多尺度技巧以降低系數矩陣的條件數，得到穩定的數值解[22].

注2多尺度技巧可以有效降低多項式特解中系數矩陣的條件數，但對于多項式直接作為基函數得到的系數矩陣，多尺度技巧不能有效降低系數矩陣的條件數，這也是多項式函數直接作為基函數求解三維Burgers 方程得到的近似解精度不高的原因之一.

1 時空多項式特解法

考慮如下的三維偏微分方程初邊值問題：

其中，Ω為三維空間中一有界連通的規則或不規則區域，f(x,t),g(x,t)和h0(x)為給定函數，?t=?/?t，L為線性或非線性偏微分算子，B為邊界算子，(0,τ]為待求時間區間.

在時空類方法中，需要將三維待求解問題轉化為相應的四維問題，新的“空間變量”為“t+(x,y,z)”，將?×{t=0} 和 ??×(0,τ)作為新四維問題的邊界.在時空類方法中，不再區分時間變量t和 空間變量x=(x,y,z)，因此，相應于方程(1)的新四維問題為

在新四維空間中，d階完全多項式基函數為

其中 Λ={0≤i+j+k+l≤d}.

空間 Qd中，線性無關的基函數個數為(d+1)(d+2)(d+3)(d+4)/24.假設方程(2)的解h(x,t)可以寫為特解 χijkl(x,t)的一個線性組合：

其中 ζijkl為待求解系數，特解 χijkl(x,t)滿足

L1為L 中的空間常系數線性算子部分.在四維空間中，關于二階線性常系數偏微分方程的特解 χijkl(x,t)有以下引理.

引理1[21]令(ω1,ω2,ω3,ω4)=(x,y,z,t)，則對于下列偏微分方程

為式(6)的一個多項式特解，其中d=i+j+k+l,L=

在區域 ? 內部任取ni個配點，邊界 ? ?×(0,τ)上任取nd1個配點，邊界? ×{t=0} 上任取nd2個配點，配點總數目為N=ni+nd1+nd2，則方程(2)的配點型強格式為

迭代求解代數方程組(8)，可得到待定系數{ζijkl|Λ}.將{ζijkl|Λ}代入式(4)即得到方程(2)的近似解.

為了比較說明多項式特解法的優點，假設方程(2)的解h(x,t) 可表示為多項式基函數?ijkl(x,t)=xiyjzktl的一個線性組合：

其中 ?ijkl為待定系數.將(x,t) 代入方程(2)，且在區域 ? 內部任取ni個配點，邊界 ? ?×(0,τ)上任取nd1個配點，邊界?×{t=0}上任取nd2個配點，配點總數目為N=ni+nd1+nd2，則方程(2)的配點型強格式為

迭代求解代數方程組(10)，得到待定系數{?ijkl|Λ}.將{?ijkl|Λ}代入式(9)即得到方程(2)的近似解.

注3為使代數方程組(8)和(10)可解，配點的數目N 需大于完全多項式空間 Qd的維數(d+1)(d+2)(d+3)(d+4)/24.

本文的計算區域為規則區域和不規則區域，所用配點有兩種方式：一是內點和邊界點均為規則配點，二是內點和邊界點均為隨機配點.以二維區域為例，圖1(a)、1(b)所示的是規則配點，圖 1(c)、1(d)所示的是隨機配點.二維非規則空間區域在新的三維時空域中隨機取點如圖2所示，規則二維空間區域在新的三維時空域中隨機取點類似.

圖1 配點圖：(a) 正方形區域規則取點；(b) 星形區域規則取點；(c) 正方形區域隨機取點；(d) 星形區域隨機取點(籃圈：內點；紅星：邊界點)Fig.1 The collocation point diagram: (a) the regular domain with regular points; (b) the irregular domain with regular points; (c) the regular domain with random points; (d) the irregular domain with random points (blue circles: interior collocation points; red stars: boundary collocation points)

圖2 時空區域中星形區域隨機取點（籃圈：內點；紅星：邊界點）Fig.2 An irregular domain with random points in the space-time domain (blue circles: interior collocation points; red stars: boundary collocation points)

2 數值實驗

數值解與精確解之間的平方根誤差定義為

考慮下列三維Burgers 方程的初邊值問題：

其中，μ=1/Re＞0為黏性系數，Re為Reynolds 數.方程(11)的精確解為

f(x,t),g(x,t)和h0(x)為基于精確解(12)的函數.

本文的研究區域為正方體區域和bumpy-shaped 型非規則區域，如圖3所示.

圖3 空間求解區域：(a) 正方體區域；(b) bumpy-shaped 區域Fig.3 The spatial solution domains: (a) the cubic domain; (b) the bumpy-shaped domain

在數值實驗中，我們將比較使用不同時刻 τ以及不同Reynolds數所得的數值模擬結果.如無特別說明，ni表示區域內部的配點數，nd表示邊界??×(0,τ)∪?×{t=0} 上的配點數，nt表示測試點數.

2.1 大τ 情形

本小節的算例采用時空多項式配點法和直接使用多項式基函數時空配點法求解Reynolds 數Re=1時的方程(11).

表1給出了 τ=1，配點ni=3 819,nd=5 220,nt=1 919時，兩種區域對不同階多項式特解的計算時間和平方根誤差.可以看出，隨著特解階的增加，平方根誤差迅速下降；特解的階越高，計算所需時間越長.同時，可以看出，立方體區域上的平方根誤差下降得比bumpy 型區域上的平方根誤差快.

表1 ni=3 819,nd=5 220,nt=1 919,τ=1時，兩種區域多項式特解的平方根誤差Table 1 The RMSE with ni=3 819,nd=5 220,nt=1 919,τ=1 for cubic and bumpy-shaped domains

表2給出了 τ=1，配點ni=3 819,nd=5 220,nt=1 919時，兩種區域對不同階多項式基函數的計算時間和平方根誤差.可以看出，隨著多項式階的增加，平方根誤差未有明顯下降；多項式基函數的階越高，計算所需時間越長.通過分析多項式基函數的系數矩陣，可以看到多尺度方法對其系數矩陣的條件數未有明顯改進，這也是多項式基函數的平方根誤差并不隨著多項式基函數階的增加而降低的原因之一.

表2 ni=3 819,nd=5 220,nt=1 919,τ=1時，兩種區域多項式基函數的平方根誤差Table 2 The RMSE with ni=3 819,nd=5 220,nt=1 919,τ=1 for cubic and bumpy-shaped domains

τ=2的數值結果如圖4所示.易見，對這兩種區域，當特解基函數的階小于或等于12 時，所得計算結果的誤差都以指數量級迅速下降.此后，再增加特解基函數的階，所得計算結果的精度不再有顯著提高，但也不會降低.當 τ=100時，兩種區域的計算結果如圖5所示.顯然，當特解基函數的階小于或等于12 時，所得計算結果都快速收斂.

圖4 τ=2,Re=1時，兩種區域的計算結果Fig.4 Results for τ=2 andRe=1

圖5 τ=100,Re=1時，兩種區域的計算結果Fig.5 Results for τ=100 andRe=1

為進一步分析所得結果，τ=2時，正方體求解區域上的系數矩陣的條件數如圖 6(a)所示.顯然，在多項式特解中，使用多尺度技巧與否，所得矩陣的條件數相差巨大.如果不使用多尺度技巧，其條件數隨著多項式特解階的增加迅速增大，很快便數值溢出：多項式特解的階等于12 時，所得數值結果發散.多項式特解的階等于20 時，條件數溢出.使用多尺度技巧，可以有效控制條件數的增長，獲得穩定的解.因此，在多項式類的解法中，多尺度技巧的使用是至關重要的.Bumpy-shaped 區域上條件數的情形與正方體求解區域上的情形類似.對于多項式基函數而言，系數矩陣的條件數使用和不使用多尺度技巧差別不大，如圖 6(b)所示，這或許是直接使用多項式基函數所得的近似解精度不高的原因之一.

圖6 τ=2時，正方體區域上系數矩陣的條件數：(a) 多項式特解系數矩陣的條件數；(b) 多項式基函數系數矩陣的條件數Fig.6 Matrix condition numbers of the cubic domain for τ=2: (a) the condition number of polynomial particular solutions; (b) the condition number of polynomial basis functions

根據本小節數值例子，可得下列結論.

1） 對于規則或不規則的計算區域，采用不同數目的配點，計算結果會有所不同，但是這些結果之間的差別很小，不足以影響到計算精度.

2） 采用不同數目配點的最大的區別在于計算時間，配點越多，所需的計算時間也越長.以正方體區域為例，ni=6 500，nd=6 000，nt=256時 所需的計算時間約為8 930 s.如果ni增加到10 000，nd增加到8 128，nt保持不變，則所需的計算時間增加4 529 s.由此可見，配點數目嚴重影響計算時間.在保證近似解精度的情況下，應選用適當數目的配點.對于非規則區域，其情形與正方體區域的情形類似，此處不再贅述.

3） 本算法中的多項式特解基函數 χijkl需要提前生成.生成1 階到30 階的多項式特解，所需總時間約為10 415 s.值得慶幸的是，特解基函數 χijkl一旦生成便可重復使用，甚至對于不同的微分方程，只要其中的線性微分算子項相同，也可使用.

2.2 不同Reynolds 數的情形

本小節中所有算例的τ均取1.在求解Burgers 方程時，的取值對所求數值解的精度影響巨大.現有研究結果表明，Reynolds 數Re越大，所得數值解的精度越低，求解難度越大.

當 μ=1，所得結果如圖7所示.正方體區域上數值解的精度達到10-11，即使在非規則區域bumpy-shaped上，數值解的精度仍然能達到10-9.當μ=0.5時，所得結果如圖8所示.正方體區域上數值解的精度達到10-8，非規則區域bumpy-shaped 上數值解的精度約為1 0-6.繼續增大Reynolds 數，μ=0.1，正方體區域上數值解的精度達到10-5，而非規則區域bumpy-shaped 上數值解的精度約為1 0-3，且數值解不再像 μ=1時那樣光滑.因此，本文數值結果進一步驗證了Reynolds 數越大，所得數值解的精度越低，求解難度越大，這與以往的研究結果一致.

圖7 τ=1,μ=1時，兩種區域的計算結果Fig.7 Results for τ=1 andμ=1

圖8 τ=1,μ=0.5時，兩種區域的計算結果Fig.8 Results for τ=1 andμ=0.5

表3給出了 τ=1，配點ni=3 819,nd=5 220,nt=1 919,μ=0.5時，兩種區域直接使用多項式基函數的計算時間和平方根誤差.可以看出，隨著多項式基函數階的增加，平方根誤差幾乎沒有降低.

表3 ni=3 819,nd=5 220,nt=1 919,μ=0.5時，兩種區域多項式基函數的平方根誤差和計算所需時間Table 3 The RMSE with ni=3 819,nd=5 220,nt=1 919,μ=0.5 for cubic and bumpy-shaped domains

3 結 論

三維Burgers 方程因其維數高，大Reynolds 數時的強非線性性，現有研究結果不多.本文提出了一種新型時空多項式特解配點法求解三維Burgers 方程，不同于以往的多項式函數直接做基函數的方法，多項式特解做基函數的優點是對控制方程中的線性導數項不進行求導，因此多項式特解的階比較高時，所得系數矩陣中元素的量級差別不大，從而算法的穩定性強，精度很高.采用配點法時，時空方法的維數比初始問題的維數增加了一維，所需配點數目較多，下一步的工作可以通過采取有效的降階技巧減少對配點的要求，將系數矩陣變為稀疏矩陣，以期得到更好的模擬結果.通過本文的工作，使學生對高維強非線性的Burgers 方程的數值求解有更加深入的理解.