應用重心插值的Burgers方程數值解法

摘要： 利用重心插值構造求解Burgers方程的無網格數值方法。在時間方向用Crank-Nicolson差分法對方程進行離散，在空間方向用重心插值切比雪夫配點法逼近函數本身及其空間導數。 對全離散數值格式進行相容性分析。數值算例表明：與經典的有限差分方法比較，重心插值配點法用較少的節點就能達到較高的精度。

關鍵詞： Burgers方程；重心插值配點法；Crank-Nicolson差分法；相容性分析

中圖分類號： O 241.82文獻標志碼： A 文章編號： 1000-5013（2025）01-0104-09

Numerical Solution Method of Burgers Equation Using Barycentric Interpolation

TENG Yuhang，LAI Yiying，HUANG Langyang

（School of Mathematical Science，Huaqiao University，Quanzhou 362021，China）

Abstract： The meshless numerical method of Burgers equation is solved using barycentric interpolation construction. The equation is discretized using the Crank-Nicolson difference method in the time direction. The function is approximated to itself and so is its spatial derivative using barycentric interpolation Chebyshev collocation method in spatial direction. The compatibility analysis of the fully discrete numerical value scheme is perfomed. Numerical experiments show that，compared with the classical finite difference method，the barycentric interpolation collocation method can achieve higher accuracy with fewer nodes.

Keywords： Burgers equation；barycentric interpolation collocation method；Crank-Nicolson difference method；compatibility analysis

Burgers方程是一類描述對流與擴散相互作用的非線性偏微分方程^［1］，是Navier-Stokes方程中忽略壓力項后的簡化形式。考慮如下Burgers方程，即

式（1）中：Ω R^d（d=1，2）為有界凸區域，邊界分段光滑。

Burgers方程在流體力學、氣體動力學、交通流等眾多領域都有著重要的作用和地位。許多學者對該問題做了大量研究，并提出了數值解法。陳蓮^［2］建立兩種求解Burgers方程的Crank-Nicolson差分法格式，并對其進行誤差分析。Wang等^［3］建立一種求解粘性Burgers方程的能量守恒型四階隱式緊差分格式。Zhang等^［4］利用線性化緊致差分格式求解具有Burgers型非線性項的二維Sobolev方程。Xu等^［5］建立兩種求解非線性Burgers方程的修正的非協調有限元全離散格式，并對其進行超收斂分析。Wang等^［6］建立一種求解Burgers方程的多區域Galerkin方法，并分析其收斂性及穩定性。Zhao等^［7］利用時空連續Galerkin方法對二維Burgers方程進行數值求解，并給出先驗誤差估計。Wang等^［8］利用弱伽遼金有限元方法求解一類時間分數階的廣義Burgers方程。

以上求解Burgers方程的數值方法都是基于網格剖分建立離散格式，Chebyshev節點的重心Lagrange插值公式可以有效克服Runge現象。很多學者將其推廣到求解各類微分方程與積分方程，如Volterra積分方程^［9］、Allen-Cahn方程^［10-13］、Black-Scholes方程^［14］、分數階微分方程^［15-16］等。?mer^［17］應用重心插值配點法數值求解了二維和三維Klein-Gordon-Schr?dinger（KGS）方程。基于此，本文對應用重心插值的Burgers方程數值解法進行研究。

1 重心Lagrange插值

1.1 重心Lagrange插值函數

設插值節點為x_j，對應的實數為y_j（j=0，1，…，n），使用多項式插值，則在次數不超過n的多項式空間中，存在唯一的重心Lagrange插值函數p（x），滿足p（x_j）=y_j，j=0，1，…，n。重心Lagrange插值函數為

p（x）=∑n/j=0w_j/（x-x_j）y_j/∑n/i=0w_i/（x-x_i）=∑n/j=0D_j（x）y_j。（2）

式（2）中：D_j（x）為重心Lagrange插值基函數；w_j為重心Lagrange插值函數的重心權，w_j為

w_j=∏n/i=0，i≠j（x_j-x_i）^-1。（3）

由于重心Lagrange插值函數的重心權只與插值節點分布有關，在利用重心插值逼近未知函數時，采用第二類 Chebyshev節點，即

x_j=cosj/nπ， j=0，1，…，n。（4）

可以得到重心Lagrange插值函數的重心權為

w_j=（-1）δ_j，δ_j=1/2， j=0或n，1， 其他。（5）

對式（2）進行求導，得到函數p（x）在插值節點x_i處的導數，p^s（x_i）為

p^s（x_i）=d^sp（x_i）/dx^s=∑n/j=0D^（s）_j（x_i）y_j=∑n/j=0D^（s）_i，jy_j。（6）

將其寫成矩陣形式即為

p^（s）=D^（s）p。（7）

式（7）中：D^（s）為關于節點x_i（i=0，1，…，n）的s階微分矩陣，D^（s）=（D^（s）_i，j）_n×n，s=1，2。

通過對重心Lagrange插值基函數求導，分別得到一、二階微分矩陣的元素^［11］為

D^（1）_i，j=ω_j/ω_i/x_i-x_j， i≠j，D^（1）_i，i=-∑n/j=0，j≠iD^（1）_i，j，D^（2）_i，j=2（D^（1）_i，iD^（1）_i，j-D^（1）_i，j/x_i-x_j）， i≠j，D^（2）_i，i=-∑n/j=0，j≠iD^（2）_i，j。（8）

1.2 逼近性質

設p（x，y）是u（x，y）的重心Lagrange插值函數，滿足p（x_m，y_n）=u（x_m，y_n），則定義誤差函數為

e（x，y）=u（x，y）-p（x，y）。（9）

引理1^［18］ 如果u（x，y）∈C^m+1（［a，b］×［c，d］），其中，m=max{m，n}，則

|e（x，y）|≤u^（m+1）_∞（C₁el_x/2m^m+C₂el_y/2nⁿ）。（10）

式（10）中：e為自然對數；l_x代表區間［a，b］的一半；l_y代表區間［c，d］的一半。

同理可以得到

e_x（x，y）≤C^*₁‖?^（m+1）_xu‖_∞el_x/2（m-1）^m-1+C₂‖?^（n+1）_yu‖_∞el_y/2nⁿ，e_x，x（x，y）≤C^**₁‖?^（m+1）_xu‖_∞el_x/2（m-2）^m-2+C₂‖?^（n+1）_yu‖_∞el_y/2nⁿ，e_y（x，y）≤C₁‖?^（m+1）_xu‖_∞el_x/2m^m+C^*₂‖?^（n+1）_yu‖_∞el_y/2^（n-1）n-1，e_y，y（x，y）≤C₁‖?^（m+1）_xu‖_∞el_x/2m^m+C^**₂‖?^（n+1）_yu‖_∞el_y/2^（n-2）n-2。（11）

式（11）中：C₁，C^*₁，C^**₁，C₂，C^*₂，C^**₂為常數。

2 二維Burgers方程的離散格式

考慮二維Burgers方程，對時間區域進行K等分，記步長為τ=T/K，T為時間計算的終點。空間區域采用第二類Chebyshev節點分別離散為m+1個計算節點a=x₀lt;x₁lt;…lt;x_m=b和n+1個計算節點 c=y₀lt;y₁lt;…lt;y_n=d，得到計算節點（x_i，y_j，t_k），0≤i≤m，0≤j≤n，0≤k≤K。

2.1 半離散格式及相容性分析

記函數u（x，y，t）在節點x₀，x₁，…，x_m上的值為u（x_i，y，t）=u_i（y，t），則函數u（x，y，t）在節點x₀，x₁，…，x_m上的重心插值函數為

u（x，y，t）=∑m/s=0L_s（x）u_s（y，t）。（12）

式（12）中：L_s（x）表示x方向的重心插值基函數。

將式（12）代入方程（1）中，并令方程在點x₀，x₁，…，x_m上精確成立，得到常微分方程組，即

∑m/s=0L_sx_i?u_s（y，t）/?t+/∑m/s=0L_sx_iu_s（y，t）∑m/s=0L′_sx_iu_s（y，t）+∑m/s=0L_sx_i?²u_s（y，t）/?y²=/

v∑m/s=0L″_sx_iu_s（y，t）+∑m/s=0L_sx_i?²u_s（y，t）/?y²。（13）

式（13）中：L′_s（x_i）=C^（1）_s，i為關于節點x₀，x₁，…，x_m的一階微分矩陣C^（1）中的元素，類似地，L″_s（x_i）=C^（2）_s，i為二階微分矩陣C^（2）中的元素。

記u_i，j（t）=u_i（y_j，t），函數u_i（y，t）在節點y₀，y₁，…，y_n的重心插值函數為

u_i（y，t）=∑n/r=0T_r（y）u_i，r（t）。（14）

式（14）中：T_r（y）表示y方向的重心插值基函數。

將式（14）代入式（13），微分矩陣形式為

dU（t）/dt+diag（U（t））（C^（1）I_n+1+I_m+1D^（1））U（t）=v（C^（2）I_n+1+I_m+1D^（2））U（t）。（15）

式（15）中：

U（t）=［U₀（t）…U_m（t）］^T；

U_i（t）=［u_i，0（t）…u_i，n（t）］^T（i=0，…，m）；

D^（k）為關于節點y₀，y₁，…，y_n的k階微分矩陣。

為了方便進行相容性分析，算子定義為

Γu（x，y，t）：=u_t+u（u_x+u_y）-v（u_x，x+u_y，y）。（16）

令u（x_m，y_n，t）為u（x，y，t）的數值解，則可以得到

Γu（x_m，y_n，t）=0。（17）

定理1 若u（x，y，t）∈C^（m+1）（［a，b］×［c，d］）×C^（0）（0，T］，其中，m=max{m，n}，假設u（x_m，y_n，t）有界，則誤差估計為

u（x，y，t）-u（x_m，y_n，t）≤C^**₁?^（m+1）_xu_∞el_x/2（m-2）^m-2+/

C^**₂‖?^（n+1）_yu‖_∞el_y/2（n-2）^n-2。（18）

證明：Γu（x，y，t）-Γu（x_m，y_n，t）=u_t（x，y，t）-u_t（x_m，y_n，t）+/

u（x，y，t）u_x（x，y，t）-u（x_m，y_n，t）u_x（x_m，y_n，t）+/

u（x，y，t）u_y（x，y，t）-u（x_m，y_n，t）u_y（x_m，y_n，t）+/

v（u_xx（x_m，y_n，t）-（x，y，t））+v（u_yy（x_m，y_n，t）-u_yy（x，y，t））=/R₁+R₂+R₃+R₄+R₅。（19）

式（19）中：

R₁=u_t（x，y，t）-u_t（x_m，y_n，t）=u_t（x，y，t）-u_t（x_m，y，t）+

u_t（x_m，y，t）-u_t（x_m，y_n，t）=e_t（x_m，y，t）+e_t（x_m，y_n，t）。（20）

由引理1，有R₁=/e_t（x_m，y，t）+e_t（x_m，y_n，t）≤/C₁‖?^（m+1）_xu‖_∞el_x/2m^m+C₂‖?^（n+1）_yu‖_∞el_y/2nⁿ。（21）

R₂的估計為

R₂=/u（x，y，t）u_x（x，y，t）-u（x_m，y_n，t）u_x（x_m，y_n，t）≤/u_x（x，y，t）u（x，y，t）-u（x_m，y_n，t）+/

u（x_m，y_n，t）u_x（x，y，t）-u_x（x_m，y_n，t）≤/

C^*₁‖?^（m+1）_xu‖_∞el_x/2（m-1）^m-1+C₂?^（n+1）_yu_∞el_y/2nⁿ。（22） 同理有

R₃=u（x，y，t）u_y（x，y，t）-u（x_m，y_n，t）u_y（x_m，y_n，t）≤

C₁‖?^（m+1）_xu‖_∞el_x/2m^m+C^*₂‖?^（n+1）_yu‖_∞el_y/2（n-1）^n-1。 （23）

類似地，對于R₄，R₅的估計分別為

R₄=ve_x，x（x_m，y，t）+ve_x，x（x_m，y_n，t）C≤

C^**₁‖?^（m+1）_xu‖_∞（el_x/2（m-2））^m-2+C₂‖?^（n+1）_yu‖_∞-el_y/2nⁿ。（24）

R₅=ve_y，y（x_m，y，t）+ve_y，y（x_m，y_n，t）≤

C₁‖?^（m+1）_xu‖_∞el_x/2m^m+C^**₂‖?^（n+1）_yu‖_∞el_y/2（n-2）^n-2。 （25）

將式（21）～（25）代入式（19），則定理得證。

2.2 全離散格式及相容性分析

定理2 設u（x，y，t）∈C^（m+1）（［a，b］×［c，d］）×C^（2）（0，T］，假設u（x_m，y_n，t）有界，使用重心插值配點法求解的數值解為u^h（x_i，y_j，t_k），則如下式子成立，即

u（x，y，t）-u^h（x_i，y_j，t_k）≤Cτ²+‖u^（m+1）‖_∞el_x/2（m-2）^m-2+el_y/2（n-2）^n-2。（26）

證明：對方程時間方向采用Crank-Nicolson差分法離散得到的數值解為u^k+1=u（x，y，t_k+1），則

δu^（k+1）/2_t+u^（k+1）/2（u^（k+1）/2_x+u^（k+1）/2_y）-v（u^（k+1）/2_x，x+u^（k+1）/2_y，y）=R^k。（27）

式（27）中：δu^（k+1）/2_t=1/τ（u（x，y，t_k+1）-u（x，y，t_k）），R^k=δu^（k+1）/2_t-u^（k+1）/2_t是時間方向的截斷誤差。

由Taylor展開，可得

R^k≤Cτ²。（28）

空間方向利用重心Lagrange插值配點法離散，有/δu^h_t（x_i，y_j，t_k）+u（x_i，y_j，t_（k+1）/2）（u_x（x_i，y_j，t_（k+1）/2）+u_y（x_i，y_j，t_（k+1）/2））-/

v（u_x，x（x_i，y_j，t_（k+1）/2）+u_y，y（x_i，y_j，t_（k+1）/2））=R^k+ε^i，j。（29）

式（29）中：ε^i，j為空間截斷誤差。

由式（27），（29），可以得到

/δu^（k+1）/2_t-δu^h_t（x_i，y_j，t_k）+u^（k+1）/2（u^（k+1）/2_x+u^（k+1）/2_y）-/

u（x_i，y_j，t_（k+1）/2）（u_x（x_i，y_j，t_（k+1）/2）+u_y（x_i，y_j，t_（k+1）/2））-/

v（u^（k+1）/2_x，x+u^（k+1）/2_y，y）+v（u_x，x（x_i，y_j，t_（k+1）/2）+u_y，y（x_i，y_j，t_（k+1）/2））=-ε^i，j。（30）

類似定理2的推導，有

ε^i，j≤C‖u^（m+1）‖_∞el_x/2（m-2）^m-2+el_y/2（n-2）^n-2。（31）

結合式（28），（31），則定理得證。

記u^k（x，y）=u（x，y，t_k），k=0，1，…，K，采用Crank-Nicolson差分法對時間離散，對方程中的非線性項線性化處理，得到時間半離散格式，即

u^k+1-u^k/τ+u^k（u^k+1_x+u^k+1_y）+u^k+1（u^k_x+u^k_y）/2=v（u^k+1_x，x+u^k+1_y，y+u^k_x，x+u^k_y，y）/2。（32）

將微分矩陣代入時間半離散格式，有Burgers方程全離散格式，即

I+τ/2（diag（U^k）A^*+diag（A^*U^k））-vτ/2B^*U^k+1=（I+vτ/2B^*）U^k。（33）

3 數值算例

為驗證重心插值配點法求解Burgers方程的有效性和合理性，通過算例對重心插值配點法穩定性進行驗證，與Crank-Nicolson差分法對比說明重心插值配點法的高精度。記U_e，U^h分別為精確解及數值解，定義在無窮范數意義下的絕對誤差和相對誤差分別為

E_∞=U_e-U^h_∞， E_r=U_e-U^h_∞/U_e∞。

上式中：·_∞表示無窮范數。

算例1 考慮d=1時的一維Burgers方程（1），取T=1，a=0，b=1，初始條件為

當v=0.01，固定空間節點n=16時，改變時間步長，可得到重心插值配點法的時間收斂階，如表1所示。由表1可知：隨著K的增加，利用重心插值配點法求出的方程數值解與精確解的絕對誤差的數量級從10^-3減到 10^-5，因此，時間收斂階為二階精度。

當固定時間步長τ=0.001時，改變空間節點數量，可得到重心插值配點法與有限差分法的誤差，如表2所示。由表2可知： 重心插值配點法在空間選取10個節點時，最大絕對誤差即可達到10^-6量級，而有限差分法空間選取128個節點時，最大絕對誤差才能達到10^-6量級，這表明重心插值配點法在空間上用較少的點可達到更高的精度。

重心插值配點法數值解與精確解（τ=0.05，n=20），如圖1所示。由圖 1 可知：兩種解吻合。

絕對誤差（τ=0.05，n=20），如圖2所示。空間收斂階（τ=0.001，v=0.01），如圖3所示。由圖3可知：重心插值配點法具有指數收斂速度。不同參數下重心插值配點法的數值解，如圖4所示。

圖4（a）為最終時間T=1，2，4，8，16的數值解圖像，圖4（b）為不同的雷諾數Re=10，20，50，100的數值解圖像。由圖4可知：隨著T，Re的增大，數值解保持穩定，說明了數值法的穩定性。

重心插值配點法的數值解與精確解的等勢圖（v=0.1），如圖5所示。

算例2 考慮如下二維Burgers方程，即

u_t+u?u/?x+?u/?y=v?²u/?x²+?²u/?y²，（x，y）∈［0，1］×［0，1］，t∈［0，1］。

初值條件和邊值條件由精確解確定，方程的精確解為

重心插值配點法時間收斂階，如表3所示。

當v=0.1，空間節點數m=n，固定空間節點數n=16時，改變時間節點數，重心插值配點法與有限差分法的誤差（τ=0.001），如表4所示。由表4可知：當n取8時，重心插值配點法的誤差具有10^-4量級，而有限差分法的誤差只有10^-3量級；重心插值配點法n=12時誤差有10^-5量級，而有限差分法在n=64時誤差只有10^-5量級。

重心插值配點法的數值解與精確解（τ=0.001，n=32），如圖6所示。由圖6可知：數值解與精確解是相符的。

絕對誤差（τ=0.001，n=32），如圖7所示。由圖7可知：數值解與精確解是相符的。

空間收斂階（τ=0.001），如圖8所示。由圖8可知：重心插值配點法在求解二維問題時，空間具有指數收斂。

算例2結果表明：重心插值配點法在二維問題中，具有高精度性與有效性，且可通過更少的節點得到相對高的精度，與理論分析相符。

5 結束語

將Crank-Nicolson差分法與重心插值切比雪夫配點法相結合，提出了一種求解Burgers的有效數值法，對方程的非線性項進行線性化處理，避免迭代帶來大計算量，并給出了全離散法的相容性分析。通過數值算例驗證數值法的有效性及高精度，通過與Crank-Nicolson差分法進行比較，重心插值法可以用較少的節點得到較高的精度。

參考文獻：

［1］ BURGERS J M.A mathematical model illustrating the theory of turbulence［J］.Advances in Applied Mechanics，1948，1：171-199.DOI：10.1016/S0065-2156（08）70100-5.

［2］ 陳蓮.一維Burgers方程的幾種有限差分解法［D］.南充：西華師范大學，2020.

［3］ WANG Xuping，ZHANG Qifeng，SUN Zhizhong.The pointwise error estimates of two energy-preserving fourth-order compact schemes for viscous Burgers equation［J］.Advances in Com-putational Mathematics，2021，47（2）：23.DOI：10.1007/S10444-021-09848-9.

［4］ ZHANG Qifeng，QIN Yifan，SUN Zhizhong.Linearly compact scheme for 2D Sobolev equation with Burgers′ type nonlinearity［J］.Numerical Algorithms，2022，91（3）：1081-1114.DOI：10.1007/S11075-022-01293-Z.

［5］ XU Chao，PEU Lifang.Unconditional superconvergence analysis of two modified finite element fully discrete schemes for nonlinear Burgers′ equation［J］.Applied Numerical Mathematics，2023，185：1-17.DOI：10.1016/j.apnum.2022.11.008.

［6］ WANG Chuan，WANG Tianjun.A multi-domain Galerkin method with numerical integration for the Burgers equation［J］.International Journal of Computer Mathematics，2023，100（5）：927-947.DOI：10.1080/00207160.2023.2171265.

［7］ ZHAO Zhihui，LI Hong.Numerical study of two-dimensional Burgers equation by using a continuous Galerkin method［J］.Computers and Mathematics with Applications，2023，149（1）：38-48.DOI：10.1016/J.CAMWA.2023.08.030.

［8］ WANG Haifeng，XU Da，ZHOU Jun，et al.Weak Galerkin finite element method for a class of time fractional generalized Burgers equation［J］.Numerical Methods for Partial Differential Equations，2021，37（1）：732-749.DOI：10.1002/num.22549.

［9］ 于孟文，張新東.重心插值配點法求解Volterra積分方程［J］.新疆師范大學學報（自然科學版），2023，42（1）：75-80.DOI：10.14100/j.cnki.1008-9659.2023.01.010.

［10］ DENG Yangfang，WENG Zhifeng.Barycentric interpolation collocation method based on Crank-Nicolson scheme for the Allen-Cahn equation［J］.AIMS Mathematics，2021，6（4）：3857-3873.DOI：10.3934/MATH.2021229.

［11］ 黃蓉，鄧楊芳，翁智峰.SAV/重心插值配點法求解Allen-Cahn方程［J］.應用數學和力學，2023，44（5）：573-582.DOI：10.21656/1000-0887.430149.

［12］ 鄧楊芳，姚澤豐，汪精英，等.二維Allen-Cahn方程的有限差分法/配點法求解［J］.華僑大學學報（自然科學版），2020，41（5）：690-694.DOI：10.11830/ISSN.1000-5013.202001001.

［13］ 翁智峰，姚澤豐，賴淑琴.重心插值配點法求解Allen-Cahn方程［J］.華僑大學學報（自然科學版），2019，40（1）：133-140.DOI：10.11830/ISSN.1000-5013.201806043.

［14］ 賴舒琴，華之維，翁智峰.重心插值配點法求解Black-Scholes方程［J］.聊城大學學報（自然科學版），2020，33（5）：1-8.DOI：10.19728/j.issn1672-6634.2020.05.001.

［15］ 黃蓉，翁智峰.時間分數階Allen-Cahn方程的重心插值配點法［J］.華僑大學學報（自然科學版），2022，43（4）：553-560.DOI：10.11830/ISSN.1000-5013.202104060.

［16］ LI Jin，SU Xiaoning，ZHAO Kaiyan.Barycentric interpolation collocation algorithm to solve fractional differential equations［J］.Mathematics and Computers in Simulation，2023，205：340-367.DOI：10.1016/J.MATCOM.2022.10.005.

［17］ ?MER O.Application of a collocation method based on linear barycentric interpolation for solving 2D and 3D Klein-Gordon-Schr?dinger （KGS） equations numerically［J］.Enginee-Ring Computations，2021，38（5）：2394-2414.DOI：10.1108/EC-06-2020-0312.

［18］ YI Shichao，YAO Linquan.A steady barycentric Lagrange interpolation method for the 2D higher order time fractional teleglaph equation with nonlocal boundary condition with error analysis［J］.Numerical Methods for Partial Differential Equations，2019，35（5）：1694-1716.DOI：10.1002/num.22371.

（責任編輯：陳志賢 英文審校：黃心中）