摘要: 利用重心插值構(gòu)造求解Burgers方程的無網(wǎng)格數(shù)值方法。在時間方向用Crank-Nicolson差分法對方程進行離散,在空間方向用重心插值切比雪夫配點法逼近函數(shù)本身及其空間導(dǎo)數(shù)。 對全離散數(shù)值格式進行相容性分析。數(shù)值算例表明:與經(jīng)典的有限差分方法比較,重心插值配點法用較少的節(jié)點就能達到較高的精度。
關(guān)鍵詞: 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)中:Ω Rd(d=1,2)為有界凸區(qū)域,邊界分段光滑。
Burgers方程在流體力學(xué)、氣體動力學(xué)、交通流等眾多領(lǐng)域都有著重要的作用和地位。許多學(xué)者對該問題做了大量研究,并提出了數(shù)值解法。陳蓮[2]建立兩種求解Burgers方程的Crank-Nicolson差分法格式,并對其進行誤差分析。Wang等[3]建立一種求解粘性Burgers方程的能量守恒型四階隱式緊差分格式。Zhang等[4]利用線性化緊致差分格式求解具有Burgers型非線性項的二維Sobolev方程。Xu等[5]建立兩種求解非線性Burgers方程的修正的非協(xié)調(diào)有限元全離散格式,并對其進行超收斂分析。Wang等[6]建立一種求解Burgers方程的多區(qū)域Galerkin方法,并分析其收斂性及穩(wěn)定性。Zhao等[7]利用時空連續(xù)Galerkin方法對二維Burgers方程進行數(shù)值求解,并給出先驗誤差估計。Wang等[8]利用弱伽遼金有限元方法求解一類時間分數(shù)階的廣義Burgers方程。
以上求解Burgers方程的數(shù)值方法都是基于網(wǎng)格剖分建立離散格式,Chebyshev節(jié)點的重心Lagrange插值公式可以有效克服Runge現(xiàn)象。很多學(xué)者將其推廣到求解各類微分方程與積分方程,如Volterra積分方程[9]、Allen-Cahn方程[10-13]、Black-Scholes方程[14]、分數(shù)階微分方程[15-16]等。?mer[17]應(yīng)用重心插值配點法數(shù)值求解了二維和三維Klein-Gordon-Schr?dinger(KGS)方程?;诖耍疚膶?yīng)用重心插值的Burgers方程數(shù)值解法進行研究。
1 重心Lagrange插值
1.1 重心Lagrange插值函數(shù)
設(shè)插值節(jié)點為xj,對應(yīng)的實數(shù)為yj(j=0,1,…,n),使用多項式插值,則在次數(shù)不超過n的多項式空間中,存在唯一的重心Lagrange插值函數(shù)p(x),滿足p(xj)=yj,j=0,1,…,n。重心Lagrange插值函數(shù)為
p(x)=∑n/j=0wj/(x-xj)yj/∑n/i=0wi/(x-xi)=∑n/j=0Dj(x)yj。(2)
式(2)中:Dj(x)為重心Lagrange插值基函數(shù);wj為重心Lagrange插值函數(shù)的重心權(quán),wj為
wj=∏n/i=0,i≠j(xj-xi)-1。(3)
由于重心Lagrange插值函數(shù)的重心權(quán)只與插值節(jié)點分布有關(guān),在利用重心插值逼近未知函數(shù)時,采用第二類 Chebyshev節(jié)點,即
xj=cosj/nπ, j=0,1,…,n。(4)
可以得到重心Lagrange插值函數(shù)的重心權(quán)為
wj=(-1)δj,δj=1/2, j=0或n,1, 其他。(5)
對式(2)進行求導(dǎo),得到函數(shù)p(x)在插值節(jié)點xi處的導(dǎo)數(shù),ps(xi)為
ps(xi)=dsp(xi)/dxs=∑n/j=0D(s)j(xi)yj=∑n/j=0D(s)i,jyj。(6)
將其寫成矩陣形式即為
p(s)=D(s)p。(7)
式(7)中:D(s)為關(guān)于節(jié)點xi(i=0,1,…,n)的s階微分矩陣,D(s)=(D(s)i,j)n×n,s=1,2。
通過對重心Lagrange插值基函數(shù)求導(dǎo),分別得到一、二階微分矩陣的元素[11]為
D(1)i,j=ωj/ωi/xi-xj, 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/xi-xj), i≠j,D(2)i,i=-∑n/j=0,j≠iD(2)i,j。(8)
1.2 逼近性質(zhì)
設(shè)p(x,y)是u(x,y)的重心Lagrange插值函數(shù),滿足p(xm,yn)=u(xm,yn),則定義誤差函數(shù)為
e(x,y)=u(x,y)-p(x,y)。(9)
引理1[18] 如果u(x,y)∈Cm+1([a,b]×[c,d]),其中,m=max{m,n},則
|e(x,y)|≤u(m+1)∞(C1elx/2mm+C2ely/2nn)。(10)
式(10)中:e為自然對數(shù);lx代表區(qū)間[a,b]的一半;ly代表區(qū)間[c,d]的一半。
同理可以得到
ex(x,y)≤C*1‖?(m+1)xu‖∞elx/2(m-1)m-1+C2‖?(n+1)yu‖∞ely/2nn,ex,x(x,y)≤C**1‖?(m+1)xu‖∞elx/2(m-2)m-2+C2‖?(n+1)yu‖∞ely/2nn,ey(x,y)≤C1‖?(m+1)xu‖∞elx/2mm+C*2‖?(n+1)yu‖∞ely/2(n-1)n-1,ey,y(x,y)≤C1‖?(m+1)xu‖∞elx/2mm+C**2‖?(n+1)yu‖∞ely/2(n-2)n-2。(11)
式(11)中:C1,C*1,C**1,C2,C*2,C**2為常數(shù)。
2 二維Burgers方程的離散格式
考慮二維Burgers方程,對時間區(qū)域進行K等分,記步長為τ=T/K,T為時間計算的終點??臻g區(qū)域采用第二類Chebyshev節(jié)點分別離散為m+1個計算節(jié)點a=x0lt;x1lt;…lt;xm=b和n+1個計算節(jié)點 c=y0lt;y1lt;…lt;yn=d,得到計算節(jié)點(xi,yj,tk),0≤i≤m,0≤j≤n,0≤k≤K。
2.1 半離散格式及相容性分析
記函數(shù)u(x,y,t)在節(jié)點x0,x1,…,xm上的值為u(xi,y,t)=ui(y,t),則函數(shù)u(x,y,t)在節(jié)點x0,x1,…,xm上的重心插值函數(shù)為
u(x,y,t)=∑m/s=0Ls(x)us(y,t)。(12)
式(12)中:Ls(x)表示x方向的重心插值基函數(shù)。
將式(12)代入方程(1)中,并令方程在點x0,x1,…,xm上精確成立,得到常微分方程組,即
∑m/s=0Lsxi?us(y,t)/?t+/∑m/s=0Lsxius(y,t)∑m/s=0L′sxius(y,t)+∑m/s=0Lsxi?2us(y,t)/?y2=/
v∑m/s=0L″sxius(y,t)+∑m/s=0Lsxi?2us(y,t)/?y2。(13)
式(13)中:L′s(xi)=C(1)s,i為關(guān)于節(jié)點x0,x1,…,xm的一階微分矩陣C(1)中的元素,類似地,L″s(xi)=C(2)s,i為二階微分矩陣C(2)中的元素。
記ui,j(t)=ui(yj,t),函數(shù)ui(y,t)在節(jié)點y0,y1,…,yn的重心插值函數(shù)為
ui(y,t)=∑n/r=0Tr(y)ui,r(t)。(14)
式(14)中:Tr(y)表示y方向的重心插值基函數(shù)。
將式(14)代入式(13),微分矩陣形式為
dU(t)/dt+diag(U(t))(C(1)In+1+Im+1D(1))U(t)=v(C(2)In+1+Im+1D(2))U(t)。(15)
式(15)中:
U(t)=[U0(t)…Um(t)]T;
Ui(t)=[ui,0(t)…ui,n(t)]T(i=0,…,m);
D(k)為關(guān)于節(jié)點y0,y1,…,yn的k階微分矩陣。
為了方便進行相容性分析,算子定義為
Γu(x,y,t):=ut+u(ux+uy)-v(ux,x+uy,y)。(16)
令u(xm,yn,t)為u(x,y,t)的數(shù)值解,則可以得到
Γu(xm,yn,t)=0。(17)
定理1 若u(x,y,t)∈C(m+1)([a,b]×[c,d])×C(0)(0,T],其中,m=max{m,n},假設(shè)u(xm,yn,t)有界,則誤差估計為
u(x,y,t)-u(xm,yn,t)≤C**1?(m+1)xu∞elx/2(m-2)m-2+/
C**2‖?(n+1)yu‖∞ely/2(n-2)n-2。(18)
證明:Γu(x,y,t)-Γu(xm,yn,t)=ut(x,y,t)-ut(xm,yn,t)+/
u(x,y,t)ux(x,y,t)-u(xm,yn,t)ux(xm,yn,t)+/
u(x,y,t)uy(x,y,t)-u(xm,yn,t)uy(xm,yn,t)+/
v(uxx(xm,yn,t)-(x,y,t))+v(uyy(xm,yn,t)-uyy(x,y,t))=/R1+R2+R3+R4+R5。(19)
式(19)中:
R1=ut(x,y,t)-ut(xm,yn,t)=ut(x,y,t)-ut(xm,y,t)+
ut(xm,y,t)-ut(xm,yn,t)=et(xm,y,t)+et(xm,yn,t)。(20)
由引理1,有R1=/et(xm,y,t)+et(xm,yn,t)≤/C1‖?(m+1)xu‖∞elx/2mm+C2‖?(n+1)yu‖∞ely/2nn。(21)
R2的估計為
R2=/u(x,y,t)ux(x,y,t)-u(xm,yn,t)ux(xm,yn,t)≤/ux(x,y,t)u(x,y,t)-u(xm,yn,t)+/
u(xm,yn,t)ux(x,y,t)-ux(xm,yn,t)≤/
C*1‖?(m+1)xu‖∞elx/2(m-1)m-1+C2?(n+1)yu∞ely/2nn。(22) 同理有
R3=u(x,y,t)uy(x,y,t)-u(xm,yn,t)uy(xm,yn,t)≤
C1‖?(m+1)xu‖∞elx/2mm+C*2‖?(n+1)yu‖∞ely/2(n-1)n-1。 (23)
類似地,對于R4,R5的估計分別為
R4=vex,x(xm,y,t)+vex,x(xm,yn,t)C≤
C**1‖?(m+1)xu‖∞(elx/2(m-2))m-2+C2‖?(n+1)yu‖∞-ely/2nn。(24)
R5=vey,y(xm,y,t)+vey,y(xm,yn,t)≤
C1‖?(m+1)xu‖∞elx/2mm+C**2‖?(n+1)yu‖∞ely/2(n-2)n-2。 (25)
將式(21)~(25)代入式(19),則定理得證。
2.2 全離散格式及相容性分析
定理2 設(shè)u(x,y,t)∈C(m+1)([a,b]×[c,d])×C(2)(0,T],假設(shè)u(xm,yn,t)有界,使用重心插值配點法求解的數(shù)值解為uh(xi,yj,tk),則如下式子成立,即
u(x,y,t)-uh(xi,yj,tk)≤Cτ2+‖u(m+1)‖∞elx/2(m-2)m-2+ely/2(n-2)n-2。(26)
證明:對方程時間方向采用Crank-Nicolson差分法離散得到的數(shù)值解為uk+1=u(x,y,tk+1),則
δu(k+1)/2t+u(k+1)/2(u(k+1)/2x+u(k+1)/2y)-v(u(k+1)/2x,x+u(k+1)/2y,y)=Rk。(27)
式(27)中:δu(k+1)/2t=1/τ(u(x,y,tk+1)-u(x,y,tk)),Rk=δu(k+1)/2t-u(k+1)/2t是時間方向的截斷誤差。
由Taylor展開,可得
Rk≤Cτ2。(28)
空間方向利用重心Lagrange插值配點法離散,有/δuht(xi,yj,tk)+u(xi,yj,t(k+1)/2)(ux(xi,yj,t(k+1)/2)+uy(xi,yj,t(k+1)/2))-/
v(ux,x(xi,yj,t(k+1)/2)+uy,y(xi,yj,t(k+1)/2))=Rk+εi,j。(29)
式(29)中:εi,j為空間截斷誤差。
由式(27),(29),可以得到
/δu(k+1)/2t-δuht(xi,yj,tk)+u(k+1)/2(u(k+1)/2x+u(k+1)/2y)-/
u(xi,yj,t(k+1)/2)(ux(xi,yj,t(k+1)/2)+uy(xi,yj,t(k+1)/2))-/
v(u(k+1)/2x,x+u(k+1)/2y,y)+v(ux,x(xi,yj,t(k+1)/2)+uy,y(xi,yj,t(k+1)/2))=-εi,j。(30)
類似定理2的推導(dǎo),有
εi,j≤C‖u(m+1)‖∞elx/2(m-2)m-2+ely/2(n-2)n-2。(31)
結(jié)合式(28),(31),則定理得證。
記uk(x,y)=u(x,y,tk),k=0,1,…,K,采用Crank-Nicolson差分法對時間離散,對方程中的非線性項線性化處理,得到時間半離散格式,即
uk+1-uk/τ+uk(uk+1x+uk+1y)+uk+1(ukx+uky)/2=v(uk+1x,x+uk+1y,y+ukx,x+uky,y)/2。(32)
將微分矩陣代入時間半離散格式,有Burgers方程全離散格式,即
I+τ/2(diag(Uk)A*+diag(A*Uk))-vτ/2B*Uk+1=(I+vτ/2B*)Uk。(33)
3 數(shù)值算例
為驗證重心插值配點法求解Burgers方程的有效性和合理性,通過算例對重心插值配點法穩(wěn)定性進行驗證,與Crank-Nicolson差分法對比說明重心插值配點法的高精度。記Ue,Uh分別為精確解及數(shù)值解,定義在無窮范數(shù)意義下的絕對誤差和相對誤差分別為
E∞=Ue-Uh∞, Er=Ue-Uh∞/Ue∞。
上式中:·∞表示無窮范數(shù)。
算例1 考慮d=1時的一維Burgers方程(1),取T=1,a=0,b=1,初始條件為
當(dāng)v=0.01,固定空間節(jié)點n=16時,改變時間步長,可得到重心插值配點法的時間收斂階,如表1所示。由表1可知:隨著K的增加,利用重心插值配點法求出的方程數(shù)值解與精確解的絕對誤差的數(shù)量級從10-3減到 10-5,因此,時間收斂階為二階精度。
當(dāng)固定時間步長τ=0.001時,改變空間節(jié)點數(shù)量,可得到重心插值配點法與有限差分法的誤差,如表2所示。由表2可知: 重心插值配點法在空間選取10個節(jié)點時,最大絕對誤差即可達到10-6量級,而有限差分法空間選取128個節(jié)點時,最大絕對誤差才能達到10-6量級,這表明重心插值配點法在空間上用較少的點可達到更高的精度。
重心插值配點法數(shù)值解與精確解(τ=0.05,n=20),如圖1所示。由圖 1 可知:兩種解吻合。
絕對誤差(τ=0.05,n=20),如圖2所示??臻g收斂階(τ=0.001,v=0.01),如圖3所示。由圖3可知:重心插值配點法具有指數(shù)收斂速度。不同參數(shù)下重心插值配點法的數(shù)值解,如圖4所示。
圖4(a)為最終時間T=1,2,4,8,16的數(shù)值解圖像,圖4(b)為不同的雷諾數(shù)Re=10,20,50,100的數(shù)值解圖像。由圖4可知:隨著T,Re的增大,數(shù)值解保持穩(wěn)定,說明了數(shù)值法的穩(wěn)定性。
重心插值配點法的數(shù)值解與精確解的等勢圖(v=0.1),如圖5所示。
算例2 考慮如下二維Burgers方程,即
ut+u?u/?x+?u/?y=v?2u/?x2+?2u/?y2,(x,y)∈[0,1]×[0,1],t∈[0,1]。
初值條件和邊值條件由精確解確定,方程的精確解為
重心插值配點法時間收斂階,如表3所示。
當(dāng)v=0.1,空間節(jié)點數(shù)m=n,固定空間節(jié)點數(shù)n=16時,改變時間節(jié)點數(shù),重心插值配點法與有限差分法的誤差(τ=0.001),如表4所示。由表4可知:當(dāng)n取8時,重心插值配點法的誤差具有10-4量級,而有限差分法的誤差只有10-3量級;重心插值配點法n=12時誤差有10-5量級,而有限差分法在n=64時誤差只有10-5量級。
重心插值配點法的數(shù)值解與精確解(τ=0.001,n=32),如圖6所示。由圖6可知:數(shù)值解與精確解是相符的。
絕對誤差(τ=0.001,n=32),如圖7所示。由圖7可知:數(shù)值解與精確解是相符的。
空間收斂階(τ=0.001),如圖8所示。由圖8可知:重心插值配點法在求解二維問題時,空間具有指數(shù)收斂。
算例2結(jié)果表明:重心插值配點法在二維問題中,具有高精度性與有效性,且可通過更少的節(jié)點得到相對高的精度,與理論分析相符。
5 結(jié)束語
將Crank-Nicolson差分法與重心插值切比雪夫配點法相結(jié)合,提出了一種求解Burgers的有效數(shù)值法,對方程的非線性項進行線性化處理,避免迭代帶來大計算量,并給出了全離散法的相容性分析。通過數(shù)值算例驗證數(shù)值法的有效性及高精度,通過與Crank-Nicolson差分法進行比較,重心插值法可以用較少的節(jié)點得到較高的精度。
參考文獻:
[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].南充:西華師范大學(xué),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].新疆師范大學(xué)學(xué)報(自然科學(xué)版),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].應(yīng)用數(shù)學(xué)和力學(xué),2023,44(5):573-582.DOI:10.21656/1000-0887.430149.
[12] 鄧楊芳,姚澤豐,汪精英,等.二維Allen-Cahn方程的有限差分法/配點法求解[J].華僑大學(xué)學(xué)報(自然科學(xué)版),2020,41(5):690-694.DOI:10.11830/ISSN.1000-5013.202001001.
[13] 翁智峰,姚澤豐,賴淑琴.重心插值配點法求解Allen-Cahn方程[J].華僑大學(xué)學(xué)報(自然科學(xué)版),2019,40(1):133-140.DOI:10.11830/ISSN.1000-5013.201806043.
[14] 賴舒琴,華之維,翁智峰.重心插值配點法求解Black-Scholes方程[J].聊城大學(xué)學(xué)報(自然科學(xué)版),2020,33(5):1-8.DOI:10.19728/j.issn1672-6634.2020.05.001.
[15] 黃蓉,翁智峰.時間分數(shù)階Allen-Cahn方程的重心插值配點法[J].華僑大學(xué)學(xué)報(自然科學(xué)版),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.
(責(zé)任編輯:陳志賢 英文審校:黃心中)