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

    一類具有時(shí)滯的Leslie-Gower捕食-食餌模型的Hopf分支

    2024-01-01 00:00:00袁海龍樊雨李一多
    關(guān)鍵詞:時(shí)滯穩(wěn)定性

    摘要: 利用Hopf分支理論, 研究一類具有時(shí)滯的Leslie-Gower捕食-食餌模型. 首先, 以時(shí)滯為分支參數(shù), 討論該模型正平衡點(diǎn)的穩(wěn)定性和Hopf分支的存在性; 其次, 根據(jù)偏泛函微分方程的規(guī)范型理論和中心流形定理, 確定Hopf分支的分支方向和分支周期解的穩(wěn)定性; 最后, 利用MATLAB進(jìn)行數(shù)值模擬.

    關(guān)鍵詞: 時(shí)滯; Leslie-Gower模型; Hopf分支; 穩(wěn)定性

    中圖分類號: O175.12" 文獻(xiàn)標(biāo)志碼: A" 文章編號: 1671-5489(2024)04-0821-10

    Hopf Bifurcation of a Class of Leslie-GowerPredator-Prey Models with Time Delay

    YUAN Hailong, FAN Yu, LI Yiduo

    (School of Mathematics amp; Data Science, Shaanxi University of Science amp; Technology, Xi’an 710021, China)

    Abstract: Using the Hopf bifurcation theory, we studied a class of Leslie-Gower predator-prey models with time delay.

    Firstly, taking time delay as the bifurcation parameter, we discussed the stability of the positive equilibrium point of the model and the existence of Hopf bifurcation.

    Secondly, according to the normal form theory and center manifold theorem for partial differential equation, we derived the direction of Hopf bifurcation and the stability of

    bifurcation periodic solutions. Finally, we used MATLAB for numerical simulations.

    Keywords: time delay; Leslie-Gower model; Hopf bifurcation; stability

    0 引 言

    捕食-食餌關(guān)系是自然界種群之間最重要和最廣泛的基本關(guān)系之一, 種群模型在研究外來物種入侵、 流行疾病傳播以及自催化化學(xué)反應(yīng)等方面都具有重要作用. 目前, 對捕食-食餌關(guān)系模型的研究已取得了豐富的成果[1-5].

    考慮到Leslie-Gower捕食-食餌模型中食餌的生長速率可能具有Allee效應(yīng), Ni等[6]研究了如下捕食-食餌模型:

    dudt=u(1-u)ub-1-βuv,tgt;0,

    dvdt=μv1-vu,tgt;0,u(0)=u0gt;0, v(0)=v0gt;0,(1)

    其中u和v分別表示食餌和捕食者的種群密度, b∈(0,1)為Allee效應(yīng)的閾值, 初始值u0和v0為非負(fù)的連續(xù)函數(shù), 參數(shù)β和μ均為正常數(shù). 對于系統(tǒng)(1), Ni等[6]主要討論了當(dāng)系統(tǒng)只有一個(gè)正平衡點(diǎn)時(shí), 正平衡點(diǎn)的穩(wěn)定流形和不穩(wěn)定流形.

    由于種群具有成熟期, 當(dāng)前的種群數(shù)量會依賴于過去某一時(shí)刻的種群狀態(tài). 因此考慮時(shí)滯影響可更準(zhǔn)確地描述系統(tǒng)的動(dòng)力學(xué)特性. 時(shí)滯效應(yīng)在自然界中普遍存在, 時(shí)滯微分方程的穩(wěn)定性及其分支問題也備受關(guān)注. 例如: Yang等[7]考慮一類具有恒定獵物避難所和時(shí)滯的捕食-食餌系統(tǒng), 研究了其正平衡點(diǎn)的局部穩(wěn)定性和圖靈不穩(wěn)定性, 證明了時(shí)滯對系統(tǒng)的影響, 得到了Hopf分支的存在性并確定了Hopf分支的性質(zhì), 同時(shí)做出了數(shù)值模擬; Chen等[8]考慮一個(gè)具有齊次Neumann邊界條件的時(shí)滯擴(kuò)散Leslie-Gower捕食-食餌系統(tǒng), 通過分析特征方程, 研究了其共存平衡點(diǎn)的穩(wěn)定性和不穩(wěn)定性以及相關(guān)的Hopf分支, 利用上下解方法, 給出了其共存平衡點(diǎn)全局漸近穩(wěn)定的參數(shù)充分條件; Du等[9]考慮具有雙時(shí)滯和擴(kuò)散的修正Leslie-Gower捕食-食餌模型, 研究了其正平衡點(diǎn)的穩(wěn)定性以及Hopf、 雙Hopf分支的存在性, 并給出了雙Hopf分支點(diǎn)附近中心流形上的正規(guī)形式及臨界點(diǎn)附近的展開, 還得到了雙Hopf分支點(diǎn)附近的復(fù)雜動(dòng)力學(xué)行為.

    常笑源等[10]在系統(tǒng)(1)的基礎(chǔ)上, 考慮如下捕食-食餌系統(tǒng):

    dudt=ru1-uK(u-m)-quv,tgt;0,

    dvdt=sv1-v(t-τ)u(t-τ),tgt;0,

    u(θ)=1(θ), v(θ)=2(θ),θ∈[-τ,0],(2)

    其中時(shí)滯τ表示捕食者的成熟期, K為食餌的最大環(huán)境容納量, r和s分別為食餌和捕食者的自然增長率, q表示食餌人均減少率的最大值, m∈(0,1)為

    Allee效應(yīng)的閾值. 對于系統(tǒng)(2), 文獻(xiàn)[10]給出了其非負(fù)平衡點(diǎn)的存在性和穩(wěn)定性, 并以時(shí)滯τ作為分支參數(shù)證明了在正平衡點(diǎn)附近產(chǎn)生Hopf分支.

    基于以上研究, 本文考慮食餌的成熟期對捕食者種群密度的影響, 在系統(tǒng)(1)中對捕食者方程引入時(shí)滯, 同時(shí)考慮其包含齊次Neumann邊界條件的反應(yīng)擴(kuò)散系統(tǒng):

    ut-d1Δu=u(1-u)ub-1-βuv,x∈Ω, tgt;0,

    vt-d2Δv=μv1-vu(t-τ),x∈Ω, tgt;0,

    νu=νv=0,x∈Ω, tgt;0,

    u(x,0)=u0(x)gt;0, v(x,0)=v0(x)≥0(不恒為0),x∈,(3)

    其中d1和d2分別對應(yīng)食餌和捕食者的擴(kuò)散系數(shù), Δ為Laplace算子, Ω∈瘙綆N為具有光滑邊界Ω的有界開集, ν為邊界Ω上的單位外法向量.

    為書寫方便, 做無量綱變換, 令

    t=tb, m=b, q=bβ, s=μ,

    仍用t表示t, 在空間域Ω∈[0,lπ]考慮系統(tǒng)(3), 且Δ=2x2, l∈瘙綆+, 從而系統(tǒng)(3)可化為

    ut-d12ux2=u(1-u)(u-m)-quv,x∈(0,lπ), tgt;0,

    vt-d22vx2=sv1-vu(t-τ),x∈(0,lπ), tgt;0,

    u(x,t)x=v(x,t)x=0,x=0,lπ, tgt;0,

    u(x,t)=u0(x,t)gt;0, v(x,t)=v0(x,t)≥0(不恒為0),x∈[0,lπ], t∈[-τ,0].(4)

    本文研究時(shí)滯效應(yīng)在偏微分系統(tǒng)中對正平衡點(diǎn)產(chǎn)生的影響, 并分析正平衡點(diǎn)的穩(wěn)定性和Hopf分支的存在性, 給出判斷Hopf分支的分支方向和分支周期解穩(wěn)定性的表達(dá)式. 結(jié)果表明:

    當(dāng)時(shí)滯參數(shù)小于某一臨界值時(shí), 正平衡點(diǎn)的穩(wěn)定性不發(fā)生改變; 當(dāng)時(shí)滯參數(shù)增大且經(jīng)過某一臨界值時(shí), 會改變正平衡點(diǎn)的穩(wěn)定性, 系統(tǒng)在正平衡點(diǎn)附近會發(fā)生振蕩并產(chǎn)生Hopf分支.

    用瘙綃和瘙綆+分別表示非負(fù)整數(shù)集和正實(shí)數(shù)集.

    1 平衡點(diǎn)的穩(wěn)定性和Hopf分支的存在性

    考慮具有時(shí)滯參數(shù)的偏微分系統(tǒng)(4), 根據(jù)計(jì)算及文獻(xiàn)[10]可知:

    1) 系統(tǒng)(4)存在兩個(gè)邊界平衡點(diǎn)E0(1,0)和Em(m,0);

    2) 當(dāng)0lt;qlt;(1-m)2時(shí), 系統(tǒng)(4)存在兩個(gè)正平衡點(diǎn)E1(u*,v*)和E2(u,v), 其中

    u*=v*=(1+m-q)+(1+m-q)2-4m2,u=v=(1+m-q)-(1+m-q)2-4m2;

    3) 當(dāng)q=(1-m)2時(shí), 系統(tǒng)(4)存在唯一一個(gè)正平衡點(diǎn)E3(m,m);

    4) 當(dāng)qgt;(1-m)2時(shí), 系統(tǒng)(4)不存在正平衡點(diǎn).

    由于當(dāng)0lt;qlt;(1-m)2時(shí), 正平衡點(diǎn)E2(u,v)是不穩(wěn)定的, 因此本文首先討論當(dāng)0lt;qlt;(1-m)2時(shí)系統(tǒng)(4)正平衡點(diǎn)E1(u*

    ,v*)的穩(wěn)定性及Hopf分支. 當(dāng)τ=0時(shí), 滿足條件0lt;qlt;(1-m)2, sgt;s*=-2u*2+(m+1)u*, 且u*∈m+12,1, 其中0

    lt;mlt;1, 則E1(u*,v*)是局部漸近穩(wěn)定的. 然后分析當(dāng)τ≠0時(shí), 對系統(tǒng)(4)的正平衡點(diǎn)E1(u*,v*)的影響, 若τ∈[0,τ00), 則系統(tǒng)(4)的正平衡點(diǎn)E1(u*,v*)

    是局部漸近穩(wěn)定的; 若τ=τjn, j∈瘙綃, 則存在n使系統(tǒng)(4)在正平衡點(diǎn)E1(u*,v*)處產(chǎn)生Hopf分支; 若τgt;τ00, 則系統(tǒng)(4)的正平衡點(diǎn)E1(u*,v*)是不穩(wěn)定的.

    定義X=C([0,lπ],瘙綆2), 在系統(tǒng)(4)中令=u-u*, =v-v*, 仍用u,v表示,, 則其在抽象空間C([-τ,0],X)中具有如下形式的抽象泛函微分方程:

    dU(t)dt=dΔU

    (t)+L(Ut)+F(Ut),(5)

    其中dΔ=(d1Δ,d2Δ), dom(dΔ)={(u,v)T: u,v∈C2([0,lπ],瘙綆), ux,vx=0, x=0,lπ}.

    已知L: C([-τ,0],X)→X, F: C([-τ,0],

    X)→X, 對于Ut=φ, U(t

    )=φ(0), φ=(φ1,φ2 )T∈C([-τ,0]), 有

    L(φ)=[-2u*2+(m+1)u*]φ1(0)-qu*φ2(0)

    sφ1(-τ)-sφ2(0)," F(φ)=F1(φ)

    F2(φ),

    其中

    F1(φ)= "(φ1(0)+u*)[1-(φ1(0)+u*)](φ1(0)+u*-m)- "q(φ1(0)+u*)(φ2(0)+v*)-[-2u*2+(m+1)u*]φ1(0)+qu*φ2(0),

    F2(φ)= "s(φ2(0)+v*)1-φ2(0)+v*(φ1(-τ)+u*)-sφ1(-τ)+sφ2(0).

    系統(tǒng)(4)在(0,0)附近的線性化系統(tǒng)為

    dU(t)dt=dΔU

    (t)+L(Ut).(6)

    根據(jù)文獻(xiàn)[11], 可得線性系統(tǒng)的特征方程為

    λy-dΔy-L(eλ·y)=0, y∈

    dom(dΔ), y≠0.(7)

    由特征值問題-ψ″=μψ, x∈(0,lπ), ψ′(0)=ψ′(lπ)=0, 可得特征值μn=n2l2(n∈瘙綃), 相應(yīng)的特

    征函數(shù)為ψn(x)=cosnlx, 將

    y=∑∞n=0cosnlxy1ny2n

    代入特征方程(7)得

    -2u*2+(m+1)u*-d1n2l2-qu*

    se-λτ-s-d2n2l2y1n

    y2n=λy1ny2n," n∈瘙綃.

    因此, 特征方程(7)的所有特征根由以下特征方程給出:

    Δn(λ,τ)=λ2+Anλ+Be-λτ+Cn=0," n∈瘙綃,(8)

    其中

    An=(d1+d2)n2l2+2u*2-(m+1)u*+s," B=qu*s,

    Cn=d1n2l2+2u*2-(m+1)u*d2n2 l2+s.

    若λ=±iσ(σgt;0)為方程(8)的一對純虛根, 則有

    σ2-Cn=Bcos(στ),σAn=Bsin(στ)," n∈瘙綃,(9)

    將式(9)中兩式平方后再相加, 可得

    σ4+(A2n-2Cn)σ2+C2n-B2=0," n∈瘙綃,(10)

    其中

    A2n-2Cn=d1n2l2+2u*2-(m+1)u*2+d2n2l2+s2,

    C2n-B2=d1n2l2+2u*2-(m+1)u*2d2n2l2+s2-q2u*2s2.(11)

    當(dāng)m+qlt;1時(shí), 有C20-B2lt;0成立, 由于 limn→∞(C2n-B2)=+∞, 故存在一個(gè)最小的N0, 使得當(dāng)0≤n≤N0時(shí), 方程(10)至多有一個(gè)正根,

    當(dāng)ngt;N0時(shí), 方程(10)沒有正根. 當(dāng)0≤n≤N0時(shí), 方程(10)有一個(gè)正根σn, 滿足

    σn=-(A2n-2Cn)+(A2n-2Cn)2-4(C2n-B2)2,(12)

    則τ的表達(dá)式可確定為

    τ=τjn=τ0n+2jπσn," j∈瘙綃,(13)

    其中

    τ0n=1σnarccosσ2n-CnB,An≥0,

    1σn2π-arccosσ2n-CnB,Anlt;0.(14)

    此時(shí)方程(8)有一對純虛根±iσn.

    令λn(τ)=αn(τ)+iβn(τ)是方程(8)的根, 且當(dāng)τ→τjn時(shí), 滿足αn(τjn)=0和βn(τjn)=σn, 下面給出橫截條件.

    引理1 對于0≤n≤N0, j∈瘙綃, 且sgt;s*, 有α′n(τjn)gt;0.

    證明: 將λn(τ)代入方程(8)并關(guān)于τ求導(dǎo), 得

    dαn dττ=τjn-1=Re

    2λneλnτ+Aneλnτ-BτBλnτ=τjn,

    由于σn和τjn滿足方程(9), 并且根據(jù)方程(12)中σn的表達(dá)式, 可得

    dαndττ=τjn-1= "Re(2λn+An)(cos(λnτ)+isin (λnτ))-BτBλn=

    (A2n-2Cn)2-4(C2n-B2)B2gt;0.

    因此橫截條件成立, 證畢.

    根據(jù)式(13), 顯然τj+1ngt;τjn, 下面給出τjn關(guān)于n的單調(diào)性.

    性質(zhì)1 若0≤n≤N0, j∈瘙綃, 且sgt;s*, 則有τjn+1gt;τjn成立.

    證明: 將式(12)變形可得

    σ2n=-(A2n-2Cn)+(A2n-2Cn)2-4(C2n-B2)2=2(A2n-2Cn)2(B2-C

    2n)2+4B2-C2n+A2n-2CnB2-C2n,

    其中A2n-2Cn和C2n-B2都在式(11)中給出. 又由于sgt;s*, 當(dāng)0≤n≤N0時(shí), A2n-2Cn關(guān)于n嚴(yán)格遞增, 同時(shí)B2-C2n關(guān)于n嚴(yán)格遞減, 所以可得

    σ2n+1lt;σ2n. 由于sgt;s*時(shí), An≥0, 根據(jù)式(14), τ0n=1σnarccosσ2n-CnB. 因此, 當(dāng)0≤n≤N0 時(shí), 有τ0n+1gt;τ0n.

    因?yàn)棣襫+1lt;σn, 結(jié)合式(13), 所以當(dāng)j≥1, 0≤n≤N0時(shí), τjn+1gt;τjn.

    綜合上述分析, 可得如下定理:

    定理1 設(shè)0lt;qlt;(1-m)2, sgt;s*, u*∈m+12,1成立, 則有:

    1) 若τ∈[0,τ00), 則系統(tǒng)(4)的正平衡點(diǎn)E1(u*,v*)是局部漸近穩(wěn)定的;

    2) 若τ=τjn, j∈瘙綃, 0≤n≤N0, 則系統(tǒng)(4)在正平衡點(diǎn)E1(u*,v*)處產(chǎn)生Hopf分支;

    3) 若τgt;τ00, 則系統(tǒng)(4)的正平衡點(diǎn)E1(u*,v*)是不穩(wěn)定的.

    2 Hopf分支的方向和穩(wěn)定性

    當(dāng)τ=τjn時(shí), 系統(tǒng)(4)在E1(u*,v*)處產(chǎn)生Hopf分支. 下面利用時(shí)滯偏微分方程的中心流形定

    理[12]和規(guī)范型理論[11,13]研究τ=τ0=τ00時(shí)Hopf分支的方向和分支周期解的穩(wěn)定性.

    令τ=τ0+μ, 則μ=0為系統(tǒng)(5)的Hopf分支值. 令t=tτ并代入方程(5), 且仍用t表示t

    , 同理, Ut=φ, U(t)=φ(0), 則可將系統(tǒng)(5)轉(zhuǎn)化為

    dU(t)dt=τ0dΔ

    U(t)+τ0L0(Ut)+G(Ut,μ),(15)

    其中,

    L0(φ)=[-2u*2+(m+1)u*]φ1(0)-qu*φ2(0)sφ1(-1)-sφ2(0),

    G(φ,μ)=μdΔφ(0)+μL0(φ)+(μ+τ

    0)F0(φ),F(xiàn)0(φ)=(φ1(0)+u*)([1-(φ1(0)+u*)](φ1(0)+u*-m)-q(φ2(0)+v*))

    -[-2u*2+(m+1)u*]φ1(0)+qu*φ2(0)s(φ2(0)+v*)1-φ2(0)+v*(φ1(-τ)+u*)-sφ1(-1)+sφ2(0),

    且φ∈C([-1,0],X).

    系統(tǒng)(5)在(0,0)處的線性化系統(tǒng)為

    dU(t)dt=τ0dΔU(t)+τ0L0(Ut).(16)

    由式(10)可知, ±iσ0τ0為系統(tǒng)(16)的一對純虛根, 其線性泛函微分方程為

    dz(t)dt=τ0L0(zt).(17)

    根據(jù)Riesz表示定理可知, 存在η(θ,μ)(θ∈[-1,0])為2×2階的有界變差函數(shù)矩陣, 且滿足

    (τ0+μ)L0(φ)=∫0-1d

    η(θ,μ)φ(θ)," φ(θ)∈C([-1,0],瘙綆2),

    其中η(θ,μ)=(τ0+μ)Eδ(θ)-(τ0+μ)Fδ(θ+1),

    E=-2u*2+(m+1)u*-qu*0-s," F=00s0,

    且對于δ(θ): [-1,0]→(X,X), 有δ(θ)=0,θ∈[-1,0),1,θ=0.

    若φ(θ)∈C1([-1,0],瘙綆2), 則A(0)定義為

    A(0)(φ(θ))=dφ(θ)dθ,θ∈[-1,0),

    ∫0-1dη(θ,0)φ(θ),θ=0.

    若ψ=(ψ1,ψ2)∈C1([-1,0],(瘙綆2)*), 則A*定義為

    A*(ψ(s))=-dψ(s)ds,s∈(0,1],

    ∫0-1ψ(-s)dη(θ,0),s=0.

    若φ(θ)∈C1([-1,0],瘙綆2), ψ(s)∈C1([-1,0],(瘙綆2)*), 引入如下雙線性形式:

    〈ψ(s),φ(θ)〉0=ψ(0)φ(0)-∫0-1∫θ0

    ψ(ξ-θ)dη(θ,0)φ(ξ)dξ,

    則A(0)和A*是雙線性形式〈ψ(s),φ(θ)〉0下的伴隨算子.

    經(jīng)驗(yàn)證可知, ±iσ0τ0是A(0)和A*的特征值, 設(shè)q(θ)=(q1,q2)Teiσ0τ0

    θ(θ∈[-1,0])和q*(s)=1(q*1,q*2)eiσ0τ0s(s∈[0,1])分

    別為A(0)和A*關(guān)于特征值iσ0τ0和-iσ0τ0對應(yīng)的特征向量, 其中

    (q1,q2)=1,-2u*+m+1q-iσ0qu

    *," (q*1,q*2)=1,-qu*(s+iσ0)s2+σ20,

    D=1-qu*(s+iσ0)s2+σ20-2u*+m+1q-iσ0qu*+τ0se-iσ0τ0.

    根據(jù)文獻(xiàn)[11]可知, AU是線性系統(tǒng)(16)的無窮小生成元, 滿足AUψ=(θ

    ). 同時(shí), 令f0=(f10,f20)T

    , 其中f10=(1,0)T, f20=(0,1)T, c=(c1,c2)T∈瘙綇

    2, 定義c·f0=c1f10+c2f20, 當(dāng)ψ(θ)∈[-1,0]時(shí), 定義(ψ·f0)(θ)=ψ(θ)·f0,

    且對u=(u1,u2), v=(u1,u2)∈X=C([0,lπ],瘙綆2), 有

    〈u,v〉=1lπ∫lπ0u1v1dx+1lπ

    ∫lπ0u2v2dx.

    因此, 〈φ,f0〉=(〈φ

    ,f10〉,〈φ,f20〉)T, 其中φ∈C([-1,0],X).

    令Φ=(q(θ),q(θ)), Ψ=(q*(s),q

    *(s))T, 則(Ψ,Φ)0=I, 其中I=10

    01. 設(shè)系統(tǒng)(17)的中心子空間為P=span{q(θ),q(θ)},

    伴隨子空間為P*=span{q*(s),q*(s)}.

    下面再由PCNC給出系統(tǒng)(17)的中心子空間, 其中

    PCNφ=Φ(Ψ,〈φ,f0〉)0·f0, φ

    ∈C, PCNC={(q(θ)z,q(θ)z)·f0: z∈瘙綇},

    C=PCNCPSC, PSC為穩(wěn)定的子空間.

    由于Hopf分支的方向和穩(wěn)定性都與μ=0有關(guān), 因此在系統(tǒng)(15)中令μ=0, 可確定中心流形

    W(z,z)=W20(θ)z22+W11(θ)z

    z+W02(θ)z22+…,

    系統(tǒng)(15)在中心流形中的流可寫為

    ut=Φ(z(t),z(t))T+W(z(t),z(t))=q(θ)z(t)+

    q(θ)z(t)+W(z(t),z(t)),

    其中

    (t)=iω0τ0z(t)+q*(0)〈G(ut,0),f0〉=iω0τ

    0z(t)+g(z,z),

    g(z,z)=q*(0)〈G(ut,0),f0〉

    =g20z22+g11zz+g02z22+g21z2z2+….(18)

    再根據(jù)G(φ,μ)的表達(dá)式, 有

    G(φ,0)=τ0F0(φ)=τ0G1G2,

    其中,G1= "(-6u*+2m-2)φ21(0)-6φ31(0)-qφ1(0)φ2(0)+O(4),

    G2=-2sv*2u*3φ21(-1)+6sv*2u*4φ31(-1)+12sv*u*4φ

    21(-1)φ2(0)+2sv*u*2φ1(-1)φ2(0)+2su*2φ1(-1)φ22(0)-2su*φ

    22(0)+O(4),(19)

    式中O(4)=O(‖(u,v)‖4).

    由式(18),(19)可得

    g20= "2q*1τ0D[(-6u*+2m-2)q21-qq1q2]+2q*2τ0D

    -2sv*2u*3q21e-2iσ0τ0+2sv*u*2q1q2e-iσ0τ0-2su*q22,

    g11= "q*1τ0D[2(-6u*+2m-2)q1q1-q(q1q2+q1q2)]

    + "q*2τ0D-4sv * 2 u*3q1q1+2sv*u*2(q1q

    2eiσ0τ0+q1q2e-iσ0τ0)-4su*q2q2,

    g02= "2q*1τ0D[(-6u*+2m-2)q21-qq1q2]+2q*2τ0D

    -2sv*2u*3q21e2iσ0τ0+2sv*u*2q

    1q2eiσ0τ0-2su*q22,

    g21= "2q*1τ0D[(-6u*+2m-2)(2q1W111(0)+q1W120(0))-2q21q21]-

    q*1τ0βD[2qq21W211(0)+qq21W220(0)+qq2W120(0)+2qq2W111(0)]-

    2q*2τ0D2sv*2u*3(2q1e-iσ0τ0W111(-1)+q1

    eiσ0τ0W120(-1))+2q*2τ0D

    6sv*2u*4(2q1q21e2iσ0τ0)+12sv*u*4(2q21q2e-2iσ0τ0)

    + "q*2τ0D2sv*u*2(2q2W111(-1)+2q2W120(-1)+q1e

    iσ0τ0W220(0)+2q1e-iσ0τ0W211(0))- "2q*2τ0D

    2su*(2q2W211(0)+q2W220(0))

    +2q*2τ0D2su*2(q1q22e-iσ0τ0),

    因此, 要得到g21, 需要先計(jì)算出W20(θ)和W11(θ).

    由于W(z(t),z(t))滿足

    = "AUW+X0G

    (Φ(z,z)T·f0+w(z,z),0)

    - "Φ(Ψ,〈X0G(

    Φ(z,z)T·f0+w(z,z),0),f0〉)0·f0

    = "A(0)W+H20 z22+H11zz+H

    02z22+…,(20)

    因此根據(jù)鏈?zhǔn)椒▌t

    =W(z,z)z·+W(z,z)z

    ·z·,

    可得

    [2iσ0τ0-AU]W20=H20,

    -AUW11=H11,

    [-2iσ0τ0-AU]W02=H02.(21)

    當(dāng)-1≤θlt;0時(shí), 有

    -Φ(Ψ,〈X0G

    (Φ(z,z)T·f0+w(z,z),0),f0〉)0·f0=H20z

    22+H11zz+H02z22+…,(22)

    H20(θ)=-[g20q(θ)+g02q(θ)]·f0,

    H11(θ)=-[g11q(θ)+g11q(θ)]·f0.

    根據(jù)式(21),(22)可得

    W20(θ)=ig20σ0τ0q(θ)·f0+ig023σ0

    τ0q(θ)·f0+E1e2iσ0τ0θ,

    W11(θ)=-ig11σ0τ0q(θ)·f0+ig11σ

    0τ0q(θ)·f0+E2.

    在式(21)中令θ=0, 再結(jié)合無窮小生成元的定義以及

    H20(0)=-g20[q(0)+g02q(0)]·f0+τ0

    (-6u*+2m-2)q21-qq1q2-2sv*2u*3q21e-2iσ0τ0+2sv*u*2q1q2e-iσ0τ0-2s

    u*q22,H11(0)=-g11[q(0)+g11 q

    (0)]·f0+τ02(-6u*+2m-2)q1q1-q(q1q2+q1q2)

    -4sv*2u*3q1q1+2sv*u*2

    (q1q2eiσ0τ0+q1q2e-iσ0τ0)-4su*q2q2,

    得E1=E11·E12,

    E2=E21·E22,(23)

    其中

    E11=2iσ0+2u*2-(m+1)u*qu*-se-2iσ0τ02iσ0+s-1,

    E12=(-6u*+2m-2)q21-qq1q2-2sv*2u*3q21e-2iσ0τ0

    +2sv*u*2q1q2e-iσ0τ0-2su*q22,

    E21=2u*2-(m+1)u*qu*-ss-1," E22=

    2(-6u*+2m-2)q1q1-q(q1q2+q1q2)-4sv*2u*3q1q1+2sv*u*

    2(q1q2eiσ0τ0+q1q2e-iσ0 τ0)-4su*q2q2.

    因此, g21表達(dá)式由此可確定.

    基于上述分析, 根據(jù)參數(shù)可計(jì)算出每個(gè)gij的值, 判斷分支方向和分支周期解穩(wěn)定性的表達(dá)式也可由計(jì)算得到:

    C1(0)=i2σ0τ00g11g20-2g112-g0223

    +g212," μ2=-Re(C1(0))Re(λ′(τ00)),β2=2Re(C1(0)),

    T2=-Im(C1(0))+μ2Im(λ′(τ00))σ0τ00.

    因此可得如下定理.

    定理2 對于系統(tǒng)(4), 有以下結(jié)論:

    1) μ2確定Hopf分支的分支方向, 當(dāng)μ2gt;0時(shí), 分支方向是超臨界的, 當(dāng)μ2lt;0時(shí), 分支方向是次臨界的;

    2) β2確定分支周期解的穩(wěn)定性, 當(dāng)β2lt;0時(shí), 分支周期解是漸近穩(wěn)定的, 當(dāng)β2gt;0時(shí), 分支周期解是不穩(wěn)定的;

    3) T2確定分支周期解的周期, 當(dāng)T2gt;0時(shí), 周期增大, 當(dāng)T2lt;0時(shí), 周期減少.

    3 數(shù)值模擬

    下面利用MATLAB工具給出具體的實(shí)例, 以驗(yàn)證上述理論分析結(jié)果.

    對于系統(tǒng)(4), 取參數(shù)m=0.29, q=0.2, s=0.185, l=2.25, d1=0.01, d2=1, 計(jì)算得到正平衡點(diǎn)為E1(u*,v*)=(0.628 8,0.628 8), 同時(shí)滿足條件

    0lt;qlt;(1-m)2, sgt;s*=-2u*2+(m+1)u*,

    且u*∈m+12,1, 其中0lt;mlt;1, 根據(jù)計(jì)算可得τ00=8.027 0.

    取τ=4.54lt;τ00, 由定理1和定理2可知, 當(dāng)τ經(jīng)過τ00, 正平衡點(diǎn)E1(u*,v*)由穩(wěn)定變?yōu)椴环€(wěn)定并產(chǎn)生Hopf分支, 如圖1~圖3所示. 取初值

    u(x,t)=u*+0.01sin(4x), v(x,t)=v*+0.01sin(4x), τ=4.54lt;8.027 0,

    由定理1和定理2可知, 系統(tǒng)(4)在正平衡點(diǎn)E1(u*,v*)處是局部漸近穩(wěn)定的, 如圖1所示.

    取τ=8.03gt;τ00時(shí), 由定理1和定理2可知, 系統(tǒng)(4)在正平衡點(diǎn)E1(u*,v*)處產(chǎn)生Hopf分支, 如圖2和圖3所示. 當(dāng)初值為

    u(x,t)=u*+0.15sin(0.01x), v(x,t)=v*+0.15sin(0.01x)

    時(shí), 系統(tǒng)(4)產(chǎn)生空間齊次周期解, 如圖2所示; 當(dāng)初值為

    u(x,t)=u*+0.014sin(0.5x)," v(x,t)=v*+0.014sin(0.5x)

    時(shí), 系統(tǒng)(4)產(chǎn)生空間非齊次周期解, 如圖3所示.

    圖1 當(dāng)τ=4.54lt;τ00時(shí)系統(tǒng)(4)在正平衡點(diǎn)處穩(wěn)定

    Fig.1 System (4) is stable at" positive equilibrium point when parameter τ=4.54lt;τ00

    圖2 當(dāng)τ=8.03gt;τ00時(shí)系統(tǒng)(4)產(chǎn)生齊次周期解

    Fig.2 System (4) generates homogeneous periodic solutions when parameter τ=8.03gt;τ00

    圖3 當(dāng)τ=8.03gt;τ00時(shí)系統(tǒng)(4)

    產(chǎn)生非齊次周期解Fig.3 System (4) generates inhomogeneous periodic solutions when parameter τ=8.03gt;τ00

    此外, 計(jì)算得到μ2=1.348 1×105gt;0, 產(chǎn)生的Hopf分支方向?yàn)槌R界, β2=-1.385 7×103lt;0, 分支周期解是漸近穩(wěn)定的, T2=-6.594 1×103lt;0, 分支周期解的周期減少.

    參考文獻(xiàn)

    [1] CHEN S S, SHI J P, WEI J J. The Effect of Delay on a Diffusive Preda

    tor-Prey System with Holling Type-Ⅱ Predator Functional Response [J]. Communications on Pure and Applied Analysis, 2013, 12(1): 481-501.

    [2] 王雅迪, 袁海龍. 時(shí)滯Lengyel-Epstein反應(yīng)擴(kuò)散系統(tǒng)的Hopf分支 [J]. 山東大學(xué)學(xué)報(bào)(理學(xué)版), 2023, 58(8): 92-103. (WANG Y D, YUAN H L. H

    opf Bifurcation Analysis in the Lengyel-Epstein Reaction Diffusion System with Time Delay [J]. Journal of Shandong University (Natural Science), 2023, 58(8): 92-103.)

    [3] 姚佳佳, 沈維. 一類食餌-捕食模型的穩(wěn)定性和Hopf分支的存在性 [J]. 吉林大學(xué)學(xué)報(bào)(理學(xué)版), 2022, 60(2): 225-230. (YAO J J, SHEN W. Sta

    bility of a Class of Prey-Predator Model and Existence of Hopf Bifurcation [J]. Journal of Jilin University (Science Edition), 2022, 60(2): 225-230.)

    [4] CHANG X Y, WEI J J. Hopf Bifurcation and Optimal Con

    trol in a Diffusive Predator-Prey System with Time Delay and Prey Harvesting [J]. Nonlinear Analysis: Modelling and Control, 2012, 17(4): 379-409.

    [5] FU S M, ZHANG H S. Effect of Hunting Cooperation on the Dynamic Behav

    ior for a Diffusive Holling Type Ⅱ Predator-Prey Model [J]. Communications in Nonlinear Science and Numerical Simulation, 2021, 99: 105807-1-105807-23.

    [6] NI W J, WANG M X. Dynamical Properties of a Leslie-

    Gower Prey-Predator Model with Strong Allee Effect in Prey [J]. Discrete and Continuous Dynamical Systems: Series B, 2017, 22(9): 3409-3420.

    [7] YANG R Z, ZHANG C R. Dynamics in a Diffusive Predator-Prey System wit

    h a Constant Prey Refuge and Delay [J]. Nonlinear Analysis: Real World Applications, 2016, 31: 1-22.

    [8] CHEN S S, SHI J P, WEI J J. Global Stability and Hopf Bifurcation in

    a Delayed Diffusive Leslie-Gower Predator-Prey System [J]. International Journal of Bifurcation and Chaos, 2012, 22(3): 1250061-1-1250061-11.

    [9] DU Y F, NIU B, WEI J J. Two Delays Induce Hopf Bifur

    cation and Double Hopf Bifurcation in a Diffusive Leslie-Gower Predator-Prey S

    ystem [J]. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2019, 29(1): 013101-1-013101-16.

    [10] 常笑源, 汪紅, 楊鴻飛, 等. 具時(shí)滯和Leslie-Gower功能反應(yīng)函數(shù)的捕食者-食餌系統(tǒng)的Hopf分支分析 [J]. 黑龍江大學(xué)自然科學(xué)學(xué)報(bào), 2020

    , 37(1): 39-44. (CHANG X Y, WANG H, YANG H F, et al. Analysis of Hopf Bifur

    cation on the Leslie-Gower Predator-Prey System with Time Delay [J]. Journal of Natural Science of Heilongjiang University, 2020, 37(1): 39-44.)

    [11] WU J H. Theory and Applications of Partial Function

    al Differential Equations [M]. Applied Mathematical Sciences. New York: Springer, 1996: 1-448.

    [12] LIN X D, SO J W H, WU J H. Center Manifolds for Partial Differential E

    quations with Delays [J]. Proceedings of the Royal Society of Edinburgh Section A: Mathematics, 1992, 122(3/4): 237-254.

    [13] HASSARD B D, KAZARINOFF N D, WAN Y H. Theory and

    Applications of Hopf Bifurcation [M]. Cambridge: Cambridge University Press, 1981: 1-499.

    (責(zé)任編輯: 趙立芹)

    猜你喜歡
    時(shí)滯穩(wěn)定性
    一類k-Hessian方程解的存在性和漸近穩(wěn)定性
    SBR改性瀝青的穩(wěn)定性評價(jià)
    石油瀝青(2021年4期)2021-10-14 08:50:44
    帶有時(shí)滯項(xiàng)的復(fù)Ginzburg-Landau方程的拉回吸引子
    非線性中立型變延遲微分方程的長時(shí)間穩(wěn)定性
    半動(dòng)力系統(tǒng)中閉集的穩(wěn)定性和極限集映射的連續(xù)性
    中立型Emden-Fowler時(shí)滯微分方程的振動(dòng)性
    一階非線性時(shí)滯微分方程正周期解的存在性
    一類時(shí)滯Duffing微分方程同宿解的存在性
    模糊微分方程的一致穩(wěn)定性
    一類離散非線性切換系統(tǒng)的穩(wěn)定性
    av在线天堂中文字幕| 国产精品二区激情视频| 亚洲精品中文字幕一二三四区| 香蕉国产在线看| 欧美国产日韩亚洲一区| 熟妇人妻久久中文字幕3abv| 99国产精品一区二区三区| 在线观看午夜福利视频| 日本 欧美在线| 在线观看日韩欧美| 男人操女人黄网站| 美女扒开内裤让男人捅视频| 99国产精品99久久久久| 曰老女人黄片| 老司机在亚洲福利影院| 日本五十路高清| 久热爱精品视频在线9| 久久香蕉精品热| 午夜两性在线视频| 老熟妇仑乱视频hdxx| 自线自在国产av| 丁香欧美五月| 在线免费观看的www视频| 国内毛片毛片毛片毛片毛片| 老司机午夜福利在线观看视频| 9191精品国产免费久久| 日韩成人在线观看一区二区三区| 18禁美女被吸乳视频| 日韩一卡2卡3卡4卡2021年| 国产精品野战在线观看| 老司机靠b影院| 午夜福利,免费看| 十八禁人妻一区二区| 成人18禁在线播放| 亚洲欧美激情在线| 91精品三级在线观看| 18美女黄网站色大片免费观看| 欧美日韩亚洲国产一区二区在线观看| 国产精品一区二区在线不卡| 岛国在线观看网站| 男人舔女人的私密视频| 国产精品免费一区二区三区在线| 如日韩欧美国产精品一区二区三区| 日本三级黄在线观看| 成人国产综合亚洲| 亚洲av美国av| 亚洲人成电影免费在线| 色播在线永久视频| 日本欧美视频一区| 一边摸一边抽搐一进一小说| 精品一区二区三区av网在线观看| 日本免费一区二区三区高清不卡 | 91国产中文字幕| 色尼玛亚洲综合影院| 日韩精品中文字幕看吧| 久久精品国产清高在天天线| 男人操女人黄网站| 黑人操中国人逼视频| 伊人久久大香线蕉亚洲五| 美女高潮到喷水免费观看| 日韩有码中文字幕| 黄色毛片三级朝国网站| 午夜福利高清视频| 叶爱在线成人免费视频播放| 久久这里只有精品19| 色av中文字幕| 久久人妻av系列| 成人精品一区二区免费| 亚洲欧美激情在线| 无遮挡黄片免费观看| 国产亚洲欧美在线一区二区| 亚洲国产欧美一区二区综合| 国产欧美日韩一区二区精品| 国产国语露脸激情在线看| 亚洲欧美精品综合久久99| 两个人看的免费小视频| 韩国av一区二区三区四区| 国语自产精品视频在线第100页| 每晚都被弄得嗷嗷叫到高潮| 国产一区在线观看成人免费| 欧美黄色片欧美黄色片| 国产一区二区三区在线臀色熟女| 一个人观看的视频www高清免费观看 | 老司机靠b影院| 午夜a级毛片| 亚洲五月天丁香| 国产人伦9x9x在线观看| 亚洲精品中文字幕在线视频| 日韩视频一区二区在线观看| 亚洲av电影不卡..在线观看| 久久精品91蜜桃| 亚洲成av片中文字幕在线观看| 久久久国产成人免费| 一级毛片精品| 欧美av亚洲av综合av国产av| 51午夜福利影视在线观看| 中出人妻视频一区二区| 午夜激情av网站| 日韩欧美一区二区三区在线观看| tocl精华| 一级毛片高清免费大全| 欧美av亚洲av综合av国产av| 午夜免费激情av| 久久热在线av| 亚洲av成人不卡在线观看播放网| 久久久久久久精品吃奶| 久久热在线av| 男女午夜视频在线观看| 免费看十八禁软件| 久久久久久久久中文| 免费看a级黄色片| 日日干狠狠操夜夜爽| 午夜成年电影在线免费观看| 国产99白浆流出| 久久香蕉国产精品| 十分钟在线观看高清视频www| www日本在线高清视频| 精品一品国产午夜福利视频| 别揉我奶头~嗯~啊~动态视频| 黄色毛片三级朝国网站| 欧美av亚洲av综合av国产av| 桃红色精品国产亚洲av| 搡老熟女国产l中国老女人| 波多野结衣av一区二区av| 制服丝袜大香蕉在线| 亚洲国产中文字幕在线视频| 色哟哟哟哟哟哟| 亚洲人成77777在线视频| 亚洲免费av在线视频| 男男h啪啪无遮挡| 国产精品久久久av美女十八| 在线观看免费日韩欧美大片| 精品久久久精品久久久| 精品国产美女av久久久久小说| 给我免费播放毛片高清在线观看| 在线国产一区二区在线| av视频免费观看在线观看| 午夜福利视频1000在线观看 | 丝袜美腿诱惑在线| 亚洲欧洲精品一区二区精品久久久| 亚洲精品国产区一区二| 大陆偷拍与自拍| 中文字幕色久视频| 国产99白浆流出| 亚洲av电影在线进入| 精品久久蜜臀av无| 女同久久另类99精品国产91| 亚洲人成77777在线视频| 成人三级做爰电影| 少妇粗大呻吟视频| 男女下面插进去视频免费观看| 成人18禁高潮啪啪吃奶动态图| 成人av一区二区三区在线看| 最近最新中文字幕大全免费视频| 窝窝影院91人妻| 午夜影院日韩av| 夜夜爽天天搞| 亚洲av片天天在线观看| 久久精品国产清高在天天线| 女人精品久久久久毛片| 欧美日韩精品网址| 色综合站精品国产| 丰满人妻熟妇乱又伦精品不卡| 亚洲精品国产区一区二| cao死你这个sao货| 国产精品野战在线观看| 国产又色又爽无遮挡免费看| 欧美日韩一级在线毛片| 欧美乱色亚洲激情| 亚洲精品美女久久av网站| 欧美丝袜亚洲另类 | 亚洲av五月六月丁香网| 在线观看www视频免费| 午夜福利一区二区在线看| а√天堂www在线а√下载| 久久这里只有精品19| 国产在线精品亚洲第一网站| 国产精品98久久久久久宅男小说| 久久久国产欧美日韩av| 久久久久久国产a免费观看| 久久午夜亚洲精品久久| 亚洲国产中文字幕在线视频| 两性午夜刺激爽爽歪歪视频在线观看 | 99国产精品一区二区蜜桃av| 国产区一区二久久| 91大片在线观看| 欧美国产精品va在线观看不卡| 久久国产精品男人的天堂亚洲| 亚洲av第一区精品v没综合| 黄色丝袜av网址大全| 精品电影一区二区在线| 视频在线观看一区二区三区| 亚洲av成人一区二区三| 人人妻人人爽人人添夜夜欢视频| 99久久精品国产亚洲精品| 少妇 在线观看| 亚洲中文av在线| 免费久久久久久久精品成人欧美视频| 亚洲伊人色综图| 夜夜爽天天搞| 亚洲专区字幕在线| 国产精品一区二区三区四区久久 | 国产高清videossex| 老熟妇仑乱视频hdxx| 亚洲av熟女| 看免费av毛片| 国产精品香港三级国产av潘金莲| 免费在线观看亚洲国产| 国产亚洲av高清不卡| 色综合站精品国产| 精品一区二区三区av网在线观看| 亚洲精品久久成人aⅴ小说| 99香蕉大伊视频| 国产亚洲精品av在线| 亚洲视频免费观看视频| 色综合欧美亚洲国产小说| 国产99白浆流出| а√天堂www在线а√下载| 看免费av毛片| 在线视频色国产色| 在线观看免费午夜福利视频| 亚洲欧美日韩另类电影网站| 久久人妻熟女aⅴ| 怎么达到女性高潮| 国产一卡二卡三卡精品| 免费久久久久久久精品成人欧美视频| 91成人精品电影| 99riav亚洲国产免费| 99久久99久久久精品蜜桃| 亚洲激情在线av| 女人精品久久久久毛片| 身体一侧抽搐| 亚洲av熟女| 看片在线看免费视频| 国产高清激情床上av| bbb黄色大片| 亚洲欧美精品综合一区二区三区| 日本免费a在线| 亚洲在线自拍视频| 日韩欧美一区二区三区在线观看| 亚洲欧美精品综合久久99| ponron亚洲| 麻豆av在线久日| 午夜福利影视在线免费观看| 成人免费观看视频高清| 国产99久久九九免费精品| av在线播放免费不卡| 日本一区二区免费在线视频| 亚洲中文字幕一区二区三区有码在线看 | 免费在线观看日本一区| 丁香欧美五月| 可以在线观看的亚洲视频| 国产精品一区二区精品视频观看| 搡老妇女老女人老熟妇| 香蕉久久夜色| 91精品三级在线观看| 欧洲精品卡2卡3卡4卡5卡区| 欧美在线黄色| 国产精品秋霞免费鲁丝片| 国产亚洲精品久久久久久毛片| 丝袜美腿诱惑在线| 99久久99久久久精品蜜桃| 岛国视频午夜一区免费看| 亚洲欧美精品综合久久99| 成人国产一区最新在线观看| 麻豆国产av国片精品| 777久久人妻少妇嫩草av网站| 国产精品爽爽va在线观看网站 | 亚洲一区二区三区不卡视频| 老汉色av国产亚洲站长工具| 人人澡人人妻人| 中文字幕人成人乱码亚洲影| 最新在线观看一区二区三区| 91字幕亚洲| 嫁个100分男人电影在线观看| 免费在线观看视频国产中文字幕亚洲| 欧美日韩中文字幕国产精品一区二区三区 | 国产乱人伦免费视频| 一个人免费在线观看的高清视频| 日日爽夜夜爽网站| 久久久水蜜桃国产精品网| svipshipincom国产片| 久久人人97超碰香蕉20202| 黄片大片在线免费观看| 一二三四在线观看免费中文在| 叶爱在线成人免费视频播放| 国产亚洲精品av在线| 天堂√8在线中文| 国产精品久久电影中文字幕| 日韩精品免费视频一区二区三区| 国产单亲对白刺激| 村上凉子中文字幕在线| 老司机午夜十八禁免费视频| 亚洲 国产 在线| 禁无遮挡网站| 亚洲色图av天堂| 亚洲九九香蕉| 亚洲色图综合在线观看| 十八禁网站免费在线| 久久久久精品国产欧美久久久| 丁香欧美五月| 99久久精品国产亚洲精品| 国产精品综合久久久久久久免费 | 色婷婷久久久亚洲欧美| 久久人妻福利社区极品人妻图片| 香蕉久久夜色| 视频区欧美日本亚洲| 亚洲第一电影网av| 999久久久国产精品视频| 久久亚洲真实| 丁香六月欧美| 涩涩av久久男人的天堂| 亚洲av熟女| 精品国内亚洲2022精品成人| 日本五十路高清| 日韩高清综合在线| 国产亚洲精品av在线| 国产精品永久免费网站| 国产亚洲av高清不卡| av有码第一页| 久久久久久亚洲精品国产蜜桃av| 国产精品99久久99久久久不卡| 中文字幕人妻丝袜一区二区| 美女 人体艺术 gogo| 少妇熟女aⅴ在线视频| 欧美老熟妇乱子伦牲交| 两人在一起打扑克的视频| 午夜激情av网站| 亚洲中文av在线| 欧美丝袜亚洲另类 | 国产精品1区2区在线观看.| 国产91精品成人一区二区三区| 美女扒开内裤让男人捅视频| av欧美777| 国产单亲对白刺激| 久久国产精品男人的天堂亚洲| 国产97色在线日韩免费| 97碰自拍视频| 18禁黄网站禁片午夜丰满| 一区二区三区高清视频在线| 日韩有码中文字幕| 韩国精品一区二区三区| 老司机午夜十八禁免费视频| 久久人人精品亚洲av| 精品国产一区二区久久| 此物有八面人人有两片| 日本vs欧美在线观看视频| 久久欧美精品欧美久久欧美| 久久这里只有精品19| 欧美+亚洲+日韩+国产| 99精品久久久久人妻精品| 黄色 视频免费看| 亚洲五月婷婷丁香| 久久久久久免费高清国产稀缺| 中国美女看黄片| 在线观看66精品国产| 啪啪无遮挡十八禁网站| 亚洲欧美日韩另类电影网站| 黄片大片在线免费观看| 久久婷婷成人综合色麻豆| 国产亚洲av嫩草精品影院| 亚洲 欧美 日韩 在线 免费| 如日韩欧美国产精品一区二区三区| 亚洲电影在线观看av| 久久性视频一级片| 欧美亚洲日本最大视频资源| 可以在线观看的亚洲视频| 999久久久国产精品视频| 免费在线观看黄色视频的| 午夜精品国产一区二区电影| 一进一出抽搐动态| 涩涩av久久男人的天堂| 亚洲欧美精品综合久久99| 国产xxxxx性猛交| 午夜福利成人在线免费观看| 少妇被粗大的猛进出69影院| 国产精品一区二区免费欧美| 亚洲熟妇中文字幕五十中出| 99精品在免费线老司机午夜| or卡值多少钱| 国产亚洲精品综合一区在线观看 | 国产欧美日韩精品亚洲av| 怎么达到女性高潮| 午夜福利免费观看在线| 一本综合久久免费| 国产高清有码在线观看视频 | svipshipincom国产片| 国产av一区二区精品久久| 多毛熟女@视频| 他把我摸到了高潮在线观看| xxx96com| 一二三四社区在线视频社区8| 欧美一级a爱片免费观看看 | 人人妻人人爽人人添夜夜欢视频| 亚洲精品av麻豆狂野| 亚洲精品粉嫩美女一区| 亚洲国产欧美网| 黑丝袜美女国产一区| 美国免费a级毛片| 午夜亚洲福利在线播放| 97人妻天天添夜夜摸| 18美女黄网站色大片免费观看| 国产亚洲av高清不卡| 青草久久国产| 999精品在线视频| 丝袜人妻中文字幕| 超碰成人久久| 淫秽高清视频在线观看| 色综合欧美亚洲国产小说| 成人特级黄色片久久久久久久| 51午夜福利影视在线观看| 一级毛片女人18水好多| 国产又色又爽无遮挡免费看| 国产国语露脸激情在线看| 一级黄色大片毛片| av在线天堂中文字幕| 久久午夜综合久久蜜桃| 午夜福利欧美成人| 中文字幕人妻熟女乱码| 黄网站色视频无遮挡免费观看| 久久精品国产亚洲av高清一级| 精品高清国产在线一区| 亚洲,欧美精品.| 自拍欧美九色日韩亚洲蝌蚪91| 中文亚洲av片在线观看爽| 精品免费久久久久久久清纯| 中文字幕高清在线视频| 深夜精品福利| 国产成人啪精品午夜网站| 亚洲成人精品中文字幕电影| 国产欧美日韩综合在线一区二区| 国产精品免费视频内射| 久久人妻熟女aⅴ| 男男h啪啪无遮挡| 91精品三级在线观看| 国产av一区在线观看免费| 午夜福利一区二区在线看| 国产精品一区二区免费欧美| 久久香蕉国产精品| 超碰成人久久| 在线观看午夜福利视频| 999精品在线视频| 精品国产一区二区久久| 日韩 欧美 亚洲 中文字幕| 日韩欧美国产在线观看| 国产一区二区三区综合在线观看| 色播在线永久视频| 日韩av在线大香蕉| 国产午夜福利久久久久久| 精品熟女少妇八av免费久了| 在线观看舔阴道视频| 久久久国产成人免费| 夜夜躁狠狠躁天天躁| 亚洲片人在线观看| 自拍欧美九色日韩亚洲蝌蚪91| 无遮挡黄片免费观看| 欧美黄色淫秽网站| 无遮挡黄片免费观看| 97超级碰碰碰精品色视频在线观看| 午夜福利高清视频| 无人区码免费观看不卡| 欧美激情久久久久久爽电影 | 亚洲va日本ⅴa欧美va伊人久久| 国产伦一二天堂av在线观看| 亚洲av美国av| www.熟女人妻精品国产| 欧美在线黄色| 欧美一级毛片孕妇| 男女午夜视频在线观看| 精品一区二区三区四区五区乱码| 精品国产超薄肉色丝袜足j| 日本 欧美在线| 久久精品国产综合久久久| 欧美黑人欧美精品刺激| 国产精品 欧美亚洲| 久久久久亚洲av毛片大全| 麻豆国产av国片精品| 国产一区二区三区视频了| 一区二区三区精品91| 午夜两性在线视频| 日韩高清综合在线| 变态另类成人亚洲欧美熟女 | 波多野结衣巨乳人妻| 女生性感内裤真人,穿戴方法视频| 啦啦啦观看免费观看视频高清 | 啦啦啦韩国在线观看视频| 国产欧美日韩精品亚洲av| 啦啦啦观看免费观看视频高清 | 热re99久久国产66热| 成人免费观看视频高清| 免费一级毛片在线播放高清视频 | 久久 成人 亚洲| 欧美黑人欧美精品刺激| 制服人妻中文乱码| 一进一出好大好爽视频| 首页视频小说图片口味搜索| 一边摸一边抽搐一进一小说| 亚洲欧美激情在线| 午夜精品在线福利| 不卡一级毛片| 国产av又大| 午夜福利高清视频| 97超级碰碰碰精品色视频在线观看| 村上凉子中文字幕在线| 99久久综合精品五月天人人| 精品人妻在线不人妻| 国产高清videossex| 黄色女人牲交| 久久久久精品国产欧美久久久| 久久国产乱子伦精品免费另类| 午夜福利成人在线免费观看| 国产精品久久久人人做人人爽| 国产精品 欧美亚洲| 99热只有精品国产| www.精华液| 婷婷六月久久综合丁香| 天堂影院成人在线观看| 久久亚洲精品不卡| 成人18禁高潮啪啪吃奶动态图| 亚洲av成人av| 日日摸夜夜添夜夜添小说| 俄罗斯特黄特色一大片| 色婷婷久久久亚洲欧美| 一个人观看的视频www高清免费观看 | 美女 人体艺术 gogo| 国产精品影院久久| 亚洲第一欧美日韩一区二区三区| 久久婷婷人人爽人人干人人爱 | 91精品三级在线观看| 丝袜美足系列| 欧美日韩亚洲综合一区二区三区_| 亚洲自偷自拍图片 自拍| 久久久久久久久久久久大奶| 一边摸一边做爽爽视频免费| 99在线视频只有这里精品首页| 成人国产一区最新在线观看| 日韩一卡2卡3卡4卡2021年| 午夜福利18| 少妇的丰满在线观看| 日韩 欧美 亚洲 中文字幕| 国产片内射在线| 国产精品,欧美在线| 国产精品一区二区免费欧美| 亚洲久久久国产精品| 精品久久久久久久人妻蜜臀av | 久久久久久亚洲精品国产蜜桃av| 村上凉子中文字幕在线| 如日韩欧美国产精品一区二区三区| 亚洲国产欧美日韩在线播放| bbb黄色大片| 最好的美女福利视频网| 可以免费在线观看a视频的电影网站| 国产真人三级小视频在线观看| 欧洲精品卡2卡3卡4卡5卡区| 宅男免费午夜| 亚洲欧美日韩高清在线视频| av有码第一页| 美国免费a级毛片| 最近最新中文字幕大全免费视频| 国产精品98久久久久久宅男小说| www.www免费av| 国产高清视频在线播放一区| 操美女的视频在线观看| 亚洲第一av免费看| 多毛熟女@视频| 日韩欧美三级三区| 黄色毛片三级朝国网站| 日本免费a在线| 久久国产精品影院| 久久久国产成人免费| 国产亚洲精品第一综合不卡| 日韩欧美在线二视频| avwww免费| 高潮久久久久久久久久久不卡| 亚洲无线在线观看| 中文字幕精品免费在线观看视频| 丝袜人妻中文字幕| 97碰自拍视频| 国产亚洲欧美精品永久| 在线免费观看的www视频| 精品国产国语对白av| av超薄肉色丝袜交足视频| 亚洲精品国产区一区二| 国产单亲对白刺激| 18禁国产床啪视频网站| 色综合欧美亚洲国产小说| 亚洲激情在线av| 国产激情欧美一区二区| 色综合站精品国产| 亚洲最大成人中文| 18美女黄网站色大片免费观看| 免费在线观看影片大全网站| 国产麻豆69| 日韩精品免费视频一区二区三区| 黄片播放在线免费| 亚洲性夜色夜夜综合| 欧美激情久久久久久爽电影 | 9热在线视频观看99| 十八禁人妻一区二区| 午夜视频精品福利| 日本vs欧美在线观看视频| 在线观看免费视频网站a站| 色尼玛亚洲综合影院| 国产一区二区三区综合在线观看| 女人爽到高潮嗷嗷叫在线视频| cao死你这个sao货| 亚洲精品国产区一区二| 欧美亚洲日本最大视频资源| 国产麻豆成人av免费视频| 丁香欧美五月| 最新美女视频免费是黄的| 色综合站精品国产| 女性被躁到高潮视频| 国产色视频综合| 亚洲精品一区av在线观看|