• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    DFR法和LDG法解線性三階KdV方程的等價(jià)性

    2024-02-13 00:00:00畢卉李曉彤
    關(guān)鍵詞:三階等價(jià)插值

    摘 要:研究了直接通量重構(gòu)方法(Direct Flux Reconstruction,簡稱DFR)和局部間斷Galerkin方法(Local discontinuous Galerkin,簡稱LDG)求解線性三階KdV方程的等價(jià)性問題。首先分別利用DFR法和LDG法對線性三階KdV方程進(jìn)行空間離散并給出兩種數(shù)值算法的空間離散格式,其次用兩種方法證明了DFR法和LDG法求解線性三階KdV方程的等價(jià)性。第一種方法借助高斯求積的性質(zhì):M點(diǎn)高斯求積具有2M-1階代數(shù)精度;第二種方法利用洛巴托多項(xiàng)式的特殊性質(zhì):M+1次洛巴托多項(xiàng)式導(dǎo)數(shù)的零點(diǎn)是M個高斯點(diǎn)。最后以M=1為例,給出了兩種數(shù)值算法求解線性三階KdV方程的離散常微分方程組,驗(yàn)證了兩種方法的等價(jià)性。

    關(guān)鍵詞:直接通量重構(gòu)法;局部間斷Galerkin法;線性三階KdV方程;高斯求積;洛巴托多項(xiàng)式

    DOI:10.15938/j.jhust.2024.05.016

    中圖分類號: O241.3

    文獻(xiàn)標(biāo)志碼: A

    文章編號: 1007-2683(2024)05-0142-07

    Equivalence Between DFR Method and LDG Method for Solving Linear Third-Order KdV Equation

    BI Hui, LI Xiaotong

    (School of Sciences, Harbin University of Science and Technology, Harbin 150080, China)

    Abstract:The equivalence between direct flux reconstruction method and local discontinuous Galerkin method in solving linear third-order KdV equation is studied. Firstly, the linear third-order KdV equation is spatially dispersed by DFR method and LDG method respectively, and the spatial discretization schemes of two numerical algorithms are given. Secondly, the equivalence of DFR and LDG methods is proved by two methods. The first proof relies on the property of Gauss quadrature: the M-point Gauss quadrature has 2M-1 order algebraic accuracy; the second proof takes advantage of the special property of the Lobatto polynomial: the zeros of the derivative of Lobatto polynomial of degree M+1 are M Gauss points. At last, taking the value of M to be 1 as an example, it is shown that the two numerical algorithms are equivalent to the ordinary differential equations used for programming when solving the linear third-order KdV equation.

    Keywords:direct flux reconstruction method; local discontinuous Galerkin method; linear third-order KdV equation; Gauss quadrature; Lobatto polynomial

    0 引 言

    1973年,Reed和Hill[1在研究中子運(yùn)輸問題時提出了一種將投影作為近似解的數(shù)值解法——間斷Galerkin(discontinuous galerkin,簡稱DG)法。該法具有高階精度和易于處理復(fù)雜幾何形狀或復(fù)雜邊界問題的優(yōu)點(diǎn),但是對于含有高階偏導(dǎo)數(shù)的偏微分方程,直接應(yīng)用DG法很難得到相容的數(shù)值格式。1998年,在Bassi和Rebay[2處理可壓縮的Navier-Stokes方程的啟發(fā)下,Cockburn和Shu[3在求解對流擴(kuò)散方程時提出了局部間斷Galerkin(local discontinuous galerkin,簡稱LDG)法。該方法的主要思想是通過引入輔助變量把原有的高階偏微分方程轉(zhuǎn)化為等價(jià)的一階偏微分方程組,再對所得到的一階偏微分方程分別使用DG法,作為DG法的推廣,可以靈活的處理含有更高階偏導(dǎo)數(shù)的偏微分方程,例如:含有三階偏導(dǎo)數(shù)的KdV方程4、含有四階偏導(dǎo)數(shù)的雙調(diào)和方程5、非線性波動方程6-7、含有對流項(xiàng)和擴(kuò)散項(xiàng)的四階線性偏微分方程8等。更多關(guān)于DG法和LDG法的研究成果可以查閱文[9-15]。

    2007年,Huynh[16在求解守恒律方程時提出了一種將插值和投影結(jié)合起來的高階精度的數(shù)值解法——通量重構(gòu)(flux reconstruction,簡稱FR)法。該方法的靈活性體現(xiàn)在先用插值法近似構(gòu)造通量函數(shù),再選擇一個校正函數(shù)將分段的不連續(xù)通量重建為全局連續(xù)通量,但是使用校正函數(shù)時需要許多不同的計(jì)算步驟,使得該方法的復(fù)雜性和計(jì)算費(fèi)用大大增加。2016年,Romero、Asthana和Jameson[17提出了一種通量重構(gòu)法的簡化方法——直接通量重構(gòu)(direct flux reconstruction,簡稱DFR)法,同時證明了DFR法和FR法求解雙曲守恒律方程的等價(jià)性。該方法是一種僅使用插值法的數(shù)值解法,更多關(guān)于FR法和DFR法的研究成果可以查閱文[18-23]。2020年,Huynh使用更簡潔的方法證明了DFR法、FR法和DG法求解一維非線性雙曲守恒律方程的等價(jià)性24,求解拋物方程、線性對流擴(kuò)散方程25的等價(jià)性也有了嚴(yán)格的數(shù)學(xué)證明。但是,求解含有更高階偏導(dǎo)數(shù)的偏微分方程是否具有等價(jià)性仍需進(jìn)一步求證。

    Korteweg-de Vries(簡稱 KdV)方程最初是用于研究一種單向運(yùn)動淺水波的偏微分方程,隨后又被用來描述超流費(fèi)米氣體、離子聲波等不同類型的物理現(xiàn)象。該方程在流體力學(xué)、水文學(xué)、物理學(xué)等領(lǐng)域中有著廣泛應(yīng)用,因此,對于KdV方程數(shù)值方法的研究一直是偏微分方程數(shù)值解法的重要課題之一。本文研究了DFR法和LDG法求解含有三階偏導(dǎo)數(shù)的線性KdV方程的等價(jià)性問題,并列舉了M=1時兩種方法對方程進(jìn)行空間離散后得到的常微分方程組,更直觀的展示兩種數(shù)值方法的等價(jià)。

    本文的結(jié)構(gòu)如下:第一節(jié)預(yù)備知識介紹了本文證明等價(jià)性時所用到的符號、線性變換、拉格朗日插值基函數(shù);第二節(jié)首先給出線性三階KdV方程的DFR法格式和LDG法格式,其次用兩種方法證明DFR法和LDG法求解線性三階KdV方程的等價(jià)性,最后以M=1為例,給出兩種數(shù)值算法用于編程的常微分方程組;最后一節(jié)給出結(jié)論和工作展望。

    1 預(yù)備知識

    1.1 網(wǎng)格剖分和有限元空間

    考慮計(jì)算區(qū)域Ω=[0,2π],并對其進(jìn)行如下的網(wǎng)格剖分:

    0=x1/2lt;x3/2lt;…lt;xN+1/2=2π

    當(dāng)j∈{1,2,…,N}時,將剖分出來的第j個區(qū)間記為Ej=[xj-1/2,xj+1/2],區(qū)間中點(diǎn)和區(qū)間長度分別記為xj=(xj-1/2+xj+1/2)/2,hj=(xj+1/2-xj-1/2)。數(shù)值解空間定義為:

    Vh={v∈L2(Ω):v|Ej∈PM-1j(Ej),j=1,2,…,N}

    其中PM-1j(Ej)為定義在區(qū)間Ej上次數(shù)不超過M-1次的多項(xiàng)式所構(gòu)成的空間,L2(Ω)為定義在Ω上平方Lebesgue可積函數(shù)所構(gòu)成的空間。對于任意的v∈Vh,定義函數(shù)v(x)在點(diǎn)xj+1/2處的左、右極限分別為v-j+1/2,v+j+1/2,函數(shù)v(x)在Ej的邊界點(diǎn)上的跳躍和平均值為[v]=v+j+1/2-v-j+1/2和{v}=(v+j+1/2+v-j+1/2)/2。

    1.2 線性變換

    設(shè)I=[-1,1],將I上次數(shù)不超過M-1次的多項(xiàng)式所構(gòu)成的空間記作PM-1(I)。對于I上任意的點(diǎn)ξ可以通過線性變換:

    x=ξhj/2+xj(1)

    得到Ej上與ξ對應(yīng)的點(diǎn)x;對于I上任意的可微函數(shù)rI(ξ),可以通過線性變換:

    rj(x)=rI(ξ(x))=rI(ξ)(2)

    得到Ej上與rI(ξ)對應(yīng)的函數(shù)rj(ξ);rj(ξ)的導(dǎo)數(shù)可以通過線性變換:

    drjdx=2hjdrIdξ(3)

    得到。

    1.3 兩組拉格朗日插值基函數(shù)

    設(shè)v(ξ,t)是[-1,1]×[0,T)上的一個二元函數(shù),規(guī)定:ξ0=-1,ξM+1=1,I上M次勒讓德多項(xiàng)式的零點(diǎn)為ξ1,ξ2,…,ξM。若將ξ1,ξ2,…,ξM作為I上M點(diǎn)高斯求積節(jié)點(diǎn),那么其也被稱為高斯點(diǎn)。若固定v的第一個變量,此時得到的v是一個一元函數(shù),如果用v0表示v(ξ0,t),那么同樣的方法可以得到v1,v2,…,vM+1

    當(dāng)n∈{1,2,…,M}時,分別以ξ1,ξ2,…,ξM與ξ0,ξ1,…,ξM+1為插值節(jié)點(diǎn)構(gòu)造兩組拉格朗日插值基函數(shù):

    n=∏Mm=1,m≠n(ξ-ξm)(ξn-ξm)(4)

    ^n=∏M+1m=0,m≠n(ξ-ξm)(ξn-ξm)(5)

    此時v在I上的M點(diǎn)與M+2點(diǎn)拉格朗日插值近似分別為

    vM=∑Mn=1vnn(6)

    vM+2=∑M+1n=0vn^n(7)

    當(dāng)j∈{1,2,….N}時,設(shè)uj是Ej×[0,T)上的一個二元函數(shù),通過式(1)~(3)可以得到Ej上M次勒讓德多項(xiàng)式的零點(diǎn)為xj,1,xj,2,…,xj,M,并規(guī)定:x0=xj-1/2,xM+1=xj+1/2。

    當(dāng)n∈{1,2,…,M}時,通過式(1)~(5)可得到分別以xj,1,xj,2,…,xj,M與xj,0,xj,2,…,xj,M+1為插值節(jié)點(diǎn)的兩組拉格朗日插值基函數(shù):

    j,n=∏Mm=1,m≠n(x-xj,m)(xj,n-xj,m)(8)

    ^j,n=∏M+1m=0,m≠n(x-xj,m)(xj,n-xj,m)(9)

    此時uj在Ej上的M點(diǎn)與M+2點(diǎn)拉格朗日插值近似分別為

    uj,M=∑Mn=1uj,nj,n(10)

    uj,M+2=∑M+1n=0uj,n^j,n(11)

    其中uj,n=uj(xj,n,t)。

    2 DFR法和LDG法解線性三階KdV方程的等價(jià)性

    2.1 線性三階KdV方程的DFR法格式

    考慮如下格式的線性三階KdV方程:

    ut+ux+uxxx=0(x,t)∈[0,2π]×[0,T)

    u(x,0)=u0(x)x∈[0,2π](12)

    假設(shè)初值u0(x)充分光滑且方程的精確解u(x,t)滿足周期性邊界條件。

    DFR法對方程式(12)進(jìn)行空間離散的具體步驟如下:

    第一步,按照本文1.1節(jié)的要求去處理方程式(12),按照1.3節(jié)的要求構(gòu)造兩組拉格朗日插值基函數(shù)。在每個Ej上用一個次數(shù)不超過M-1次的多項(xiàng)式uj去表達(dá)DFR法的數(shù)值解,由式(10)得:

    uj=∑Mn=1uj,nj,n(13)

    用一個次數(shù)不超過M+1次的多項(xiàng)式Fj近似表達(dá)-(uj+ujxx)。由式(11)得:

    Fj=∑M+1n=0Fj,n^j,n(14)

    Fj的具體構(gòu)造如下:

    Fj,0=F(xj,0,t)=-(j+j)|xj,0

    Fj,n=F(xj,n,t)=-(uj,n+qj,n) n∈{1,2,…,M}

    Fj,M+1=F(xj,M+1,t)=-(j+j)|xj,M+1

    其中qj是定義在Ej上的用來近似表達(dá)pjx的不超過M-1次的多項(xiàng)式,pj是定義在Ej上的用來近似表達(dá)ujx的不超過M-1次的多項(xiàng)式,pj與qj的具體構(gòu)造如下:

    pj=∑Mn=1pj(xj,n,t)j,n=∑Mn=1pj,nj,n

    qj=∑Mn=1qj(xj,n,t)j,n=∑Mn=1qj,nj,n

    pj,n=ujx(xj,n,t)+2hjwn[(j-u-j)-j,n|xj+1/2-(j-u+j)+j,n|xj-1/2

    qj,n=pjx(xj,n,t)+2hjwn[(j-p-j)-j,n|xj+1/2-(j-p+j)+j,n|xj-1/2

    其中j,j,j為數(shù)值通量,選取方法如下:

    j=u-j,j=p+j,j=q+j(15)

    第二步,令uj和Fj在Ej上的xj,1,xj,2,…,xj,M點(diǎn)處滿足方程(12)式:

    (uj)t(xj,n,t)=(Fj)x(xj,n,t) n∈{1,2,…,M}(16)

    根據(jù)拉格朗日插值基函數(shù)的性質(zhì),j,n僅在點(diǎn)xj,n處取值為1,其余點(diǎn)處取值為0可將(16)式化簡為

    (uj,n)t=(Fj)x(xj,n,t) n∈{1,2,…,M}(17)

    式(17)即為DFR法對方程式(12)的空間離散格式,它是一個由M個常微分方程所組成的常微分方程組,利用Runge-Kutta法求解式(17)可以得到Ej上的DFR法數(shù)值解uj。將u1,u2,…,uN連起來就可以得到[0,2π]×[0,T)上的DFR法數(shù)值解uh。因?yàn)閡h是在每個Ej上分別求得的,所以uh僅能保證在每個Ej上連續(xù),不能保證在整個定義域上連續(xù)。

    2.2 線性三階KdV方程的LDG法格式

    LDG法對方程(12)進(jìn)行空間離散的具體步驟如下:

    第一步,引入輔助變量p,q將含有三階偏導(dǎo)數(shù)的方程(12)等價(jià)的改寫為僅含有一階偏導(dǎo)數(shù)的偏微分方程組:

    ut=-(ux+qx)

    q=px

    p=ux(18)

    第二步,按照本文1.1節(jié)的要求處理式(18),在每個Ej上,先用一個次數(shù)不超過M-1次的多項(xiàng)式uj表示數(shù)值解,其次用一個次數(shù)不超過M-1次的多項(xiàng)式pj表示ujx,最后用一個次數(shù)不超過M-1次的多項(xiàng)式qj表示pjx。要求uj滿足:

    Ejujtj,ndx-∫Ej(uj+qj)j,nxdx=

    -(j+j)-j,n|xj+1/2+(j+j)+j,n|xj-1/2

    Ejqjj,ndx+∫Ejpjj,nxdx=

    j-j,n|xj+1/2-j+j,n|xj-1/2 n∈{1,2,…,M}

    Ejpjj,ndx+∫Ejujj,nxdx=

    j-j,n|xj+1/2-j+j,n|xj-1/2(19)

    其中j,j,j為數(shù)值通量,選取方法如下:

    j=u-j,j=p+j,j=q+j(20)

    本文選用以M點(diǎn)高斯點(diǎn)作為插值節(jié)點(diǎn)的拉格朗日插值基函數(shù)作為投影空間基函數(shù),這種方法稱為節(jié)點(diǎn)局部間斷有限元(nodalLDG法)法。

    式(19)即為nodalLDG法對方程(12)的空間離散格式,它是一個由M個常微分方程所組成的常微分方程組,利用Runge-Kutta法求解式(19)可以得到Ej上的LDG法數(shù)值解uj,被寫作式(13)的形式。將u1,u2,…uN連起來,就可以得到[0,2π]×[0,T)上的LDG法數(shù)值解uh。與DFR法數(shù)值解一樣,uh僅能保證在每個Ej上連續(xù),不能保證在整個定義域上連續(xù)。

    2.3 DFR法和LDG法求解線性三階KdV方程的等價(jià)性證明1

    定理1 DFR法和LDG法求解線性三階KdV方程(12)所得到的數(shù)值解等價(jià)。

    證明:由本文2.1節(jié)、2.2節(jié)對DFR法和LDG法解線性三階KdV方程的格式介紹可以看出要證明兩種數(shù)值解法等價(jià)只需證明:當(dāng)j∈{1,2,…,N}時,在每個Ej上,對于n∈{1,2,…,M},式(17)和式(19)等價(jià)即可。

    對于n∈{1,2,…,M},先對ujt和Fjx分別乘以j,n再在Ej上積分得到:

    Ejujtj,ndx=∫EjFjxj,ndx(21)

    下面我們分別證明式(21)和式(17)、式(19)都是等價(jià)的。對式(21)左右兩側(cè)應(yīng)用M點(diǎn)高斯求積,由高斯求積的性質(zhì),M點(diǎn)高斯求積具有2M-1階代數(shù)精度可以得到:

    wnhj2(uj,n)t=wnhj2(Fj)x(xj,n,t)(22)

    其中:wn為Ej上的M點(diǎn)高斯求積系數(shù)。對式(17)左右兩側(cè)同時乘以wnhj/2也可以得到式(22),即證得式(21)和式(17)等價(jià)。

    下證式(21)和式(19)等價(jià)。首先需要證明Fj與-(uj+ujxx)在Ej上的xj,1,xj,2,…,xj,M點(diǎn)處取值相等。對于n∈{1,2,…,M},對式(19)第三式和式(19)第二式的左側(cè)第二項(xiàng)使用分部積分可以得到:

    Ejpjj,ndx-∫Ejujxj,ndx=

    (j-u-j)-j,n|xj+1/2-(j-u+j)+j,n|xj-1/2(23)

    Ejqjj,ndx-∫Ejpjxj,ndx=

    (j-p-j)-j,n|xj+1/2-(j-p+j)+j,n|xj-1/2(24)

    對式(23)左側(cè)應(yīng)用M點(diǎn)高斯求積可得:

    wnhj2pj(xj,n)=∫Ejpjj,ndx=∫Ejujxj,ndx+

    (j-u-j)-j,n|xj+1/2-(j-u+j)+j,n|xj-1/2=

    wnhj2ujx(xj,n)+(j-u-j)-j,n|xj+1/2

    (j-u+j)+j,n|xj-1/2=wnhj2pj,n

    對式(24)左側(cè)應(yīng)用M點(diǎn)高斯求積可得:

    wnhj2qj(xj,n)=∫Ejqjj,ndx=∫Ejpjxj,ndx+

    (j-p-j)-j,n|xj+1/2-(j-p+j)+j,n|xj-1/2=

    wnhj2pjx(xj,n)+(j-p-j)-j,n|xj+1/2

    (j-p+j)+j,n|xj-1/2=wihj2qj,n

    于是有:

    Fj,n=Fj(xj,n,t)=-(uj,n+qj,n)=-

    (uj,n+pj,nx)=-(uj,n+pj,nx(xj,n))=-

    (uj,n+qj,n(xj,n))

    以上過程就證明了Fj與-(uj+ujxx)在Ej上的xj,1,xj,2,…,xj,M點(diǎn)處取值相等,上述分析也能夠說明DFR法與LDG法對qj與pj的定義是等價(jià)的。然后對式(21)右側(cè)使用分部積分可以得到:

    Ejujtj,ndx+∫EjFjj,nxdx=

    Fj,M+1j,n|xj+1/2-Fj,0+j,n|xj-1/2(25)

    對式(19)第一式和式(25)應(yīng)用M點(diǎn)高斯求積,由高斯求積的性質(zhì)和Fj的定義可知式(19)第一式和式(25)相等。綜上可知式(21)和式(19)第一式等價(jià)。

    以上證明了式(21)和式(17)、式(19)都是等價(jià)的,故式(17)和式(19)等價(jià)。證畢。

    這樣就完成了DFR法和LDG法求解線性三階KdV方程的第一種等價(jià)證明,證明的關(guān)鍵在于M點(diǎn)高斯求積具有2M-1階代數(shù)精度。接下來,通過先構(gòu)造多項(xiàng)式f^j使其滿足:在PM-2j上的投影與-(uj+qj)在PM-2j上的投影相同,在Ej邊界處的通量與-(uj+qj)在Ej邊界處的通量取值相同,再利用洛巴托多項(xiàng)式的特殊性質(zhì),證明FR方法分別與DFR法、LDG法等價(jià),給出第二種等價(jià)證明。

    2.4 DFR法和LDG法求解線性三階KdV方程的等價(jià)性證明2

    證明:定理1的第二種證明主要分兩步:

    第一步,對于每個Ej構(gòu)造一個多項(xiàng)式f^j,要求f^j滿足:

    PM-2j(f^j)=PM-2j[-(uj+qj)](26)

    f^j(xj,0)=-(j+j)|xj,0(27)

    f^j(xj,M+1)=-(j+j)|xj,M+1(28)

    第二步證明:

    f^jx(xj,n)=Fjx(xj,n),n∈{1,2,…,M}(29)

    具體證明過程如下:本文2.3節(jié)中已經(jīng)證明了DFR法和LDG法對于qj與pj的定義是等價(jià)的,故不再重復(fù)證明。首先令:

    f^j=-(uj+qj)+(-(j+j)|xj,0+

    (uj+qj)(xj,0))RMj,R+(-(j+j)|xj,M+1+

    (uj+qj)(xj,M+1))RMj,L

    其中,當(dāng)M≥1時,RMj,L是定義在Ej上的標(biāo)準(zhǔn)M次左拉登多項(xiàng)式,RMj,R是定義在Ej上的標(biāo)準(zhǔn)M次右拉登多項(xiàng)式。

    由左拉登多項(xiàng)式和右拉登多項(xiàng)式的正交性,RMj,L和RMj,R在Ej上任意不超過M-2次多項(xiàng)式構(gòu)成的空間PM-2j投影為0,可計(jì)算f^j在PM-2j的投影證得式(30)成立。由左拉登多項(xiàng)式的性質(zhì),RMj,L在點(diǎn)xj,0處取值為0,在點(diǎn)xj,M+1處取值為1,右拉登多項(xiàng)式的性質(zhì),RMj,R在xj,0處取值為1,在點(diǎn)xj,M+1處取值為0,可以計(jì)算f^j在點(diǎn)xj,0,xj,M+1處的取值證得式(27)和式(28)成立。

    對于n∈{1,2,…,M},觀察下式:

    Ejujtj,ndx+∫Ejf^jj,nxdx=

    f^j-j,n|xj+1/2f^j+j,n|xj-1/2(30)

    式(19)第一式和式(30)左側(cè)第一項(xiàng)相同。由式(26)可得:式(19)第一式和式(30)左側(cè)第二項(xiàng)相等。由式(27)和式(28)可得:式(19)第一式和式(30)右側(cè)兩項(xiàng)一一對應(yīng)相等。即式(19)第一式和式(30)等價(jià)。

    對(30)式左側(cè)第二項(xiàng)使用分部積分可得到:

    Ejujtj,ndx=∫Ejf^jxj,ndx(31)

    對式(31)左右兩側(cè)同時應(yīng)用M點(diǎn)高斯求積,再對兩側(cè)同時乘以wihj/2得到:

    (uj,n)t=f^jx(xj,n)(32)

    顯然式(32)與式(31)等價(jià)。然后令:

    QM+1j=Fj-f^j(33)

    由式(25)左側(cè)第二項(xiàng)和式(26)可得:

    PM-2j(QM+1j)=0(34)

    將QM+1j寫作勒讓德多項(xiàng)式的展開形式:

    QM+1j=∑M+1n=0CnLnj(35)

    其中,當(dāng)n≥0時,Lnj是定義在Ej上的標(biāo)準(zhǔn)n次勒讓德多項(xiàng)式,Cn是勒讓德多項(xiàng)式展開系數(shù)。由(34)式得:

    PM-2j(QM+1j)=PM-2j(∑M+1n=0CnLnj)=

    PM-2j(∑M-2n=1CnLnj)+PM-2j(∑M+1n=M-1CnLnj)=0

    由勒讓德多項(xiàng)式的正交性,Lnj在Ej上任意不超過M-1次多項(xiàng)式組成的空間PM-1j投影為0,可得到:

    C0=C1=…=CM-2=0

    由Fj和f^j左、右端點(diǎn)的定義可以得到:

    CM=0,CM-1=-CM+1

    于是式(33)可以被寫為如下形式:

    QM+1j=cLoM+1j(36)

    其中,當(dāng)M≥2時,LoM+1j是定義在Ej上的標(biāo)準(zhǔn)M+1次洛巴托多項(xiàng)式,c是一個常數(shù)。

    由洛巴托多項(xiàng)式的性質(zhì),M+1次洛巴托多項(xiàng)式導(dǎo)數(shù)的零點(diǎn)是M次勒讓德多項(xiàng)式的零點(diǎn),可以得到:

    f^jx(xj,n)-Fjx(xj,n)=f^j-Fj)x(xj,n)=

    (QM+1j)x(xj,n)=(cLoM+1j)x(xj,n)=0

    n∈{1,2,…,M}

    再由式(32)可得到:

    (uj,n)t=f^jx(xj,n)=Fjx(xj,n

    n∈{1,2,…,M}(37)

    式(37)說明:式(17)和式(32)是等價(jià)的。因?yàn)槭剑?9)第一式和式(30)是等價(jià)的,所以式(19)第一式和式(17)也是等價(jià)的。證畢。

    2.5 兩種算法用到的常微分方程組

    針對方程(12),以M=1為例說明DFR法和LDG法用于編程的常微分方程組是等價(jià)的。

    考慮在每個Ej上,使用DFR法處理方程(12)如下:

    uj=uj,1

    Fj=Fj,0^j,0+Fj,1^j,1+Fj,2^j,2

    ^j,0=(x-xj,1)(x-xj,2)(xj,0-xj,1)(xj,0-xj,2

    ^j,1=(x-xj,0)(x-xj,2)(xj,1-xj,0)(xj,1-xj,2

    ^j,2=(x-xj,0)(x-xj,1)(xj,2-xj,0)(xj,2-xj,1

    Fj,0=-(uj-1,1+qj,1

    Fj,1=-(uj,1+qj,1

    Fj,2=-(uj,1+qj+1,1

    因此:

    (Fj)x(xj,1)=Fj,0^j,0)x(xj,1)+

    Fj,1^j,1)x(xj,1)+Fj,2^j,2)x(xj,1)=

    (xj,1-xj,2)(xj,0-xj,1)(xj,0-xj,2)(-(uj-1,1+qj,1))+

    (xj,1-xj,0)(xj,2-xj,0)(xj,2-xj,1)(-(uj,1+qj+1,1))

    最終得到的常微分方程組為

    (uj,1)t=(xj,1-xj,2)(xj,0-xj,1)(xj,0-xj,2)(-(uj-1,1+qj,1))+(xj,1-xj,0)(xj,2-xj,0)(xj,2-xj,1)(-(uj,1+qj+1,1))

    qj,1=1(xj,2-xj,0)(pj+1,1-pj,1

    pj,1=1(xj,2-xj,0)(-(uj-1,1-uj,1))

    使用LDG法處理方程(16)最終得到的常微分方程組如下:

    (uj,1)t=1(xj,2-xj,0)(-(uj,1+qj+1,1))+

    1(xj,2-xj,0)(-(uj-1,1+qj,1))

    pj,1=1(xj,2-xj,0)(-(uj-1,1-uj,1))

    qj,1=1(xj,2-xj,0)(pj+1,1-pj,1

    由于M=1時,勒讓德多項(xiàng)式的零點(diǎn)xj,1=(xj,2+xj,0)/2。對比以上兩組常微分方程組可知,兩方程組相等。

    3 結(jié) 論

    本文研究了DFR法和LDG法求解高階偏微分方程的等價(jià)性問題。針對線性三階KdV方程,用兩種方法證明DFR法和LDG法所求得的數(shù)值解等價(jià),進(jìn)一步完善了求解偏微分方程時使用插值理論和投影理論的聯(lián)系。以上工作是針對線性方程所做出的研究,這使得研究過程中所用到的高斯求積結(jié)果具有高精度,但對于非線性方程這一點(diǎn)是難以滿足的。關(guān)于含有更高階偏導(dǎo)數(shù)的非線性方程的插值法與投影法是否具有等價(jià)性有待證明。

    參 考 文 獻(xiàn):

    [1] REED W H, HILL T R. Triangular Mesh Methods for the Neutron Transport Equation[R]. Los Alamos Report LA-UR-73-479: 1973.

    [2] BASSI F, REBAY S. A High-Order Accurate Discontinuous Finite Element Method for the Numerical Solution of the Compressible Navier-Stokes Equations[J].Journal of Computational Physics, 1997, 131(2): 267.

    [3] COCKBURN B, SHU C W. The Local Discontinuous Galerkin Method for Time-Dependent Convection Diffusion Systems[J]. SIAM Journal on Numerical Analysis, 1998, 35(6): 2440.

    [4] YAN J, SHU C W. A Local Discontinuous Galerkin Method for KdV Type Equations[J]. SIAM Journal on Numerical Analysis, 2002, 40(2): 769.

    [5] YAN J, SHU C W. Local Discontinuous Galerkin Methods for Partial Differential Equations with HigherOrder Derivatives[J]. Journal of Scientific Computing,2002, 17(1/4): 27.

    [6] SHU C W, XU Y. Local Discontinuous Galerkin Methods for Three Classes of Nonlinear Wave Equations[J]. Journal of Computational Mathematics, 2004,22 (2) :250.

    [7] XU Y, SHU C W. Local Discontinuous Galerkin Methods for Two Classes of Two-dimensional Nonlinear Wave Equations[J]. Physica D, 2005, 208 (1/2):21.

    [8] MENG X, SHU C W, WU B Y. Superconvergence of the Local Discontinuous Galerkin Method for Linear Fourth-Order Time-Dependent" Problems in One Space Dimension[J]. IMA Journal of Numerical Analysis, 2012, 32(4): 1294.

    [9] 曹外香,張智民. 解一維雙曲守恒律方程和拋物方程的間斷有限元法的逐點(diǎn)和區(qū)間平均值誤差估計(jì)[J]. 中國科學(xué):數(shù)學(xué) 2015, 45(8): 1115.

    CAO WaiXiang , ZHANG ZhiMin. Point-wise and Cell Average Error Estimates of the DG and LDG Methods for 1D Hyperbolic and Parabolic Equations[J].Scientia Sinica Mathematica,2015,45(8):1115.

    [10]CHENG Y, MENG X, ZHANG Q. Application of Generalized Gauss-Radau Projections for the Local Discontinuous Galerkin Method for Linear ConvectionDiffusion Equations[J]. Math. Comp, 2016, 86(305): 1233.

    [11]畢卉,陳莎莎. 四階線性方程局部間斷 Galerkin 方法的誤差估計(jì)[J]. 哈爾濱理工大學(xué)學(xué)報(bào), 2021,26(4): 159.

    BI Hui,CHEN Shasha.Error Estimates for Local Discontinuous Galerkin Methods for Linear Fourth-order Equations[J].Journal of Harbin University of Science and Technology,2021,26(4):159.

    [12]GUO W, HUANG J T, TAO Z J, et al. An Adaptive Sparse Grid Local Discontinuous Galerkin Method for Hamilton-Jacobi Equations in High Dimensions[J].Journal of Computational Physics 2021, 436: 170.

    [13]WEI L L, LI W B. Local Discontinuous Galerkin Approximations to Variable-Order Time-Fractional Diffusion Model Based on the Caputo-Fabrizio FractionalDerivative[J]. Mathematics and Computers in Simulation 2021, 188:280.

    [14]SHU C W. Discontinuous Galerkin Method for TimeDependent Problems: Survey and Recent Developments. In Recent Developments in Discontinuous Galerkin Finite Element Methods for Partial Differential Equations[M]. Lecture Notes in The IMA Volumes in Mathematics and its Applications, Cham: Springer, 2014, 157: 25.

    [15]SHU C W. Discontinuous Galerkin Methods for Time-Dependent Convection Dominated Problems: BasicsRecent Developments and Comparison with Other Methods[M]. In Building: Connections and Challenges in Modern Approaches to Numerical Partial Differential Equations. Springer International Publishing, 2016:369.

    [16]HUYNH H T. A Flux Reconstruction Approach to High-Order Schemes Including Discontinuous Galerkin Methods[C]. 18th AIAA Computational Fluid Dynamics Conference. 2007: 4079.

    [17]ROMERO J, ASTHANA K, JAMESON A. A Simplified Formulation of the Flux Reconstruction Method[J]. Journal of Scientific Computing, 2016, 67(1): 351.

    [18]ALEXANDRE E, MARTIN V. A Posteriori Error Estimation Based on Potential and Flux Reconstruction for the Heat Equation[J]. SIAM Journal on Numerical Analysis,2010, 48(1): 198.

    [19]WITHERDEN F D, VINCENT, P E, JAMESON, A.:High-order flux reconstruction schemes. In: Abgrall, R, SHU C W. (eds.) Handbook of Numerical Analysis, Chapter 10, vol. 17, pp. 227-263. Elsevier, Amsterdam (2016).

    [20]WANG Z J, HUYNH H T. : A Review of Flux Reconstruction or Correction Procedure Via Reconstruction Method for the Navier-Stokes Equations. Mech. Eng. Rev. 3(1), 1-16 (2016)

    [21]ROMERO J, ASTHANA K, JAMESON A. A Simplified Formulation of the Flux Reconstruction Method[J]. Journal of Scientific Computing, 2016, 67(1): 351.

    [22]ROMERO J, WITHERDEN F D, JAMESON A. A Direct Flux Reconstruction Scheme for Advection-Diffusion Problems on Triangular Grids[J]. Journal of Scientific Computing, 2017, 73(2/3): 1115.

    [23]YOU H J, KIM C. Direct Reconstruction Method for Discontinuous Galerkin Methods on Higher-Order Mixed-Curved Meshes I. Volume Integration[J]. Journal of Computational Physics, 2019,395: 223.

    [24]HUYNH H T. Discontinuous Galerkin Via Interpolation: The Direct Flux Reconstruction Method[J]. Journal of Scientific Computing 2020, 82(3): 28.

    [25]畢卉,劉磊. DFR法與DG法解拋物方程和對流擴(kuò)散方程等價(jià)性[J]. 哈爾濱理工大學(xué)學(xué)報(bào),2022(6):152.

    BI Hui,LIU Lei.Equivalence Between DFR Method and DG Method for Solving Parabolic Equation and Convection-diffusion Equation[J].Journal of Harbin University of Science and Technology,2022,27(6):152.

    (編輯:溫澤宇)

    基金項(xiàng)目: 國家自然科學(xué)基金青年基金(12201157);黑龍江省自然科學(xué)基金聯(lián)合引導(dǎo)項(xiàng)目(LH2020A015).

    作者簡介:李曉彤(1997—),女,碩士研究生.

    通信作者:畢 卉(1982—),女,博士,教授,博士研究生導(dǎo)師,E-mail:bihui@hrbust.edu.cn.

    猜你喜歡
    三階等價(jià)插值
    三階非線性微分方程周期解的非退化和存在唯一性
    基于Sinc插值與相關(guān)譜的縱橫波速度比掃描方法
    n次自然數(shù)冪和的一個等價(jià)無窮大
    中文信息(2017年12期)2018-01-27 08:22:58
    一種改進(jìn)FFT多譜線插值諧波分析方法
    基于四項(xiàng)最低旁瓣Nuttall窗的插值FFT諧波分析
    三類可降階的三階非線性微分方程
    收斂的非線性迭代數(shù)列xn+1=g(xn)的等價(jià)數(shù)列
    三階微分方程理論
    Blackman-Harris窗的插值FFT諧波分析與應(yīng)用
    環(huán)Fpm+uFpm+…+uk-1Fpm上常循環(huán)碼的等價(jià)性
    免费高清视频大片| 白带黄色成豆腐渣| 久久热在线av| 国产69精品久久久久777片 | 日韩欧美国产一区二区入口| 两个人的视频大全免费| 香蕉久久夜色| 看免费av毛片| 一级毛片女人18水好多| 国产精品电影一区二区三区| 欧美国产日韩亚洲一区| 国产激情久久老熟女| 老司机午夜福利在线观看视频| 可以在线观看的亚洲视频| av福利片在线| 麻豆国产av国片精品| 久久精品91无色码中文字幕| 久久久久久久精品吃奶| 久久久久久久精品吃奶| 欧美成人性av电影在线观看| av视频在线观看入口| 亚洲欧美一区二区三区黑人| 亚洲狠狠婷婷综合久久图片| 中文字幕人妻丝袜一区二区| 婷婷六月久久综合丁香| 午夜免费成人在线视频| 国内少妇人妻偷人精品xxx网站 | 制服丝袜大香蕉在线| 大型黄色视频在线免费观看| 999久久久精品免费观看国产| 此物有八面人人有两片| 久久精品亚洲精品国产色婷小说| 在线观看舔阴道视频| 精品久久久久久久人妻蜜臀av| 免费一级毛片在线播放高清视频| 久久久精品大字幕| 两个人看的免费小视频| 精品午夜福利视频在线观看一区| 女警被强在线播放| 香蕉国产在线看| a在线观看视频网站| 69av精品久久久久久| 在线观看美女被高潮喷水网站 | 色老头精品视频在线观看| 在线十欧美十亚洲十日本专区| 欧美乱色亚洲激情| 亚洲最大成人中文| 久久久久免费精品人妻一区二区| netflix在线观看网站| 在线观看舔阴道视频| 熟妇人妻久久中文字幕3abv| 欧美成人午夜精品| 欧美最黄视频在线播放免费| 正在播放国产对白刺激| 成人18禁在线播放| 国产又黄又爽又无遮挡在线| 亚洲欧美日韩高清在线视频| 一级毛片女人18水好多| 国产精品电影一区二区三区| 人人妻人人看人人澡| 搞女人的毛片| www日本在线高清视频| 国产熟女xx| 国产亚洲精品一区二区www| 国产精品乱码一区二三区的特点| 少妇熟女aⅴ在线视频| 又黄又爽又免费观看的视频| 免费无遮挡裸体视频| 一二三四在线观看免费中文在| 国产精品国产高清国产av| 亚洲国产精品合色在线| 黄色成人免费大全| 在线a可以看的网站| 精品不卡国产一区二区三区| 亚洲 欧美 日韩 在线 免费| 国产精品久久久久久久电影 | 国产成人精品久久二区二区91| 岛国在线免费视频观看| 亚洲一码二码三码区别大吗| 免费搜索国产男女视频| 午夜福利在线观看吧| 欧美在线一区亚洲| a级毛片在线看网站| 日韩av在线大香蕉| 久久久久九九精品影院| 成人手机av| 欧美在线一区亚洲| 欧美日韩中文字幕国产精品一区二区三区| 国产三级黄色录像| 长腿黑丝高跟| 女生性感内裤真人,穿戴方法视频| 国产精品九九99| 亚洲18禁久久av| 国产精品九九99| 天天添夜夜摸| 黄色视频,在线免费观看| 99国产精品一区二区蜜桃av| 一边摸一边抽搐一进一小说| 亚洲一卡2卡3卡4卡5卡精品中文| 久99久视频精品免费| 久久这里只有精品19| 国产97色在线日韩免费| 久久久精品大字幕| 亚洲一区二区三区不卡视频| 无限看片的www在线观看| 亚洲欧美日韩东京热| 国产午夜精品论理片| 精品国产亚洲在线| 日韩成人在线观看一区二区三区| 日韩精品免费视频一区二区三区| 老汉色av国产亚洲站长工具| 亚洲黑人精品在线| 成人手机av| 亚洲免费av在线视频| svipshipincom国产片| 精品久久久久久久久久久久久| 一级毛片高清免费大全| 免费无遮挡裸体视频| 亚洲成人久久爱视频| 琪琪午夜伦伦电影理论片6080| 亚洲av日韩精品久久久久久密| 精品国产乱码久久久久久男人| 国产精品久久久久久人妻精品电影| 99riav亚洲国产免费| 99在线视频只有这里精品首页| bbb黄色大片| 久久久久国产一级毛片高清牌| 观看免费一级毛片| 91在线观看av| 又紧又爽又黄一区二区| 亚洲av成人不卡在线观看播放网| 日韩成人在线观看一区二区三区| 美女黄网站色视频| 丝袜美腿诱惑在线| 亚洲aⅴ乱码一区二区在线播放 | 亚洲自偷自拍图片 自拍| 曰老女人黄片| 老熟妇仑乱视频hdxx| 一a级毛片在线观看| 啦啦啦韩国在线观看视频| 一本大道久久a久久精品| 在线观看www视频免费| 亚洲真实伦在线观看| 久久婷婷人人爽人人干人人爱| 日韩欧美三级三区| 久久久久久大精品| 禁无遮挡网站| 成人国语在线视频| 黄色a级毛片大全视频| 特级一级黄色大片| 女同久久另类99精品国产91| 色哟哟哟哟哟哟| 国产在线观看jvid| 国产亚洲精品第一综合不卡| 久久精品夜夜夜夜夜久久蜜豆 | 一本一本综合久久| 狂野欧美激情性xxxx| 国产精品久久电影中文字幕| 日本一本二区三区精品| 香蕉久久夜色| 国产精品亚洲一级av第二区| 大型黄色视频在线免费观看| 国产高清激情床上av| av在线天堂中文字幕| 青草久久国产| 很黄的视频免费| 校园春色视频在线观看| 真人做人爱边吃奶动态| 麻豆成人av在线观看| 亚洲成人久久性| av国产免费在线观看| 亚洲精品美女久久av网站| 1024视频免费在线观看| 欧美日韩一级在线毛片| 日韩欧美在线二视频| netflix在线观看网站| 在线a可以看的网站| 国产亚洲精品第一综合不卡| 深夜精品福利| 日日干狠狠操夜夜爽| 国产精品美女特级片免费视频播放器 | 少妇熟女aⅴ在线视频| 亚洲中文字幕一区二区三区有码在线看 | 天堂影院成人在线观看| 精品高清国产在线一区| 又紧又爽又黄一区二区| 2021天堂中文幕一二区在线观| 久久99热这里只有精品18| 一个人免费在线观看的高清视频| 久久久久九九精品影院| 国产日本99.免费观看| 国产aⅴ精品一区二区三区波| 日本熟妇午夜| 欧美一区二区精品小视频在线| 男女做爰动态图高潮gif福利片| 淫秽高清视频在线观看| 久久九九热精品免费| 人妻久久中文字幕网| 国产1区2区3区精品| 精品不卡国产一区二区三区| 村上凉子中文字幕在线| 美女高潮喷水抽搐中文字幕| 夜夜爽天天搞| 很黄的视频免费| 看片在线看免费视频| 999久久久精品免费观看国产| av欧美777| 欧美一区二区国产精品久久精品 | 正在播放国产对白刺激| 两性午夜刺激爽爽歪歪视频在线观看 | 日韩 欧美 亚洲 中文字幕| 久久久水蜜桃国产精品网| 母亲3免费完整高清在线观看| 成人国语在线视频| 真人做人爱边吃奶动态| 三级国产精品欧美在线观看 | 国产主播在线观看一区二区| 中文字幕熟女人妻在线| 久久九九热精品免费| 午夜久久久久精精品| 成年版毛片免费区| 五月伊人婷婷丁香| 此物有八面人人有两片| 免费在线观看影片大全网站| 99久久精品热视频| 久久精品国产亚洲av高清一级| 别揉我奶头~嗯~啊~动态视频| 精品国产乱子伦一区二区三区| 国产精品98久久久久久宅男小说| 午夜福利欧美成人| 亚洲国产看品久久| 日日夜夜操网爽| 国产精品国产高清国产av| 亚洲精品在线观看二区| 成人国产一区最新在线观看| 亚洲成人免费电影在线观看| 一边摸一边做爽爽视频免费| 久久国产乱子伦精品免费另类| 19禁男女啪啪无遮挡网站| 在线观看免费日韩欧美大片| 成人亚洲精品av一区二区| 亚洲专区国产一区二区| 日本精品一区二区三区蜜桃| 18美女黄网站色大片免费观看| 日韩欧美国产在线观看| 国产视频一区二区在线看| 香蕉av资源在线| 色播亚洲综合网| 成人三级做爰电影| 两性午夜刺激爽爽歪歪视频在线观看 | 久久99热这里只有精品18| 欧美成人性av电影在线观看| 欧美国产日韩亚洲一区| 国产乱人伦免费视频| 美女午夜性视频免费| 亚洲av第一区精品v没综合| 久久草成人影院| 中文字幕高清在线视频| 我的老师免费观看完整版| 日韩欧美一区二区三区在线观看| 国产成人精品久久二区二区91| 性欧美人与动物交配| 欧美高清成人免费视频www| 精品电影一区二区在线| 99在线人妻在线中文字幕| 18禁黄网站禁片午夜丰满| 一级毛片女人18水好多| 亚洲国产高清在线一区二区三| 宅男免费午夜| 一边摸一边抽搐一进一小说| 淫妇啪啪啪对白视频| 欧美三级亚洲精品| 国产精品,欧美在线| 高清毛片免费观看视频网站| 国产在线精品亚洲第一网站| 少妇粗大呻吟视频| 亚洲av成人av| 色精品久久人妻99蜜桃| 最新美女视频免费是黄的| 国产精品亚洲一级av第二区| 国产私拍福利视频在线观看| 国产精品久久久久久人妻精品电影| 在线播放国产精品三级| 麻豆成人av在线观看| 宅男免费午夜| 欧美色视频一区免费| 国产成人精品久久二区二区免费| 精品电影一区二区在线| 视频区欧美日本亚洲| 99热6这里只有精品| 成人精品一区二区免费| 亚洲全国av大片| 91成年电影在线观看| aaaaa片日本免费| 国产成+人综合+亚洲专区| 男女午夜视频在线观看| 国产伦人伦偷精品视频| 特大巨黑吊av在线直播| 亚洲最大成人中文| 少妇人妻一区二区三区视频| 久久精品国产亚洲av高清一级| 免费看日本二区| 亚洲va日本ⅴa欧美va伊人久久| 亚洲欧美日韩无卡精品| 欧美zozozo另类| xxx96com| 国产av在哪里看| 亚洲精品中文字幕一二三四区| 日韩欧美在线二视频| 成人特级黄色片久久久久久久| 国产欧美日韩精品亚洲av| 国产成人av教育| 天天一区二区日本电影三级| 欧美黑人精品巨大| 青草久久国产| 亚洲av成人一区二区三| 精华霜和精华液先用哪个| 亚洲aⅴ乱码一区二区在线播放 | 99国产极品粉嫩在线观看| 叶爱在线成人免费视频播放| 婷婷丁香在线五月| 精品不卡国产一区二区三区| АⅤ资源中文在线天堂| 久久久久久亚洲精品国产蜜桃av| 国产亚洲精品一区二区www| 丁香欧美五月| 丝袜人妻中文字幕| 久久 成人 亚洲| 无遮挡黄片免费观看| 中亚洲国语对白在线视频| 欧美一级a爱片免费观看看 | 国产成人aa在线观看| 亚洲自偷自拍图片 自拍| 久久精品夜夜夜夜夜久久蜜豆 | 最新美女视频免费是黄的| 老熟妇仑乱视频hdxx| 日韩欧美免费精品| 亚洲第一电影网av| 中文亚洲av片在线观看爽| 亚洲精品粉嫩美女一区| e午夜精品久久久久久久| 亚洲最大成人中文| 18禁国产床啪视频网站| 日韩有码中文字幕| 一本大道久久a久久精品| 精品国产亚洲在线| 一a级毛片在线观看| 最近在线观看免费完整版| 国产精品亚洲美女久久久| 国产蜜桃级精品一区二区三区| 香蕉国产在线看| 国产亚洲精品一区二区www| 丰满的人妻完整版| 亚洲成人久久爱视频| 激情在线观看视频在线高清| 免费在线观看完整版高清| 午夜老司机福利片| 国产精品国产高清国产av| 美女黄网站色视频| bbb黄色大片| 老熟妇乱子伦视频在线观看| 国产精品一区二区三区四区久久| 99热6这里只有精品| 在线观看66精品国产| 亚洲熟妇中文字幕五十中出| 精品一区二区三区四区五区乱码| 又黄又粗又硬又大视频| 男女做爰动态图高潮gif福利片| 制服诱惑二区| 香蕉久久夜色| 黄色女人牲交| 亚洲成人久久爱视频| 欧美乱色亚洲激情| 亚洲美女视频黄频| 99在线视频只有这里精品首页| 一本综合久久免费| 久久精品国产亚洲av香蕉五月| bbb黄色大片| 777久久人妻少妇嫩草av网站| videosex国产| 精品久久久久久成人av| 三级国产精品欧美在线观看 | 国产精品日韩av在线免费观看| 又爽又黄无遮挡网站| 小说图片视频综合网站| АⅤ资源中文在线天堂| 久久久国产成人精品二区| 欧美成人午夜精品| 老汉色av国产亚洲站长工具| aaaaa片日本免费| 国产精品一区二区精品视频观看| 日韩欧美一区二区三区在线观看| 1024手机看黄色片| 欧美 亚洲 国产 日韩一| 亚洲男人天堂网一区| 婷婷丁香在线五月| 亚洲最大成人中文| 美女黄网站色视频| 久久久久性生活片| 欧美成人午夜精品| 免费观看人在逋| 久久久精品大字幕| 少妇人妻一区二区三区视频| 男插女下体视频免费在线播放| 啦啦啦韩国在线观看视频| 亚洲国产精品sss在线观看| 亚洲成av人片免费观看| 91在线观看av| 久久亚洲精品不卡| 欧美日韩一级在线毛片| 亚洲专区字幕在线| 亚洲av成人av| 国产69精品久久久久777片 | 不卡一级毛片| 国产黄色小视频在线观看| 中文在线观看免费www的网站 | 亚洲人成伊人成综合网2020| 国产精品1区2区在线观看.| 舔av片在线| 久久国产精品人妻蜜桃| 色在线成人网| 真人一进一出gif抽搐免费| 国产亚洲精品一区二区www| 草草在线视频免费看| 一进一出抽搐动态| 麻豆成人av在线观看| 国产一区二区三区在线臀色熟女| 黄色丝袜av网址大全| 国产男靠女视频免费网站| 一边摸一边抽搐一进一小说| 在线观看舔阴道视频| 国产欧美日韩精品亚洲av| 国产激情久久老熟女| 波多野结衣高清无吗| 1024视频免费在线观看| 中文字幕高清在线视频| 一进一出好大好爽视频| 成人国语在线视频| 麻豆久久精品国产亚洲av| 国产精品 国内视频| 精品久久久久久久人妻蜜臀av| 国产激情欧美一区二区| aaaaa片日本免费| 国产一区二区三区视频了| 日韩av在线大香蕉| 成人欧美大片| 国产私拍福利视频在线观看| 欧美3d第一页| 国产精品综合久久久久久久免费| 久久亚洲真实| 久久精品亚洲精品国产色婷小说| 国产av不卡久久| 欧美国产日韩亚洲一区| ponron亚洲| www.www免费av| 欧美性长视频在线观看| 十八禁网站免费在线| 色尼玛亚洲综合影院| 美女大奶头视频| 18禁黄网站禁片午夜丰满| 亚洲人成网站在线播放欧美日韩| 免费看美女性在线毛片视频| 黄色毛片三级朝国网站| 国产精品,欧美在线| 精品欧美一区二区三区在线| 久久精品人妻少妇| 在线观看午夜福利视频| 日韩国内少妇激情av| 露出奶头的视频| 日本一本二区三区精品| 黄色女人牲交| 色av中文字幕| 亚洲国产精品999在线| 三级毛片av免费| 黄色 视频免费看| 丝袜人妻中文字幕| 又黄又爽又免费观看的视频| 国产一区二区三区在线臀色熟女| 亚洲在线自拍视频| 久久国产乱子伦精品免费另类| 三级毛片av免费| 视频区欧美日本亚洲| 国产人伦9x9x在线观看| 国产不卡一卡二| 久久久久久亚洲精品国产蜜桃av| ponron亚洲| 禁无遮挡网站| av免费在线观看网站| 久久精品国产亚洲av高清一级| 亚洲国产欧美网| 久久久水蜜桃国产精品网| 欧美大码av| 久久99热这里只有精品18| 亚洲av电影不卡..在线观看| 此物有八面人人有两片| 午夜精品久久久久久毛片777| 亚洲精品在线美女| 这个男人来自地球电影免费观看| 欧美日韩精品网址| 精品国产乱码久久久久久男人| 亚洲精品av麻豆狂野| 国产私拍福利视频在线观看| 国产单亲对白刺激| 日韩精品青青久久久久久| 久久香蕉激情| 精品一区二区三区视频在线观看免费| 亚洲精华国产精华精| 免费无遮挡裸体视频| 两个人视频免费观看高清| 一本大道久久a久久精品| 天天添夜夜摸| 一区二区三区国产精品乱码| 天堂动漫精品| 亚洲人成网站高清观看| 国产精品av视频在线免费观看| 欧美一级毛片孕妇| 国产三级在线视频| aaaaa片日本免费| 男人舔奶头视频| 国产真实乱freesex| 在线观看www视频免费| 国产精品久久视频播放| 老司机深夜福利视频在线观看| 国产精品av久久久久免费| 亚洲欧美日韩高清在线视频| 丰满的人妻完整版| 亚洲av电影不卡..在线观看| 日韩成人在线观看一区二区三区| 亚洲人成伊人成综合网2020| 国产三级中文精品| 女警被强在线播放| 一进一出抽搐gif免费好疼| 久9热在线精品视频| 欧美成人午夜精品| 在线观看午夜福利视频| 亚洲专区国产一区二区| 国模一区二区三区四区视频 | 757午夜福利合集在线观看| 国内毛片毛片毛片毛片毛片| 国产伦一二天堂av在线观看| or卡值多少钱| 欧美乱码精品一区二区三区| 欧美日韩亚洲综合一区二区三区_| 国产激情久久老熟女| 国产一区在线观看成人免费| 国产v大片淫在线免费观看| 搡老岳熟女国产| 亚洲成a人片在线一区二区| 最好的美女福利视频网| aaaaa片日本免费| 色在线成人网| 老司机午夜福利在线观看视频| 成人三级黄色视频| 亚洲第一欧美日韩一区二区三区| 国产在线观看jvid| 久久草成人影院| 天天躁狠狠躁夜夜躁狠狠躁| 国产成人欧美在线观看| 一级a爱片免费观看的视频| 国产v大片淫在线免费观看| 一级a爱片免费观看的视频| 欧美日韩亚洲国产一区二区在线观看| 18禁黄网站禁片免费观看直播| 久久草成人影院| 日韩 欧美 亚洲 中文字幕| 精品午夜福利视频在线观看一区| 亚洲成人精品中文字幕电影| 婷婷丁香在线五月| av在线播放免费不卡| 高清在线国产一区| 欧美黄色片欧美黄色片| 精品国内亚洲2022精品成人| 午夜福利免费观看在线| 欧美最黄视频在线播放免费| 黄色视频不卡| 中文字幕av在线有码专区| 九九热线精品视视频播放| 欧美日本视频| 免费观看人在逋| 色老头精品视频在线观看| 最近最新免费中文字幕在线| 熟女电影av网| 99riav亚洲国产免费| 午夜两性在线视频| 国产av又大| 国产精品亚洲一级av第二区| 成人手机av| 好男人电影高清在线观看| 欧美色欧美亚洲另类二区| 国产精品久久久久久亚洲av鲁大| 夜夜夜夜夜久久久久| 亚洲国产精品sss在线观看| 校园春色视频在线观看| 亚洲激情在线av| 亚洲欧美日韩高清在线视频| 琪琪午夜伦伦电影理论片6080| 人人妻人人看人人澡| 搡老妇女老女人老熟妇| 国产黄片美女视频| www.www免费av| 两个人的视频大全免费| 国内精品久久久久久久电影| 丁香欧美五月| 精品久久久久久成人av| 国产成人系列免费观看| 精品午夜福利视频在线观看一区| 久久亚洲精品不卡| 色尼玛亚洲综合影院| 日韩欧美三级三区| 久久午夜综合久久蜜桃| 日韩三级视频一区二区三区| 国产熟女xx| 亚洲黑人精品在线| 亚洲熟女毛片儿| av福利片在线|