Loading [MathJax]/jax/output/SVG/jax.js
Research article

Dynamic properties and numerical simulations of a fractional phytoplankton-zooplankton ecological model

  • Received: 17 January 2025 Revised: 19 May 2025 Accepted: 28 May 2025 Published: 10 June 2025
  • It is always a difficult and hot topic to effectively explore the dynamic behavior of phytoplankton-zooplankton models. Thus, we investigated dynamic properties and a numerical solution of a fractional-order phytoplankton-zooplankton ecological model (PZEM) that incorporates the effects of toxic substances and additional food transmission in the environment. First, stability, Turing instability, Hopf bifurcation, and weakly nonlinear analysis were analyzed for the PZEM. Second, a new high-precision numerical method was developed for the fractional PZEM without diffusion terms. We compared the method with other methods to determine the effectiveness of the present method. A discretization method was established for the PZEM with diffusion term. Finally, numerical simulation verified the feasibility of the theory. Numerical simulations showed the chaotic attractor and some novel pattern dynamical behaviors of the PZEM.

    Citation: Shuai Zhang, Haolu Zhang, Yulan Wang, Zhiyuan Li. Dynamic properties and numerical simulations of a fractional phytoplankton-zooplankton ecological model[J]. Networks and Heterogeneous Media, 2025, 20(2): 648-669. doi: 10.3934/nhm.2025028

    Related Papers:

    [1] Yimamu Maimaiti, Zunyou Lv, Ahmadjan Muhammadhaji, Wang Zhang . Analyzing vegetation pattern formation through a time-ordered fractional vegetation-sand model: A spatiotemporal dynamic approach. Networks and Heterogeneous Media, 2024, 19(3): 1286-1308. doi: 10.3934/nhm.2024055
    [2] Jiaxin Zhang, Wei Zhang, Xiaoyu Li . Numerical simulation of chaotic dynamics in a fractional-order vibration model with Grünwald-Letnikov fractional derivative. Networks and Heterogeneous Media, 2025, 20(2): 625-647. doi: 10.3934/nhm.2025027
    [3] Toshiyuki Ogawa, Takashi Okuda . Oscillatory dynamics in a reaction-diffusion system in the presence of 0:1:2 resonance. Networks and Heterogeneous Media, 2012, 7(4): 893-926. doi: 10.3934/nhm.2012.7.893
    [4] Yue Qiu, Sara Grundel, Martin Stoll, Peter Benner . Efficient numerical methods for gas network modeling and simulation. Networks and Heterogeneous Media, 2020, 15(4): 653-679. doi: 10.3934/nhm.2020018
    [5] Antonio DeSimone, Martin Kružík . Domain patterns and hysteresis in phase-transforming solids: Analysis and numerical simulations of a sharp interface dissipative model via phase-field approximation. Networks and Heterogeneous Media, 2013, 8(2): 481-499. doi: 10.3934/nhm.2013.8.481
    [6] Yicheng Liu, Yipeng Chen, Jun Wu, Xiao Wang . Periodic consensus in network systems with general distributed processing delays. Networks and Heterogeneous Media, 2021, 16(1): 139-153. doi: 10.3934/nhm.2021002
    [7] Khalid Boushaba . A multi layer method applied to a model of phytoplankton. Networks and Heterogeneous Media, 2007, 2(1): 37-54. doi: 10.3934/nhm.2007.2.37
    [8] Gary Bunting, Yihong Du, Krzysztof Krakowski . Spreading speed revisited: Analysis of a free boundary model. Networks and Heterogeneous Media, 2012, 7(4): 583-603. doi: 10.3934/nhm.2012.7.583
    [9] Riccardo Bonetto, Hildeberto Jardón Kojakhmetov . Nonlinear diffusion on networks: Perturbations and consensus dynamics. Networks and Heterogeneous Media, 2024, 19(3): 1344-1380. doi: 10.3934/nhm.2024058
    [10] Simone Göttlich, Elisa Iacomini, Thomas Jung . Properties of the LWR model with time delay. Networks and Heterogeneous Media, 2021, 16(1): 31-47. doi: 10.3934/nhm.2020032
  • It is always a difficult and hot topic to effectively explore the dynamic behavior of phytoplankton-zooplankton models. Thus, we investigated dynamic properties and a numerical solution of a fractional-order phytoplankton-zooplankton ecological model (PZEM) that incorporates the effects of toxic substances and additional food transmission in the environment. First, stability, Turing instability, Hopf bifurcation, and weakly nonlinear analysis were analyzed for the PZEM. Second, a new high-precision numerical method was developed for the fractional PZEM without diffusion terms. We compared the method with other methods to determine the effectiveness of the present method. A discretization method was established for the PZEM with diffusion term. Finally, numerical simulation verified the feasibility of the theory. Numerical simulations showed the chaotic attractor and some novel pattern dynamical behaviors of the PZEM.



    Due to the difficulty in measuring the biomass of plankton, mathematical models of plankton populations are important methods for understanding the physical and biological processes of plankton. Various phytoplankton-zooplankton models (PZMs) have been studied: An NPZD model [1], a stochastic PZM [2], a plankton-nutrient system [3], a marine phytoplankton-zooplankton ecological model (PZEM) [4], Trophic phytoplankton-zooplankton models (PZM) [5], a bioeconomic PZM with time delays [6], a two-zooplankton one-phytoplankton PZM [7], a phytoplankton and two zooplankton PZM with toxin-producing delay [8], planktonic animal and plant PZM [9], and other PZMs [10, 11, 12]. In [12], Wang et al. investigated the following PZEM:

    {dPdt=plant growthrP(1PK)plant losse1PZ1+e1h1P+e2h2A,dZdt=n1e1PZ+n2e2AZ1+e1h1P+e2h2Aanimal gainτZnatural deathσPZh+Panimal death. (1.1)

    Fractional-order reaction-diffusion equation can describe the dynamic behavior of real systems more accurately. In recent years, fractional-order models have made remarkable progress in fields such as biomedicine, ecology, and public health. In biomedicine, fractional-order models have been used to study tumor growth [13] and hepatitis B virus (HBV) treatment [14], revealing their advantages in describing the complex dynamic behaviors of biological systems. In ecology, fractional-order models have demonstrated their potential in the study of ecosystem complexity by analyzing the dynamic behaviors of plant-herbivore interactions [15]. In public health, the fractional model has been applied to analyze the transmission of smoking behavior [16] and the correlation between human papilloma virus (HPV) and cervical cancer [17], highlighting their utility in public health research. These studies not only enrich the theoretical application of fractional models but also have made significant progress in numerical analysis, providing new tools and methods for solving complex biomedical and ecological problems [13, 17]. For example, in ecology, fractional PZEMs can more accurately describe the complex dynamics of ecosystems, including population fluctuations and interactions. Li et al. [18] studied an established fractional-order delayed zooplankton-phytoplankton model. Javidi and Ahmad [19] studied dynamic analysis of time fractional-order phytoplankton-toxic PZM. Kumar et al. [20] studied a fractional plankton-oxygen modeling. The proposed fractional-order PZEM extends traditional integer-order models by incorporating memory effects through fractional derivatives, which better capture long-term ecological interactions. This model accounts for toxic substances and additional food transmission, factors critical to understanding algal blooms, and population collapse. Models often neglect cross-diffusion and fractional dynamics, limiting their ability to simulate complex spatiotemporal patterns. Our work bridges this gap, offering insights into chaotic attractors and pattern formation under realistic ecological constraints. In this paper, we investigate the following fractional-order PZEM:

    {α1Ptα1=plant growthrP(1PK)plant losse1PZ1+e1h1P+e2h2A+plant diffusiond112P,α2Ztα2=n1e1PZ+n2e2AZ1+e1h1P+e2h2Aanimal gainτZnatural deathσPZh+Panimal death+d212P+d312Zplant and animal diffusion, (1.2)

    where P and Z represent the population densities of phytoplankton and zooplankton, respectively, at time t. 2 stands for the Laplacian operator of the two-dimensional plane, Ω is a bounded connected region with smooth boundaries Ω, and d11 and d31 represent the diffusion coefficients of phytoplankton and zooplankton, respectively. These parameters characterize the stochastic movement of individuals within a population from areas of high concentration to those of lower density. d21 is a cross-diffusion coefficient that indicates the effect of the presence of phytoplankton on the movement of zooplankton populations. αi(i=1,2) is the fractional-order orders, 0<αi<1 and α1Ptα1,α2Ztα2 are the Grünwald-Letnikov (GL) fractional derivative. t is the time variable. The meanings of the parameters in the model are shown in Table 1.

    Table 1.  Parameter and biological significance in the model.
    Parameter Ecological Significance
    r The intrinsic growth rate of phytoplankton
    K The carrying capacity of the environment for phytoplankton
    h1 The time zooplankton need to process one unit of phytoplankton
    h2 The time zooplankton need to process one unit of supplementary food
    e1 The capacity of zooplankton to prey on phytoplankton
    e2 The capacity of zooplankton to forage for supplementary food
    n1 The nutritional value of phytoplankton
    n2 The nutritional value of supplementary food
    A The biomass of supplementary food
    τ The mortality rate of zooplankton
    σ The proportion of toxins released by phytoplankton
    h The half-saturation constant

     | Show Table
    DownLoad: CSV

    Several numerical schemes have been used to solve the Fractional-order reaction-diffusion equation, the Galerkin method [21, 22], the Fourier spectral method [23, 24], the finite difference method [25, 26], the operational matrix method [27, 28], and so on [29, 30, 31, 32, 33, 34]. In this paper, we investigate dynamic properties and numerical solutions of the fractional-order PZEM (1.2). The major contributions of this paper are as follows:

    (a) The concept of the phytoplankton-zooplankton ecological model is first extended to fractional-order PZEM that incorporates the effects of toxic substances and additional food transmission in the environment. This model leverages the memory effect of fractional-order derivatives to more effectively capture biological significance compared to traditional integer-order models.

    (b) Stability, Turing instability, Hopf bifurcation, and weakly nonlinear property are analyzed for the PZEM. a new high-precision numerical method is developed for the fractional PZEM without diffusion term, and a discretization method is established for the PZEM with a diffusion term.

    (c) Numerical simulations show some novel chaotic attractor and pattern dynamical behaviors of the PZEM.

    The paper is organized as follows: In Section 2, we give the stability and Hopf bifurcation analysis. In Section 3, we detail weakly nonlinear analysis. In Section 4, we give the numerical algorithm and the numerical simulation results of the system. Finally, in Section 5, we the summarize paper.

    In this section, we investigate the stability of the equilibrium point and the emergence of Hopf bifurcation. Let

    ˉu=e1h1P,ˉv=e1rZ,t=rT,κ=e1h1K,α=n1h2n2h1,τ=n1rh1,ζ=n2e2h1n1A,m=τr,μ=τr,G=e1hh1.

    Equation (1.2) is converted to the following form:

    {α1tα1ˉu=ˉu(1ˉuκ)ˉuˉv1+αζ+ˉu+d12ˉu,α2tα2ˉv=τ(ˉu+ζ)ˉv1+αζ+ˉumˉvμˉuˉvG+ˉu+d22ˉu+d32ˉv, (2.1)

    where, d1=e1h1d11,d2=e1h1d21,d3=e1rd31.

    We analyze the stability of Eq (2.1) without the diffusion term.

    {α1tα1ˉu=ˉu(1ˉuκ)ˉuˉv1+αζ+ˉu,α2tα2ˉv=τ(ˉu+ζ)ˉv1+αζ+ˉumˉvμˉuˉvG+ˉu. (2.2)

    Solving the following equation:

    {ˉu(1ˉuκ)ˉuˉv1+αζ+ˉu=0,τ(ˉu+ζ)ˉv1+αζ+ˉumˉvμˉuˉvG+ˉu=0. (2.3)

    We can get the following equilibrium points of Eq (2.2):

    E0=(0,0),E1=(κ,0),E=(ˉu,ˉv), (2.4)

    where

    ˉv=(1+αζ+ˉu)(1ˉuκ).

    u represents the positive solution to the quadratic equation (2.5)

    ˉD1ˉu2+ˉD2ˉu+ˉD3=0, (2.5)

    where

    ˉD1=τmμ,ˉD3=τζGmGmαζG,ˉD2=τζ+τGmαζmGμαζμm.

    Solving Eq (2.5), we can get:

    {ˉu1,2=(τζ+τHmαζmHμαζμm)2(τmμ)±(τζ+τHmαζmHμαζμm)24(τmμ)(τζHmHmαζH)2(τmμ),ˉv=(1+αζ+ˉu)(1ˉuκ). (2.6)

    The Jacobian matrix associated with E=(ˉu,ˉv) is formulated as follows

    J=[a11a12a21a22], (2.7)

    where

    a11=12ˉuκˉv(1+αζ)(1+αζ+ˉu)2, a12=ˉu1+αζ+ˉu,a21=τˉv(1+αζζ)(1+αζ+ˉu)2μˉvG(G+ˉu)2, a22=τ(ˉu+ζ)1+αζ+ˉumμˉuG+ˉu.

    When (α1,α2)=(1,1), the characteristic equation at the equilibrium point is as follows:

    λ2+Tr0λ+0=0, (2.8)

    where,

    Tr0=a11+a22,0=a11a22a21a12.

    Let α=0.55,τ=4.1,ζ=0.72,G=8.43,μ=0.79,κ=2.38, and m=2.45. The system's equilibrium point and its associated eigenvalues are presented in Table 2.

    Table 2.  The calculation results of the equilibrium point, eigenvalue, and argument of the system.
    Equilibrium point Eigenvalues of Jacobian matrix Stability
    0 0 –0.3354 1 Unstability
    2.3800 0 –1 0.7421 Unstability
    0.3131 1.4843 0.0138–0.4838i 0.0138+0.4838i Unstability

     | Show Table
    DownLoad: CSV

    It can be seen from Table 3 that all three equilibrium points are unstable points, and the system may generate chaos. Next, we introduce a high-precision numerical method for numerical simulation.

    Table 3.  Comparison of methods for problem (2.9) when β=0.9 and t[0,1].
    Method α=0.4 α=0.7 α=1
    h=0.01 h=0.05 h=0.01 h=0.05 h=0.01 h=0.05
    Closed-from solution [35] 3.2307e-02 5.4302e-02 4.5734e-02 2.5713e-01 2.0716e-01 4.5136e-01
    Present method 8.9463e-05 5.6724e-04 9.1684e-05 7.3524e-04 2.1527e-04 4.5136e-03

     | Show Table
    DownLoad: CSV

    In simulations of this section, let α=1.1,τ=4.4,ζ=0.8,G=8,n=1.75,κ=5.3, and m=2.5 and initial conditions u0=1 and v0=1. We set the time step size h=0.01 and the time as T=600 for simulation. Figure 1 shows numerical results at (α1,α2)=(0.895,0.994).

    Figure 1.  Comparative numerical result of the two methods of the model (α,τ,ζ,G,n,κ,m)=(1.1,4.4,0.8,8,1.75,5.3,2.5), α1=0.895, and α2=0.994.

    We consider the following fractional order GL initial value problem

    {Dαty(t)=tβ,y(0)=0. (2.9)

    Numerical simulations are performed for Eq (2.9). Compared with the closed-from solution, the simulation results of the high-precision numerical method are in good agreement with the simulation results of other methods, and the accuracy is higher. This verifies the effectiveness of the high-precision numerical method.

    We first consider the following GL differential equation:

    dαdtαy(t)=f(t,y(t)). (2.10)

    According to [36], we can get the high-precision numerical method [36, 37, 38, 39]:

    y(tk)=hαf(tk,y(tk1))mi=1ϑ(α,p)jy(tkj), (2.11)

    where,

    {ϑ(α,p)0=g0,k=0,ϑ(α,p)k=1g0k1i=0(1i1+αk)giϑ(α,p)ki,k=1,2,....,p1,ϑ(α,p)k=1g0pi=0(1i1+αk)giϑ(α,p)ki,k=p,p+1,p+2,.... (2.12)
    (g0g1g2gp)=(1111123p+112232(p+1)212p3p(p+1)p)1(012p). (2.13)

    Where m=[tk/h]+1, h is the step size. Therefore, a high-precision numerical method for the system (2.2) is the following

    {ˉuk=hα1(ˉuk1(1ˉuk1κ)ˉuk1ˉvk11+αζ+ˉuk1)+ˉu0mj=1ϑ(α1,p)jˉukj,ˉvk=hα2(τ(ˉuk+ζ)ˉvk11+αζ+ˉukmˉvk1μˉukˉvk1G+ˉuk)+ˉv0mj=1ϑ(α2,p)jˉvkj. (2.14)

    According to [36, 40], the least common order of the system is identified as 0.9819. Let (α1,α2)=(1.05,0.95),α=0.55,τ=4.1,ζ=0.72,G=8.43,μ=0.79,κ=2.38, and m=2.45, and the characteristic equation of the system is given by

    λ401066728823λ212478990071λ19+10396924437077=0, (2.15)

    with characteristic roots λ1,2=1.0358±0.0439i, as |arg(λ1,2)|=0.0423<π40. This indicates that the instability of the system may generate chaos. It can be seen from Figure 2 that the system generates chaos. The chaotic phenomenon of the system gradually disappears over time and tends to a stable state.

    Figure 2.  Numerical simulation of phase diagrams and time series plot of system (2.2) at (α1,α2)=(1.05,0.95)..

    We set different orders (α1,α2)=(1.1,0.95) to observe the influence of the order on the dynamic behavior of the system at parameters α=0.55,τ=4.1,ζ=0.72,G=8.43,n=0.79,κ=2.38, and m=2.45. The characteristic equation of the system is given by

    λ411066728823λ222478990071λ19+10396924437077=0, (2.16)

    with characteristic roots λ1,2=1.0021±0.2687i, as |arg(λ1,2)|=0.2620<π40. It indicates that the instability of the system may generate chaos. It can be seen from Figure 3 that the system is in a quasi-periodic state.

    Figure 3.  Numerical simulation of phase diagrams and time series plot of system (2.2) at (α1,α2)=(1.1,0.95)..

    To gain a clearer view of the system's dynamic behavior, Figure 4 shows the phase diagram of the system under different α orders.

    Figure 4.  Comparison of the phytoplankton-zooplankton phase diagram of the system (2.2) at different fractional-order derivatives.

    Next, we explore the potential for a Hopf bifurcation at E. A Hopf bifurcation is critical as it marks a transition in the model's stability and the emergence of periodic solutions. To achieve this, the characteristic roots of system (2.2) must be purely complex.

    The solution to system (2.2) is obtained as follows:

    λ1,2=Tr0±Tr204det02. (2.17)

    When α=1, system (2.2) undergoes a destabilizing Hopf bifurcation under the conditions tr0=0 and det0>0. Given that the stability of system (2.2) is influenced by the fractional derivative, this derivative can be considered a parameter in the Hopf bifurcation. We now determine the conditions for the Hopf bifurcation for system (2.2) around parameter E, α=αh:

    (1) At the equilibrium point E, the Jacobian matrix possesses a pair of complex conjugate eigenvalues λ1,2=ai+ibi, which transition to being purely imaginary at α=αh;

    (2) m(αh)=0 where m(α)=απ2min1i2|arg(λi)|;

    (3) m(α)α|α=αh0.

    We now demonstrate that E experiences a Hopf bifurcation as α crosses the value αh.

    Proof. Given that tr204det0<0 and tr0>0, the eigenvalues form a complex conjugate pair with a positive real component. Thus,

    0<arg(λ12)=tan1(4det0tr20tr0)<π2, (2.18)

    and απ2>|tan1(4det0tr20tr0)| for some α. Let αhπ2=|tan1(4det0tr20tr0)|, get αh=2πtan1(4det0tr20tr0). Moreover, m(α)α|α=αh=π20. Therefore, all Hopf conditions satisfy.

    Above, we analyze the dynamical behavior of fractional system (2.2). The expression of the solution of system (2.1) involves more complex special functions, such as the Mittag-Leffler function and other hypergeometric functions. Although the form of the solution exists, its expression is very large. Therefore, in this section, we reveal various spatiotemporal behaviors around the Turing bifurcation threshold d2 and derive the amplitude equations associated with system (2.1) at α1=α2=1 using multi-scale and weakly nonlinear analysis. The solution of system (2.1) is written in the following form:

    (ˉuˉv)=(ˉuˉv)+3j=1(AˉujAˉvj)eiqjr+c.c., (3.1)

    where (Aˉuj,Aˉvj)T represents the magnitude of the wave vector qj, fulfilling the condition that |qj|=qT.

    Adding a perturbation to the equilibrium point E(u,v), such that ˉu=ˉu+u, ˉv=ˉv+v, and substituting them into system (2.1), and performing a Taylor expansion, we can get the following form:

    {ut=a11u+a12v+k20u2+k02v2+k11uv+k30u3+k03v3+k21u2v+k12uv2+d12u,vt=a11u+a22v+m20u2+m02v2+m11uv+m30u3+m03v3+m21u2v+m12uv2+d22u+d32v, (3.2)

    where

    k02=0,k03=0,k12=0,m02=0,m03=0,m12=0,k11=1+αζ(1+αζ+ˉu)2,m11=(1+αζζ)τ(1+αζ+ˉu)2μG(G+ˉu)2,k21=1+αζ(1+αζ+ˉu)3,m21=τ(1+αζζ)(1+αζ+ˉu)3+μG(G+ˉu)3,k30=(1+αζ)(1ˉuκ)(1+αζ+ˉu)3,m30=τ(1ˉuκ)(1+αζζ)(1+αζ+ˉu)3μGˉv(G+ˉu)4,k20=1κ+(1+αζ)(1ˉuκ)(1+αζ+ˉu)2,m20=τ(1+αζζ)(1ˉuκ)(1+αζ+ˉu)2+μGˉv(G+ˉu)3.

    Let U=(u,v)T and system (3.2) can be simplified as

    Ut=LU+N(U,U), (3.3)
    L=(a11+d12a12a21+d22a22+d32),

    where

    N=(k20u2+k02v2+k11uv+k30u3+k03v3+k21u2v+k12uv2m20u2+m02v2+m11uc+m30u3+m03v3+m21u2v+m12uv2).

    Here, L denotes a linear operator and N represents a nonlinear operator.

    We expand parameter d2 with a sufficiently small parameter ε, so we can obtain:

    dT2=εd(1)2+ε2d(2)2+ε3d(3)2+O(ε3). (3.4)

    We develop variable U and nonlinear term N, respectively, the small parameter ε:

    U=(uv)=ε(u1v1)+ε2(u2v2)+ε3(u3v3)+O(ε3), (3.5)
    N=ε2N2+ε3N3+O(ε4), (3.6)

    where

    N2=(k20u21+k11u1v1m20u21+m11u1v1),N3=(2k20u1u2+k11u1v2+k11u2v1+k30u31+k21u21v12m20u1u2+m11u1v2+m11u2v1+m30u31+m21u21v1).

    We decompose operator L into the following form:

    L=Lc+(d2dT2)M, (3.7)

    where

    Lc=(d12+a11a12dT22+a21d32+a22),M=(b11b12b21b22).

    Using the multiple-scale approach, we differentiate the system's time scales into T1=εt,T2=ε2t,T3=ε3t, which are mutually independent. Consequently, the time derivative can be articulated as:

    t=εT1+ε2T2+ε3T3+O(ε3). (3.8)

    Substituting Eqs (3.4)–(3.8) into Eq (3.3), we can obtain the following equation:

    ε:Lc(u1v1)=0, (3.9)
    ε2:Lc(u2v2)=T1(u1v1)d(1)2M(u1v1)N2, (3.10)
    ε3:Lc(u3v3)=T1(u2v2)+T2(u1v1)d(1)2M(u2v2)d(2)2M(u1v1)N3. (3.11)

    The 1st order of (O(ε)):

    Lc(u1v1)=0. (3.12)

    The solution of Eq (3.12) is

    (u1v1)=(φ1)3j=1(Wjeiqjr)+c.c., (3.13)

    where φ=a12a11d1q2c,|qj|=qc,qc=qT(dT2). Additionally, Wj denotes the amplitude of eiqjr.

    The 2nd order of (O(ε2)):

    Lc(u2v2)=T1(u1v1)d(1)2M(u1v1)N2=(FuFv), (3.14)

    Fju and Fjv(j=1,2,3) are the coefficients with respect to eiqjr in Fu and Fv, respectively.

    Using the Fredholm solvability condition for Eq (3.14), the vector function on the right-hand side must be perpendicular to the zero eigenvalue of L+c. The zero eigenvector of Lc is:

    (1ψ)eiqjr+c.c., j=1,2,3, (3.15)

    with ψ=a11d1q2ca21dT2q2c. From the orthogonal condition:

    (1,ψ)(FjuFjv)=0,j=1,2,3. (3.16)

    Utilizing the Fredholm solvability condition as stated in Eq (3.16), we arrive at the following:

    {(φ+ψ)W1T1=d(1)2η0W1+2(η1+ψη2)¯W2¯W3,(φ+ψ)W2T1=d(1)2η0W2+2(η1+ψη2)¯W1¯W3,(φ+ψ)W3T1=d(1)2η0W3+2(η1+ψη2)¯W1¯W2, (3.17)

    where

    {η0=(φb11+b22)+ψ(φb21+b22),η1=12k20φ2+k11φ,η2=12m20φ2+m11φ.

    Next, higher-order disturbance terms are introduced:

    (u2v2)=(U0V0)+3j=1(UjVj)eiqjr+3j=1(UjjVjj)e2iqjr+(U12V12)ei(q1q2)r+(U23V23)ei(q2q3)r+(U31V31)ei(q3q1)r+c.c.. (3.18)

    We substitute formulas (3.13) and (3.18) into formula (3.10). By solving the equations corresponding to different modes, it is known that the coefficients have the following forms:

    (U0V0)=(u00v00)(|W1|2+|W2|2+|W3|2), Uj=φYj,(UjjVjj)=(u11v11)W2j,(UijVij)=(ˉuˉv)Wi¯Wj,

    where

    (ˉu00ˉv00)=2a11a22a21a12(a22η1a12η2a11η2a21η1), (3.19)
    (ˉu11ˉv11)=(a12η2(a224d3q2c)η1(a114d1q2c)(a224d3q2c)(a214dT2q2c)a12(a214dT2q2c)η1(a114d1q2c)η2(a114d1q2c)(a224d3q2c)(a214dT2q2c)a12), (3.20)
    (ˉuˉv)=2(a12η2(a223d3q2c)η1(a113d1q2c)(a223d3q2c)(a213dT2q2c)a12(a213dT2q2c)η1(a113d1q2c)η2(a113d1q2c)(a223d3q2c)(a213dT2q2c)a12). (3.21)

    From the orthogonal condition, we obtain:

    {(φ+ψ)(Y1T1+W1T2)=η0(d(2)2W1+d(1)2Y1)[(Q1+ψR1)|W1|2+(Q2+ψR2)(|W2|2+|W3|2)]W1+2(η1+ψη2)(¯W2¯Y3+¯W3¯Y2),(φ+ψ)(Y2T1+W2T2)=η0(d(2)2W2+d(1)2Y2)[(Q1+ψR1)|W2|2+(Q2+ψR2)(|W1|2+|W3|2)]W2+2(η1+ψη2)(¯W1¯Y3+¯W3¯Y1),(φ+ψ)(Y3T1+W3T2)=η0(d(2)2W3+d(1)2Y3)[(Q1+ψR1)|W3|2+(Q2+ψR2)(|W2|2+|W1|2)]W3+2(η1+ψη2)(¯W2¯Y1+¯W1¯Y2), (3.22)

    where

    Q1=(2φk20+k11)(u00+u11)φk11(v00+v11)3k30φ33k21φ2,Q2=(2φk20+k11)(u00+u)φk11(v00+v)6k30φ36k21φ2,R1=(2φm20+m11)(u00+u11)φm11(v00+v11)3m30φ33m21φ2,R2=(2φm20+m11)(u00+u)φm11(v00+v)6m30φ36m21φ2.

    The amplitude Aj is the coefficient of eiqjr at each level, so

    Aj=εWj+ε2Yj+O(ε3). (3.23)

    Substituting Eqs (3.17) and (3.22) into Eq (3.23), we get the following amplitude equations:

    {τ0Z1t=ηZ1+χˉZ2ˉZ3[g1|Z1|2+g2(|Z2|2+|Z3|2)]Z1,τ0Z2t=ηZ2+χˉZ1ˉZ3[g1|Z2|2+g2(|Z1|2+|Z3|2)]Z2,τ0Z3t=ηZ3+χˉZ1ˉZ2[g1|Z3|2+g2(|Z1|2+|Z2|2)]Z3, (3.24)

    where

    η=d2dT2dT2,τ0=φ+ψdT2[(φb11+b12)+ψ(φb21+b22)],χ=k20φ2+2k11φ+m20φ2ψ+2m11φψdT2[(φb11+b12)+ψ(φb21+b22)],g1=Q1+ψR1dT2[(φb11+b12)+ψ(φb21+b22)],g2=Q2+ψR2dT2[(φb11+b12)+ψ(φb21+b22)].

    Since each amplitude Aj=ωjeiθj(j=1,2,3) in Eq (3.24) can be decomposed into mode ωj=|Aj| and phase angle θj, substituting Aj into Eq (3.24) to separate the real and imaginary parts yields the following equation:

    {τ0θt=χω21ω22+ω21ω23+ω22ω23ω1ω2ω3sinθ,τ0ω1t=ηω1+χω2ω3cosθ[g1ω31+g2(ω22+ω23)ω1],τ0ω1t=ηω1+χω2ω3cosθ[g1ω31+g2(ω22+ω23)ω1],τ0ω1t=ηω1+χω2ω3cosθ[g1ω31+g2(ω22+ω23)ω1], (3.25)

    where, θ=θ1+θ2+θ3. From Eq (3.25), it can be deduced that the solution is stable under the conditions χ>0,ψ=0 and χ<0,ψ=π. The solutions of Eq (3.25) are shown in Table 4.

    Table 4.  The relationship between pattern shape and steady-state solutions.
    Conditions Solution Pattern shape
    ω1=ω2=ω3=0 Stationary state
    ω1=ηg1,ω2=ω3=0 Strip pattern
    η>η1=χ24(g1+2g2) ω1=ω2=ω3=|χ|±χ2+4(g1+2g2)η2(g1+2g2) Hexagon pattern
    g2>g1,η>g1ω21 ω1=|χ|g2g1,ω2=ω3=ηg1ω21g1+g2 Mixed state

     | Show Table
    DownLoad: CSV

    In this section, we conduct numerical simulations to verify the theoretical analysis and observe its pattern dynamic behavior. In order to observe the dynamic behavior, we introduce the numerical method of the PZEM model (2.1). If 0<αi<1(i=1,2), the time derivative is expressed by formula (2.11), and the space derivative is expressed by the following formula (4.2). If αi=1(i=1,2), we employ the Euler discretization approach for conducting numerical simulations in a two-dimensional domain Ω=[0,Lx]×[0,Ly]. We select Lx=250,Ly=250, a time increment Δt=0.2, and a spatial increment Δh=0.69. We denote unpq=u(xp,yq,nΔt) and vnpq=v(xp,yq,nΔt). The model (2.1) is discretized using the Euler method, as outlined below:

    {ˉun+1pqˉunpqΔt=ˉunpq(1ˉunpqκ)ˉunpqˉvnpq1+αζ+ˉunpq+d12ˉunpq,ˉvn+1pqˉvnpqΔt=τ(ˉunpq+ζ)ˉvnpq1+αζ+ˉunpqmˉvnpqμˉunpqˉvnpqG+ˉunpq+d22ˉunpq+d32ˉvnpq, (4.1)

    where

    {2ˉupq=ˉup+1,q+1+ˉup1,q1+ˉup+1,q1+ˉup1,q+1+4(ˉup+1,q+ˉup1,q+ˉup,q+1+ˉup,q1)20ˉupq6h2,2ˉvpq=ˉvp+1,q+1+ˉvp1,q1+ˉvp+1,q1+ˉvp1,q+1+4(ˉvp+1,q+ˉvp1,q+ˉvp,q+1+ˉvp,q1)20ˉvpq6h2. (4.2)

    Next, we select the parameters shown in Table 5 and use Eq (4.1) to conduct numerical simulations.

    Table 5.  Parameter values.
    α τ ζ G n κ m κc d1 d2 d3
    0.55 4.1 0.72 8.43 0.79 2.38 2.45 1.91 0.0036 0.00158 0.0141

     | Show Table
    DownLoad: CSV

    The parameter values are shown in Table 5, and the calculation results are as follows:

    E=(0.3131,1.4843),χ=0.4312,η=0.2461,τ0=0.3270,g1=1.2869,g2=1.7464.

    The initial conditions are specified as follows:

    u(x,y,0)=u(1+0.1(rand0.5)),v(x,y,0)=v(1+0.1(rand0.5)).

    The outcomes of the numerical simulation indicate that speckle and mixed-structure patterns emerge in the graphical representation with this particular set of parameters, as illustrated in Figure 5.

    Figure 5.  Density distribution pattern of system (2.1) at α=0.55,τ=4.1,ζ=0.72,G=8.43,n=0.79,κ=2.38,m=2.45,d1=0.0036,d2=0.00158,d3=0.0141, and α1=α2=1.

    Now, we use the parameters in Table 5 and symmetric initial conditions to perform numerical simulations by changing parameters d1, d2, and d3. Numerical simulation was carried out with other parameters remaining unchanged, and the results showed that there is a symmetric mixed structure solution. The subsequent figure illustrates the intricate pattern evolution of the symmetric mixed pattern across parameters. We select the following symmetric initial conditions:

    u(x,y,0)={u,x,y(80,120),u0.001,other,v(x,y,0)={v,x,y(80,120),v0.001,other.

    Figure 6 shows the spatial distribution patterns of phytoplankton and zooplankton at different diffusion coefficients, d2. The diffusion coefficient d reflects the influence of phytoplankton on the distribution of zooplankton. When d2=0.003, phytoplankton and zooplankton form highly clustered patterns, with zooplankton closely following the prey. When d2=0.005, phytoplankton and zooplankton are more dispersed, and the zooplankton pattern is more complex. This indicates that the diffusion ability of phytoplankton affects the hunting efficiency of zooplankton, and zooplankton tend to migrate to areas with a high density of phytoplankton. When d2=0.004, the distribution densities of phytoplankton and zooplankton are between the distribution densities of phytoplankton and zooplankton at d2=0.003 and d2=0.005.

    Figure 6.  Population density distribution pattern of phytoplankton and zootoplankton with with different parameters d1=0.002 and d3=0.018.

    Let d1=0.002,d2=0.004, and d3=0.015 through the observation of the distribution pattern of plankton population density in Figure 7. At t=6000, the system presents an incomplete pattern structure. Phytoplankton begins to form an irregular pattern distribution, while zooplankton shows an aggregation pattern following predation. At t=8000, the pattern structure gradually forms a stable state. At t=10000, the system reaches a stable state and forms a predator aggregation area that is misaligned with the plant pattern. We find that as time goes by, t, the distribution pattern of plankton population density becomes more complex. This indicates that the system presents a complex competitive process over time. This reflects the process by which the ecosystem reaches a stable state through dynamic adjustment.

    Figure 7.  The distribution patterns of phytoplankton and zooplankton population densities at different time periods.

    Let parameter d1=0.002,d2=0.004, and d3=0.017, the patterns of phytoplankton and zooplankton at different times are shown in Figure 8a1a4 and Figure 8c1c4. For diffusion coefficients at d1=0.002,d2=0.006, and d3=0.017, the patterns of phytoplankton and zooplankton at different times are shown in Figure 8b1b4 and Figure 8d1d4. In Figure 8, the density of phytoplankton gradually increases over time. It can be seen from the pattern that the difference between the high-density area and the low-density area significantly increases. The internal structure of the pattern becomes more complex, showing more details and layers. The distribution pattern of zooplankton shows a significantly misaligned predation aggregation area with that of plants, and the density of predation hotspots continues to increase over time. It can be seen from the pattern that the spatial distribution of phytoplankton and zooplankton shows complex dynamic behaviors.

    Figure 8.  Population density distribution patterns of phytoplankton and zootoplankton with parameters d1=0.002 and d3=0.017.

    Through the analysis of the patterns of phytoplankton and zooplankton, we see that, as time goes by, the patterns of zooplankton evolve from simple patterns to complex ones, demonstrating the complexity of this ecosystem. The diffusion coefficient is highly sensitive to the generation of patterns. A smaller diffusion coefficient leads to more regular and symmetrical patterns, while a larger diffusion coefficient results in more complex and irregular patterns, significantly affecting the distribution pattern of plankton. The initial conditions and diffusion coefficient significantly impact the long-term ecological behavior of the model. Moreover, initial conditions affect the formation of the pattern, while the parameters affect the distribution structure of the ecosystem.

    We combine theoretical analysis with numerical simulation to deeply explore the dynamic behavioral characteristics of the fractional-order PZEM. In theoretical analysis, we conduct stability, Turing instability, Hopf bifurcation, and nonlinear analysis for the PZEM. In numerical methods, we develop a high-precision numerical method for fractional-order PZEM without diffusion terms and verify the superiority of this method through comparison with other methods. Moreover, an effective discretization method is established for the model with diffusion terms. The numerical simulation results not only verify the correctness of the theoretical analysis, but also demonstrate complex dynamic behaviors at different diffusion coefficients. Moreover, the system presents significantly different spatial pattern evolution characteristics. The results reveal that the diffusion coefficient is crucial to the generation of patterns. A bigger diffusion coefficient will lead to more complex and irregular pattern structures, which directly affects the spatial distribution patterns of plankton populations. Furthermore, the initial conditions and the diffusion coefficient have a decisive influence on the long-term ecological behavior of the model: The initial conditions affect the formation process of the pattern, while the diffusion coefficient affects the dynamic evolution of the pattern structure. In future work, we will explore the response mechanism of the system at the coupling effect of multiple factors, as well as the effect of dynamic behavior for fractional-order.

    Conceptualization, Methodology, Software, Data, Formal analysis and Funding acquisition, Writing-original draft and writing-review and editing: Shuai Zhang, Hao Lu Zhang, Yu Lan Wang and Zhi Yuan Li. All authors have read and agreed to the published version of the manuscript.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This paper is supported by Doctoral research start-up Fund of Inner Mongolia University of Technology (DC2300001252) and the Natural Science Foundation of Inner Mongolia (2024LHMS06025).

    The authors declare that there are no conflicts of interest regarding the publication of this article.



    [1] J. C. Macdonald, H. Gulbudak, Forward hysteresis and Hopf bifurcation in an Npzd model with application to harmful algal blooms, J. Math. Biol., 87 (2023), 45. https://doi.org/10.1007/s00285-023-01969-7 doi: 10.1007/s00285-023-01969-7
    [2] H. Liu, C. J. Dai, H. G. Yu, Q. Guo, J. B. Li, A. M. Hao, et al. Dynamics of a stochastic non-autonomous phytoplankton-zooplankton system involving toxin-producing phytoplankton and impulsive perturbations, Math. Comput. Simul., 203 (2023), 368–386. https://doi.org/10.1016/j.matcom.2022.06.012 doi: 10.1016/j.matcom.2022.06.012
    [3] K. Dehingia, A. Das, E. Hincal, K. Hosseini, The effect of time delay on the dynamics of a plankton-nutrient system with refuge, Braz. J. Phys., 55 (2025), 28. https://doi.org/10.1007/s13538-024-01670-0 doi: 10.1007/s13538-024-01670-0
    [4] M. Xu, Y. Liu, Z. Zhao, K. Fu, X. Lv, Advancing parameter estimation with characteristic finite difference method (CFDM) for a marine ecosystem model by assimilating satellite observations: Spatial distributions, Front. Mar. Sci., 9 (2022), 997537. https://doi.org/10.3389/fmars.2022.997537 doi: 10.3389/fmars.2022.997537
    [5] A. W. Omta, E. A. Heiny, H. Rajakaruna, D. Talmy, M. J. Follows, Trophic model closure influences ecosystem response to enrichment, Ecol. Modell., 475 (2023), 110183. https://doi.org/10.1016/j.ecolmodel.2022.110183 doi: 10.1016/j.ecolmodel.2022.110183
    [6] M. N. Huda, Q. Q. A'yun, S. Wigantono, H. Sandariria, I. Raming, A. Asmaidi, Effects of harvesting and planktivorous fish on bioeconomic phytoplankton-zooplankton models with ratio-dependent response functions and time delays, Chaos Soliton Fract., 173 (2023), 113736. https://doi.org/10.1016/j.chaos.2023.113736 doi: 10.1016/j.chaos.2023.113736
    [7] K. Chakraborty, K. Das, Modeling and analysis of a two-zooplankton one-phytoplankton system in the presence of toxicity, Appl. Math. Model., 39 (2015), 1241–1265. https://doi.org/10.1016/j.apm.2014.08.004 doi: 10.1016/j.apm.2014.08.004
    [8] Y. Lv, J. Cao, J. Song, R. Yuan, Y. Pei, Global stability and Hopf-bifurcation in a zooplankton-phytoplankton model, Nonlinear Dyn., 76 (2014), 345–366. https://doi.org/10.1007/s11071-013-1130-2 doi: 10.1007/s11071-013-1130-2
    [9] V. P. Dubey, J. Singh, A. M. Alshehri, S. Dubey, D. Kumar, Numerical investigation of fractional model of phytoplankton-toxic phytoplankton-zooplankton system with convergence analysis, Int. J. Biomath., 15 (2022), 2250006. https://doi.org/10.1142/S1793524522500061 doi: 10.1142/S1793524522500061
    [10] P. Panja, T. Kar, D. K. Jana, Impacts of global warming on phytoplankton-zooplankton dynamics: A modelling study, Environ. Dev. Sustain., 26 (2024), 13495–13513. https://doi.org/10.1007/s10668-023-04430-3 doi: 10.1007/s10668-023-04430-3
    [11] I. M. Suthers, A. J. Richardson, D. Rissik, The importance of plankton, in Plankton: A Guide to Their Ecology and Monitoring for Water Quality, Australia: CSIRO Publishing, 2019, 1–13. https://doi.org/10.1080/00288330.2019.1625497
    [12] W. Wang, S. Liu, D. Tian, D. Wang, Pattern dynamics in a toxin-producing phytoplankton-zooplankton model with additional food, Nonlinear Dyn., 94 (2018), 211–228. https://doi.org/10.1007/s11071-018-4354-3 doi: 10.1007/s11071-018-4354-3
    [13] J. Alzabut, R. Dhineshbabu, A. G. M. Selvam, J. F. Gomez-Aguilar, H. Khan, Existence, uniqueness and synchronization of a fractional tumor growth model in discrete time with numerical results, Results Phys., 54 (2023), 107030. https://doi.org/10.1016/j.rinp.2023.107030 doi: 10.1016/j.rinp.2023.107030
    [14] W. AI-Sadi, Z. Wei, S. Q. T. Abdullah, A. Alkhazzan, J. F. Gomez-Aguilar, Dynamical and numerical analysis of the hepatitis B virus treatment model through fractal-fractional derivative, Math. Methods Appl. Sci., 48 (2024), 639–657. https://doi.org/10.1002/mma.10348 doi: 10.1002/mma.10348
    [15] M. S. Shabbir, Q. Din, M. De la Sen, J. F. Gomez-Aguilar, Exploring dynamics of plant-herbivore interactions: Bifurcation analysis and chaos control with Holling type-Ⅱ functional response, J. Math. Biol., 88 (2024), 8. https://doi.org/10.1007/s00285-023-02020-5 doi: 10.1007/s00285-023-02020-5
    [16] H. Khan, J. lzabut, J. F. Gomez-Aguilar, A. Alkhazan, Essential criteria for existence of solution of a modified-ABC fractional order smoking model, Ain Shams Eng. J., 15 (2024), 102646. https://doi.org/10.1016/j.asej.2024.102646 doi: 10.1016/j.asej.2024.102646
    [17] N. Raza, A. Raza, Y. Chahlaoui, J. F. Gomez-Aguilar, Numerical analysis of HPV and its association with cervical cancer using Atangana-Baleanu fractional derivative, Model. Earth Syst. Environ., 11 (2025), 60. https://doi.org/10.1007/s40808-024-02243-5 doi: 10.1007/s40808-024-02243-5
    [18] P. Li, R. Gao, C. Xu, Y. Li, A. Akgul, D. Baleanu, Dynamics exploration for a fractional-order delayed zooplankton-phytoplankton system. Chaos Soliton Fract., 166 (2023), 112975. https://doi.org/10.1016/j.chaos.2022.112975
    [19] M. Javidi, B. Ahmad, Dynamic analysis of time fractional order phytoplankton-toxic phytoplankton-zooplankton system, Ecol. Model., 318 (2015), 8–18. https://doi.org/10.1016/j.ecolmodel.2015.06.016 doi: 10.1016/j.ecolmodel.2015.06.016
    [20] P. Kumar, V. Suat Erturk, R. Banerjee, M. Yavuz, V. Govindaraj, Fractional modeling of plankton-oxygen dynamics under climate change by the application of a recent numerical algorithm, Phys. Scr., 96 (2021), 124044. https://doi.org/10.1088/1402-4896/ac2da7 doi: 10.1088/1402-4896/ac2da7
    [21] C. P. Li, Z. Q. Li, Z. Wang, Mathematical analysis and the local discontinuous Galerkin method for Caputo-Hadamard fractional partial differential equation, J. Sci. Comput., 85 (2020), 41. https://doi.org/10.1007/s10915-020-01353-3 doi: 10.1007/s10915-020-01353-3
    [22] C. P. Li, Z. Q. Li, The blow-up and global existence of solution to Caputo-Hadamard fractional partial differential equation with fractional Laplacian, J. Nonlinear Sci., 31 (2021), 80. https://doi.org/10.1007/s00332-021-09736-y doi: 10.1007/s00332-021-09736-y
    [23] C. Han, Y. L. Wang, Z. Y. Li, Novel patterns in a class of fractional reaction-diffusion models with the Riesz fractional derivative, Math. Comput. Simul., 202 (2022), 149–163. https://doi.org/10.1016/j.matcom.2022.05.037 doi: 10.1016/j.matcom.2022.05.037
    [24] C. Han, Y. L. Wang, Z. Y. Li, A high-precision numerical approach to solving space fractional Gray-Scott model, Appl. Math. Lett., 125 (2022), 107759. https://doi.org/10.1016/j.aml.2021.107759 doi: 10.1016/j.aml.2021.107759
    [25] C. P. Li, Q. Yi, A. Chen, Finite difference methods with non-uniform meshes for nonlinear fractional differential equations, J. Comput. Phys., 316 (2016), 614–631. https://doi.org/10.1016/j.jcp.2016.04.039 doi: 10.1016/j.jcp.2016.04.039
    [26] X. H. Wang, H. L. Zhang, Y. L. Wang, Z. Y. Li, Dynamic properties and numerical simulations of the fractional Hastings–Powell model with the Grünwald-Letnikov differential derivative, Int. J. Bifurcation Chaos, 35 (2025), In press.
    [27] M. Singh, S. Das, Rajeev, S. H. Ong, Novel operational matrix method for the numerical solution of nonlinear reaction-advection-diffusion equation of fractional order, Comput. Appl. Math., 41 (2022), 306. https://doi.org/10.1007/s40314-022-02017-8 doi: 10.1007/s40314-022-02017-8
    [28] M. Kashif, M. Singh, T. Som, E. M. Craciun, Numerical study of variable order model arising in chemical processes using operational matrix and collocation method, J. Comput. Sci., 80 (2024), 102339. https://doi.org/10.1016/j.jocs.2024.102339 doi: 10.1016/j.jocs.2024.102339
    [29] P. Roul, Design and analysis of efficient computational techniques for solving a temporal-fractional partial differential equation with the weakly singular solution, Math. Methods Appl. Sci., 47 (2024), 2226–2249. https://doi.org/10.1002/mma.9744 doi: 10.1002/mma.9744
    [30] T. Kumari, P. Roul, An efficient computational technique for solving a time-fractional reaction-subdiffusion model in 2D space, Comput. Math. Appl., 160 (2024), 191–208. https://doi.org/10.1016/j.camwa.2024.02.018 doi: 10.1016/j.camwa.2024.02.018
    [31] Anjuman, A. Y. T. Leung, S. Das, Two-dimensional time-fractional nonlinear drift reaction-diffusion equation arising in electrical field, Fractal Fract., 8 (2024), 456. https://doi.org/10.3390/fractalfract8080456 doi: 10.3390/fractalfract8080456
    [32] R. Sharma, A numerical study for nonlinear time-space fractional reaction-diffusion model of fourth-order, J. Comput. Nonlinear Dyn., 20 (2025), 021001. https://doi.org/10.1115/1.4067065 doi: 10.1115/1.4067065
    [33] K. M. Owolabi, S. Jain, E. Pindza, E. Mare, Comprehensive numerical analysis of time-fractional reaction-diffusion models with applications to chemical and biological phenomena, Mathematics, 12 (2024), 3251. https://doi.org/10.3390/math12203251 doi: 10.3390/math12203251
    [34] M. Sohaib, K. M. Furati, A. Shah, Space fractional Allen-Cahn equation and its applications in phase separation: A numerical study, Commun. Nonlinear Sci. Numer. Simul., 137 (2024), 108173. https://doi.org/10.1016/j.cnsns.2024.108173 doi: 10.1016/j.cnsns.2024.108173
    [35] I. Petráš, Fractional-Order Nonlinear Systems: Modeling, Analysis and Simulation, Heidelberg: Springer Science & Business Media, 2011. https://doi.org/10.1007/978-3-642-18101-6
    [36] X. L. Gao, H. L. Zhang, X. Y. Li, Research on pattern dynamics of a class of predator-prey model with interval biological coefficients for capture, , AIMS Math., 9 (2024), 18506–18527. https://doi.org/10.3934/math.2024901 doi: 10.3934/math.2024901
    [37] D. Y. Xue, L. Bai, Numerical algorithms for Caputo fractional-order differential equations, Int. J. Control, 90 (2016), 1201–1211. https://doi.org/10.1080/00207179.2016.1158419 doi: 10.1080/00207179.2016.1158419
    [38] D. Y. Xue, C. N. Zhao, Y. Q. Chen, A modified approximation method of fractional order system, in 2006 International Conference on Mechatronics and Automation, IEEE, Luoyang, (2006), 1043–1048. https://doi.org/10.1109/ICMA.2006.257769
    [39] D. Y. Xue, Fractional Calculus and Fractional-Order Control, Beijing: Science Press, 2018.
    [40] W. H. Deng, J. H. Lü, Design of multidirectional multiscroll chaotic attractors based on fractional differential systems via switching control, Chaos, 16 (2006), 043120. https://doi.org/10.1063/1.2401061 doi: 10.1063/1.2401061
  • This article has been cited by:

    1. XiongFei Chi, Ting Wang, ZhiYuan Li, Chaotic dynamics analysis and control of fractional-order ESER systems using a high-precision Caputo derivative approach, 2025, 33, 2688-1594, 3613, 10.3934/era.2025161
    2. Xian-He Lei, Yutong Cui, Harmonic balance method for MEMS oscillators, 2025, 13, 2296-424X, 10.3389/fphy.2025.1597421
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(165) PDF downloads(24) Cited by(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog