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

Deep belief improved bidirectional LSTM for multivariate time series forecasting

  • Received: 28 June 2023 Revised: 04 August 2023 Accepted: 07 August 2023 Published: 17 August 2023
  • Multivariate time series (MTS) play essential roles in daily life because most real-world time series datasets are multivariate and rich in time-dependent information. Traditional forecasting methods for MTS are time-consuming and filled with complicated limitations. One efficient method being explored within the dynamical systems is the extended short-term memory networks (LSTMs). However, existing MTS models only partially use the hidden spatial relationship as effectively as LSTMs. Shallow LSTMs are inadequate in extracting features from high-dimensional MTS; however, the multilayer bidirectional LSTM (BiLSTM) can learn more MTS features in both directions. This study tries to generate a novel and improved BiLSTM network (DBI-BiLSTM) based on a deep belief network (DBN), bidirectional propagation technique, and a chained structure. The deep structures are constructed by a DBN layer and multiple stacked BiLSTM layers, which increase the feature representation of DBI-BiLSTM and allow for the model to further learn the extended features in two directions. First, the input is processed by DBN to obtain comprehensive features. Then, the known features, divided into clusters based on a global sensitivity analysis method, are used as the inputs of every BiLSTM layer. Meanwhile, the previous outputs of the shallow layer are combined with the clustered features to reconstitute new input signals for the next deep layer. Four experimental real-world time series datasets illustrate our one-step-ahead prediction performance. The simulating results confirm that the DBI-BiLSTM not only outperforms the traditional shallow artificial neural networks (ANNs), deep LSTMs, and some recently improved LSTMs, but also learns more features of the MTS data. As compared with conventional LSTM, the percentage improvement of DBI-BiLSTM on the four MTS datasets is 85.41, 75.47, 61.66 and 30.72%, respectively.

    Citation: Keruo Jiang, Zhen Huang, Xinyan Zhou, Chudong Tong, Minjie Zhu, Heshan Wang. Deep belief improved bidirectional LSTM for multivariate time series forecasting[J]. Mathematical Biosciences and Engineering, 2023, 20(9): 16596-16627. doi: 10.3934/mbe.2023739

    Related Papers:

    [1] Yuanshi Wang, Donald L. DeAngelis . A mutualism-parasitism system modeling host and parasite with mutualism at low density. Mathematical Biosciences and Engineering, 2012, 9(2): 431-444. doi: 10.3934/mbe.2012.9.431
    [2] Jordi Ripoll, Jordi Font . Numerical approach to an age-structured Lotka-Volterra model. Mathematical Biosciences and Engineering, 2023, 20(9): 15603-15622. doi: 10.3934/mbe.2023696
    [3] Suqing Lin, Zhengyi Lu . Permanence for two-species Lotka-Volterra systems with delays. Mathematical Biosciences and Engineering, 2006, 3(1): 137-144. doi: 10.3934/mbe.2006.3.137
    [4] Wenrui Li, Qimin Zhang, Meyer-Baese Anke, Ming Ye, Yan Li . Taylor approximation of the solution of age-dependent stochastic delay population equations with Ornstein-Uhlenbeck process and Poisson jumps. Mathematical Biosciences and Engineering, 2020, 17(3): 2650-2675. doi: 10.3934/mbe.2020145
    [5] Guichen Lu, Zhengyi Lu . Permanence for two-species Lotka-Volterra cooperative systems with delays. Mathematical Biosciences and Engineering, 2008, 5(3): 477-484. doi: 10.3934/mbe.2008.5.477
    [6] Junjing Xiong, Xiong Li, Hao Wang . The survival analysis of a stochastic Lotka-Volterra competition model with a coexistence equilibrium. Mathematical Biosciences and Engineering, 2019, 16(4): 2717-2737. doi: 10.3934/mbe.2019135
    [7] Xixia Ma, Rongsong Liu, Liming Cai . Stability of traveling wave solutions for a nonlocal Lotka-Volterra model. Mathematical Biosciences and Engineering, 2024, 21(1): 444-473. doi: 10.3934/mbe.2024020
    [8] Harindri Chaudhary, Mohammad Sajid, Santosh Kaushik, Ali Allahem . Stability analysis of chaotic generalized Lotka-Volterra system via active compound difference anti-synchronization method. Mathematical Biosciences and Engineering, 2023, 20(5): 9410-9422. doi: 10.3934/mbe.2023413
    [9] Tuan A. Phan, Benjamin J. Ridenhour, Christopher H. Remien . Resilience of a stochastic generalized Lotka–Volterra model for microbiome studies. Mathematical Biosciences and Engineering, 2025, 22(6): 1517-1550. doi: 10.3934/mbe.2025056
    [10] Xiaomeng Ma, Zhanbing Bai, Sujing Sun . Stability and bifurcation control for a fractional-order chemostat model with time delays and incommensurate orders. Mathematical Biosciences and Engineering, 2023, 20(1): 437-455. doi: 10.3934/mbe.2023020
  • Multivariate time series (MTS) play essential roles in daily life because most real-world time series datasets are multivariate and rich in time-dependent information. Traditional forecasting methods for MTS are time-consuming and filled with complicated limitations. One efficient method being explored within the dynamical systems is the extended short-term memory networks (LSTMs). However, existing MTS models only partially use the hidden spatial relationship as effectively as LSTMs. Shallow LSTMs are inadequate in extracting features from high-dimensional MTS; however, the multilayer bidirectional LSTM (BiLSTM) can learn more MTS features in both directions. This study tries to generate a novel and improved BiLSTM network (DBI-BiLSTM) based on a deep belief network (DBN), bidirectional propagation technique, and a chained structure. The deep structures are constructed by a DBN layer and multiple stacked BiLSTM layers, which increase the feature representation of DBI-BiLSTM and allow for the model to further learn the extended features in two directions. First, the input is processed by DBN to obtain comprehensive features. Then, the known features, divided into clusters based on a global sensitivity analysis method, are used as the inputs of every BiLSTM layer. Meanwhile, the previous outputs of the shallow layer are combined with the clustered features to reconstitute new input signals for the next deep layer. Four experimental real-world time series datasets illustrate our one-step-ahead prediction performance. The simulating results confirm that the DBI-BiLSTM not only outperforms the traditional shallow artificial neural networks (ANNs), deep LSTMs, and some recently improved LSTMs, but also learns more features of the MTS data. As compared with conventional LSTM, the percentage improvement of DBI-BiLSTM on the four MTS datasets is 85.41, 75.47, 61.66 and 30.72%, respectively.



    The Lotka-Volterra (LV) system is often applied to depict the evolutionary process in population dynamic, physics, and economics[1,2,9,13]. Particularly, a cooperative LV system explains the interactions and benefits together among species in social animals and in human society, such as bees and flowers, algal-fungal associations of lichens, etc [4,5,6,7]. In literature, many studies focused on the deterministic cooperative LV system [8,9,10,11,12,13]. In [9], the authors considered the permanence in a two-species delay LV cooperative system. Also, Lu et al. [10] showed the certain and general delay effects on the performance of a LV cooperative system. Furthermore, in [11], the authors considered a non-autonomous cooperative LV system and established sufficient conditions to keep the system permanent. Xu et al. [12] studied the discrete LV cooperative system with periodic boundary conditions from the aspect of stability. Since environmental noise is nonnegligible in the ecosystem, some scholars introduced the stochastic noise into the cooperative LV model and revealed the effects of stochastic noise on the model (see e.g. [2,3]). In particular, Zuo et al. [2] introduced the white noise to the intrinsic growth rates of the certain delay cooperative LV model, then they established the sufficient conditions for the persistence and obtained the existence and uniqueness of the stationary distribution. In [3], Liu studied a stochastic food-limited LV cooperative model and formulated the sufficient conditions of p-moment persistence and extinction. However, all these works mentioned above mainly focused on the stochastic ordinary differential equations (SODEs) rather than the stochastic partial differential equations (SPDEs).

    Age-dependent theory of populations was first studied by Lotka in 1907 [14]. From then on, age-dependent has been introduced in many models to describe the dynamical behavior of species, such as the predator-prey LV model, population model and epidemic model [16,17,20,21,22,23,24,25,26,31,32,37,38]. Zhang and Liu [16] investigated the age-dependent in the predator-prey model and analysed the non-trivial periodic oscillation phenomenon. Solis and Ku-Carrillo [17] studied a predator-prey LV model with age-dependent and presented the theoretical results which preserve the properties of the solutions of the model. Chekroun and Kuniya [34] concerned with the global asymptotic behavior of an infection age-structured SIR epidemic model with diffusion in a general n-dimensional bounded spatial domain. Then they showed that the basic reproduction R0 depended on the shape of the spatial domain in the 2-dimensional case. Lu and Wei [35] formulated a stochastic epidemic model with age of vaccination and adopted a generalized incidence rate to make the model more realistic. They got the sufficient conditions for extinction and two types of permanence. Yang and Wang [36] introduced an age-dependent predation system with prey harvesting and provided the explicit formulae which can determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions.

    However, as far as we know, most of the studies in literature focused on the predator-prey LV and population models instead of cooperative LV system. Based on the importance of the cooperative LV system and to fill the gap of lacks of studying the SPDEs in these systems, in this study, we introduce the age-dependent into a stochastic cooperative LV system and study the properties of this system.

    It is extremely difficult to get the explicit solutions of the stochastic age-dependent system, so in this study, we mainly focus on the numerical solutions of the system. The numerical method we use in this study is the Euler-Maruyama (EM) method. There are various advantages of the EM method, chief among them are the simple algebraic expression, low computational cost, and good convergence rate, etc [25,26,27,28]. However, the EM method usually requires the system to have Lipschitz continuous and linear growth condition on both the drift and diffusion coefficients. Since our system does not meet the Lipschitz continuous condition, we could not apply the EM algorithm directly. Instead, we use the theory of equilibrium point to avoid the restrictions on the coefficients and apply the EM method to the stochastic age-dependent cooperative LV system skilfully to obtain the approximate solution. This is one of the innovations of this paper.

    We provide below a brief summary of, and comments on, our results.

    ● We introduce the age-dependent and stochastic white noise into a cooperative LV system and study the existence and uniqueness of the solutions. In particular, we have proven that a unique positive global solution exists for this stochastic age-dependent cooperative LV system;

    ● Using the positive equilibrium, we apply EM method to the stochastic age-dependent cooperative LV system in an appropriate region and study the numerical solution of the system;

    ● The pth-moment boundedness and strong convergence of the EM scheme for the stochastic age-dependent cooperative LV system under certain conditions are proved;

    ● Numerical experiments are presented, and they support our theoretical results.

    The construction of this paper is designed as follows. In Section 2, we establish the stochastic age-dependent cooperative LV system and present some basic preliminaries. In Section 3, we introduce the EM approximation method to the stochastic age-dependent cooperative LV system. The pth-moment estimations of this scheme have been provided. Section 4 shows that the scheme is convergent and the strong convergence order is provided. In Section 5, some numerical examples are presented to verify our theories in Section 4. Section 6 gives the conclusions of this paper.

    The typical cooperative LV system can be described in the form:

    {dx(t)dt=x(t)(α11x(t)+α12y(t)+r1),dy(t)dt=y(t)(α21x(t)α22y(t)+r2), (2.1)

    where x(t), y(t) represent the density of the two cooperative species, r1,r2 are the intrinsic growth rates of the two species, α11, α22 denote the intraspecific competition rates, α12,α21 are the interspecific cooperation rates. All the rates are assumed to be positive constants [2].

    Since birth rate is sensitive to age, according to [18,19], we introduce the age-dependent fertility rate into Eq (2.1) and obtain the age-dependent cooperative LV system as follows:

    {X(t,a)t+X(t,a)a=X(t,a)(α11(a)X(t,a)+α12(a)Y(t,a)+r1(a)),(t,a)Q,Y(t,a)t+Y(t,a)a=Y(t,a)(α21(a)X(t,a)α22(a)Y(t,a)+r2(a)),(t,a)Q,X(t,0)=A0γ(t,a)X(t,a)da,t[0,T],Y(t,0)=A0β(t,a)Y(t,a)da,t[0,T],X(0,a)=X0(a), Y(0,a)=Y0(a),a[0,A], (2.2)

    where X(t,a) and Y(t,a) denote the density of the two species at time t, age a. Q=(0,T)×(0,A).

    Now, we consider the effect of random environment noise in Eq (2.2) and assume ri(a)ri(a)+σi˙ω(t), (i=1,2), where w(t) is a standard Brownian motion. So, we have the stochastic age-dependent cooperative LV system.

    {dtX(t,a)=X(t,a)adt+X(t,a)(α11(a)X(t,a)+α12(a)Y(t,a)+r1(a))dt+σ1X(t,a)dw(t),(t,a)Q,dtY(t,a)=Y(t,a)adt+Y(t,a)(α21(a)X(t,a)α22(a)Y(t,a)+r2(a))dt+σ2Y(t,a)dw(t),(t,a)Q,X(t,0)=A0γ(t,a)X(t,a)da,t[0,T],Y(t,0)=A0β(t,a)Y(t,a)da,t[0,T],X(0,a)=X0(a), Y(0,a)=Y0(a),a[0,A], (2.3)

    where dtX and dtY are the differential of X(t,a) and Y(t,a) relative to t, i.e., dtX=Xtdt, and dtY=Ytdt. Features of the parameters in Eq (2.3) are showed in Table 1.

    Table 1.  Symbols and their meaning in Eq (2.3).
    Parameters Meanings
    ri(a) the intrinsic growth rates of the species at age a, i=1,2
    αi,i(a) the intraspecific competition rates at age a, i=1,2
    α12(a) the interspecific cooperation rate at age a
    α21(a) the interspecific cooperation rate at age a
    γ(t,a) the fertility rate of females of X(t,a) at age a
    β(t,a) the fertility rate of females of Y(t,a) at age a
    σi the size of uncertainty from the population growth, i=1,2

     | Show Table
    DownLoad: CSV

    The integral form of Eq (2.3) is given by

    {X(t)=X0t0X(s)ads+t0X(s)(α11(a)X(s)+α12(a)Y(s)+r1(a))ds+t0σ1X(s)dw(s),Y(t)=Y0t0Y(s)ads+t0Y(s)(α21(a)X(s)α22(a)Y(s)+r2(a))ds+t0σ2Y(s)dw(s),X(t,0)=A0γ(t,a)X(t,a)da,t[0,T],Y(t,0)=A0β(t,a)Y(t,a)da,t[0,T], (2.4)

    where X(t):=X(t,a), Y(t):=Y(t,a), X0:=X(0,a), and Y0:=Y(0,a).

    Remark 2.1. Compared with Eq (2.1), the advantages of Eq (2.3) are summarized as follows:

    The intrinsic growth rates r1,r2 in Eq (2.1) are constant. However, the intrinsic growth rate should not be a certain number due to the random noise in the natural environment. Therefore, we describe the influence of the intrinsic growth rate of the environmental noise through a standard Brownian motion, that is, ri(a)ri(a)+σi˙ω(t),(i=1,2) in Eq (2.3).

    The density of the species is not only affected by interspecific relationships, but also by the reproductive capacity of females. Eq (2.1) does not take the latter into account. So we introduce the female fertility into the Eq (2.3) and can be expressed as X(t,0)=A0γ(t,a)X(t,a)da,Y(t,0)=A0β(t,a)Y(t,a)da, which makes the Eq (2.3) more practical than model Eq (2.1).

    The density of the two species are the function of time t in model Eq (2.1). While in Eq (2.3), the density of the two species are the binary function of time t and age a. This is because the species density are not only affected by the time, but also affected by the change of birth rate and death rate which are caused by the age a.

    Throughout the paper, let (Ω,F,{Ft}t0,P) be a complete probability space. Here, the filtration {Ft}t0 satisfies the usual conditions (that is, it is increasing and right continuous with F0 containing all P-null sets), E denotes the expectation corresponding to P. For a set A, its indicator function is denoted by 1A={1,xA,0,xA. We also denote by Rn+ the positive cone in Rn, that is Rn+={xRn:xi>0 for all 1in}. If xRn, its norms are denoted as |x|ι={|x1|+|x2|++|xn|,ι=1,(ni=1x2i)12,ι=2.

    Let V{φ|φLp([0,A]),φaLp([0,A]), where φa is generalized partial derivatives with respect to age a}. H=L2([0,A]) satisfy that VHHV, where V is the dual space of V, H is the dual space of H. The norm in H is denoted as ||. The duality product between V,V is written as ,, the scalar product in H is denoted by (,). C=C([0,T];H) is the space of all continuous functions from [0,T] into H. Ip([0,T];V) denotes the space of all V-valued processes (Pt)t[0,T], LpV=Lp([0,T];V). W:=(Ip([0,T];V)L2(Ω;C([0,T];H)))×(Ip([0,T];V)L2(Ω;C([0,T];H))).

    In order to analyze the stochastic age-dependent cooperative LV system Eq (2.3), we make the following basic assumptions:

    (A1) limαAAαγ(t,a)X(t,a)da=λx>0; limαAAαβ(t,a)Y(t,a)da=λy>0.

    (A2) The fertility rate of females γ(t,a),β(t,a)C([0,T]×[0,A];H).

    (A3) αij(a)C([0,A];R+),i,j{1,2} and satisfy supa[0,A]α12(a)α21(a)<infa[0,A]α11(a)α22(a).

    (A4) ri(a)C([0,A];R) and ri(a) is nondecreasing with <ri(0)<0<ri(A)<, for i=1,2.

    (A5) supa[0,A]α12(a)α21(a)<1 and r1(0)(1α12α21)+α12(r2(A)+α21r1(A))<0; r2(0)(1α12α21)+α21(r1(A)+α12r2(A))<0.

    Remark 2.2. The biological background significance of assumptions (A1)-(A5) are described as follows:

    (A1) When the age a approaches to point A for species X, the number of new individuals born at time t is λx, which is a positive constant. This means that the species X at age A has not lost its reproductive capacity.

    (A2) γ(t,a) and β(t,a) are the fertility rate of females in species X and Y, respectively. From biological point of view, the fertility rate of females in X and Y are both continuous functions, i.e., there is no sudden external noise that may increase or decrease the fertility rate of the species within a short time.

    (A3) The product of interspecific cooperation rate is less than the intraspecific competition rate. This assumption is mainly used to prove the existence of global positive solution of the system.

    (A4) ri(0) represents the intrinsic growth rate at age 0. "ri(0)<0" means that the birth rate is less than the death rate. "0<ri(A)<" means that the birth rate is greater than the death rate at age A. Combined with the actual biological background, "0<ri(A)<" also indicates that the species will not overreproduce at age A.

    (A5) The value of α12 and α21 are estimated by the actual data based on the least-square method. From biological point of view, this means that the intrinsic growth rate of species X at age 0 is less than α12(r2(A)+α21r1(A))1α12α21.

    Theorem 2.1. Under the assumptions (A1)–(A3), for any given initial value (X0(a0),Y0(a0))R2+, there exists a unique global solution (X(t),Y(t)) of Eq (2.3) on W. Moreover, the solution is nonnegative with probability 1.

    Proof. It is easy to see that for any given initial value (X0(a0),Y0(a0))R2+, there is a unique local solution (X(t,a),Y(t,a))W when t[0,τe), where τe is the explosion time of Eq (2.3). Now we need to proof that τe=, a.s. Choosing a positive constant k0>1 such that for any initial value (X0(a0),Y0(a0)), we have (X0(a0),Y0(a0))[1k0,k0]×[1k0,k0]. Then, for each kk0, a[0,A], we define the stopping time τk as τk=inf{t[0,τe):(X,Y)(1k,k)×(1k,k)}.

    We can deduce that τk is increasing when k and limkτk=τ with ττe. Assuming that there are two constants T>0 and ε(0,1) such that P{τT}>ε, then, for the constant k0, there is an integer k1 such that

    P{τkT}ε, for all kk1. (2.5)

    Now, for any positive constants cx,cy, we define a function V:WR+ as

    V(X,Y)=cx(XlogX1)+cy(YlogY1).

    we can check that V(X,Y)0. By the Itô formula [15],

    LV(X,Y)=cx11X,Xa+X(α11X+α12Y+r1)+12(cxσ21+cyσ22)+cy11Y,Ya+Y(α21Xα22Y+r2)Ψdiag(cx,cy)ˉRˉC(AΨ+ˉR)+12Ψ(diag(cx,cy)A+Adiag(cx,cy))Ψ+cxA01XdaX+cyA01YdaY+12(cxσ21+cyσ22)cxlogλx+cylogλy[cxlog|A0γ(t,a)X(t,a)da|+cylog|A0β(t,a)Y(t,a)da|]+12(cxσ21+cyσ22)+Ψdiag(cx,cy)ˉRˉC(AΨ+ˉR)+12Ψ(diag(cx,cy)A+Adiag(cx,cy))Ψ

    where Ψ=(X(t,a)Y(t,a)), ˉR=(r1(a)r2(a)), diag(cx,cy)=(cx00cy), A=(α11(a)α12(a)α21(a)α22(a)), ˉC=(cxcy).

    By assumption (A3), since A is a nonsingular M-matrix, we deduce

    LV(X,Y)cxlogλx+cylogλy[cxlog|A0γ(t,a)X(t,a)da|+cylog|A0β(t,a)Y(t,a)da|]+12(cxσ21+cyσ22)+Ψdiag(cx,cy)ˉRˉC(AΨ+ˉR)cxlogλx+cylogλy[cxlog|A0γ(t,a)X(t,a)da|+cylog|A0β(t,a)Y(t,a)da|]+{ˉRdiag(cx,cy)ˉCA}Ψ+12(cxσ21+cyσ22)r1cxr2cy.

    According to assumptions (A1) and (A2), we have

    LV(x,y)|ˉRdiag(cx,cy)ˉCA||Ψ|+CM1(1+|Ψ|). (2.6)

    Now, let M2=cxcycxcy,M3=cxcy, then

    M2V(X,Y)=cxcycxcy[cx(X1logX)+cy(Y1logY)]V(X,Y).

    Since |Ψ|=(X2+Y2)12X+Y, we have

    |Ψ|2(X1logX)+2(Y1logY)+44+2M3V(X,Y). (2.7)

    By Eqs (2.6) and (2.7), we derive that

    LV(X,Y)M4(1+V(X,Y)),

    where M4=M1(52M3). Then, we have

    EV(X(tτk),Y(tτk))V(X0,Y0)+Etτk0M4(1+V(X(s),Y(s)))dsM5+M4t0EV(X(sτk),Y(sτk))ds,

    where M5=V(X0,Y0)+M4T.

    Applying the Gronwall inequality, we obtain

    EV(X(τkT),Y(τkT))M5eM4T.

    Since P(Ωk)=P{τkT}ε for kk1, for every ωΩk,(t,a)Q, X(τk,ω) and Y(τk,ω) equal either k or 1k, hence

    V(X(τk,ω),Y(τk,ω))[cx(k1logk)+cy(k1logk)][cx(1k1+logk)+cy(1k1+logk)]

    and

    M5eM4TEV(X(τkT),Y(τkT))E[IΩk(ω)V(X(τkT),Y(τkT))]ε[cx(k1logk)+cy(k1logk)][cx(1k1+logk)+cy(1k1+logk)].

    Taking k, we have the contradiction that

    >M5eM4T.

    The proof is completed.

    In this part, we apply the EM approximate method to the stochastic age-dependent cooperative LV Equation (2.3) and explore the moment boundedness of the scheme for Eq (2.3).

    We begin with introducing the following notations in the Eq (2.3).

    G1x(X,Y):=σ1X(t,a), G1y(X,Y):=σ2Y(t,a),
    H1x(X,Y):=r1(a)X(t,a), H2x(X,Y):=X(t,a)[α11(a)X(t,a)+α12(a)Y(t,a)],
    H1y(X,Y):=r2(a)Y(t,a), H2y(X,Y):=Y(t,a)[α21(a)X(t,a)α22(a)Y(t,a)].

    Equation (2.3) has four equilibria: E0(0,0), E1(r1(A)α11,0), E2(0,r2(A)α22), and E(x,y), where

    x=α22(a)r1(A)+α12(a)r2(A)α11(a)α22(a)α12(a)α21(a), y=α11(a)r2(A)+α21(a)r1(A)α11(a)α22(a)α12(a)α21(a).

    Now, let us first present some lemmas.

    Remark 3.1. We note Hix(Ψj):=Hix(Xj,Yj),Hiy(Ψj):=Hiy(Xj,Yj),Hix(ψj):=Hix(xj,yj),Hiy(ψj):=Hiy(xj,yj), (i=1,2),(j=1,2) as follows.

    Lemma 3.1. Under assumptions (A4) and (A5), for any ψ1,ψ2W, we have

    |H1x(ψ1)H1x(ψ2)||G1x(ψ1)G1x(ψ2)|L1|ψ1ψ2|, (3.1)
    |H1y(ψ1)H1y(ψ2)||G1y(ψ1)G1y(ψ2)|L2|ψ1ψ2|, (3.2)
    |H2x(ψ1)H2x(ψ2)|ρ1|ψ1ψ2|1212ρ1|ψ1ψ2|2, (3.3)

    and

    |H2y(ψ1)H2y(ψ2)|ρ2|ψ1ψ2|1212ρ2|ψ1ψ2|2, (3.4)

    where L1, L2 are constants and ρ1=2xr1(0)+α12(x+y), ρ2=2yr2(0)+α21(x+y), ψi(t,a)=(xi(t,a),yi(t,a)), (i=1,2).

    The proof of this lemma is similar to that in Yang et al. [30]. We skip it here.

    Remark 3.2. From Lemma 3.1, we can easily derive that there exist constants K1 and K2 such that

    |H1x(ψ)||G1x(ψ)|K1(1+|ψ|), (3.5)
    |H1y(ψ)||G1y(ψ)|K2(1+|ψ|), (3.6)

    and

    xH2x(ψ)212ρ1(1+|ψ|2), yH2y(ψ)212ρ2(1+|ψ|2), (3.7)

    where ψ(t,a)=(x(t,a),y(t,a)) is the numerical solution corresponding to the exact solution Ψ of Equation (2.3), H1x(ψ)=H1x(x,y),H1x(ψ)=H1x(x,y)

    Lemma 3.2. For any p>2, we have the following inequalities.

    x[H1x(ψ)+H2x(ψ)]+p12|G1x(ψ)|2[2K1+K21(p1)+212ρ1](1+|ψ|2),

    and

    y[H1y(ψ)+H2y(ψ)]+p12|G1y(ψ)|2[2K2+K22(p1)+212ρ2](1+|ψ|2).

    Proof. From Eqs (3.1), (3.5) and (3.7), we have

    x[H1x(ψ)+H2x(ψ)]+p12|G1x(ψ)|2K1|x|(1+|ψ|)+p12K21(1+|ψ|)2+212ρ1(1+|ψ|2)K1|ψ|(1+|ψ|)+(p1)K21(1+|ψ|2)+212ρ1(1+|ψ|2)K1|ψ|+K1|ψ|2+K21(p1)(1+|ψ|2)+212ρ1(1+|ψ|2),

    by the definition of |ψ| and the positive of |ψ|2|ψ|+1, we derive that

    x[H1x(ψ)+H2x(ψ)]+p12|G1x(ψ)|22K1(1+|ψ|2)+K21(p1)(1+|ψ|2)+212ρ1(1+|ψ|2)=[K21(p1)+2K1+212ρ1](1+|ψ|2).

    Similarly, using Eqs (3.2), (3.6) and (3.7), we have

    y[H1y(ψ)+H2y(ψ)]+p12|G1y(ψ)|2K2|ψ|(1+|ψ|)+p12K22(1+|ψ|)2+212ρ2(1+|ψ|2)K2|ψ|+K2|ψ|2+K22(p1)(1+|ψ|2)+212ρ2(1+|ψ|2)[K22(p1)+2K2+212ρ2](1+|ψ|2).

    This completes the proof.

    Now, we investigate the approximate solution of Eq (2.3) by using the EM method. Given the fixed time T and time step Δ(0,1), let tk=kΔ for k=0,1,2,,[TΔ], where [TΔ] denotes the integer part of TΔ. So, the discrete time EM approximate solution to Eq (2.3) is defined as

    {xk+1=xk+{xka+r1(a)xkα11(a)xkxk+α12(a)xkyk}Δ+σ1xkΔwk,yk+1=yk+{yka+r2(a)yk+α21(a)xkykα22(a)ykyk}Δ+σ2ykΔwk,

    where Δwk=wtk+1wtk, X(0)=X0(a), and Y(0)=Y0(a).

    The corresponding continuous EM approximate solution of Eq (2.3) is formed as

    {x(t)=x0t0ˉx(s)ads+t0{r1(a)ˉx(s)α11(a)ˉx(s)ˉx(s)+α12(a)ˉx(s)ˉy(s)}ds+t0σ1ˉx(s)dw(s),y(t)=y0t0ˉy(s)ads+t0{r2(a)ˉy(s)+α21(a)ˉx(s)ˉy(s)α22(a)ˉy(s)ˉy(s)}ds+t0σ2ˉy(s)dw(s), (3.8)

    where

    ˉx(t)=[TΔ]k=0xk1[tk,tk+1)(t), ˉy(t)=[TΔ]k=0yk1[tk,tk+1)(t)

    are called the step processes. We can easily see that ˉx(t)=xk and ˉy(t)=yk for t[tk,tk+1) when k=0,1,2,,[TΔ]. Without loss of generality, we denote ψ(t)=(x(t),y(t)) and ˉψ(t)=(ˉx(t),ˉy(t)).

    Theorem 3.1. Under assumptions (A1)–(A5) and p>2, there exists a constant C>0 such that

    sup0tTE|ψ(t)|pC,

    where C only dependents on T and p.

    Proof. For any p>2,t(0,T), applying the Itô formula to the first equation of (3.8), we have

    E|x(t)|p|x0|pEt0p|x(s)|p2x(s)a,x(s)ds+Et0p(p1)2|x(s)|p2|G1x(ˉx(s),ˉy(s))|2ds+Et0p|x(s)|p2{|x(s)|[H1x(ˉx(s),ˉy(s))+H2x(ˉx(s),ˉy(s))]}dsEt0p|x(s)|p2{|ˉx(s)|[H1x(ˉx(s),ˉy(s))+H2x(ˉx(s),ˉy(s))]+p12|G1x(ˉx(s),ˉy(s))|2}ds+Et0p|x(s)|p2{|x(s)ˉx(s)|[H1x(ˉx(s),ˉy(s))+H2x(ˉx(s),ˉy(s))]}ds+Et0p|x(s)|p2x(s)a,x(s)ds.

    Since

    x(s)a,x(s)=A0x(s)dax(s)=12[A0γ(s,a)x(s,a)da]212A0γ2(s,a)daA0x2(s,a)da12Aˇγ2(s,a)|x(s,a)|2,

    where ˇγ means the maximum value of continuous function γ(t,a). Then we have

    E|x(t)|p|x0|pEt0p|x(s)|p2{|ˉx(s)|[H1x(ˉx(s),ˉy(s))+H2x(ˉx(s),ˉy(s))]+p12|G1x(ˉx(s),ˉy(s))|2}ds+Et0p|x(s)|p2{|x(s)ˉx(s)|[H1x(ˉx(s),ˉy(s))+H2x(ˉx(s),ˉy(s))]}ds+12Aˇγ2pEt0|x(s)|pds. (3.9)

    Let

    Q1=Et0p|x(s)|p2{|ˉx(s)|[H1x(ˉx(s),ˉy(s))+H2x(ˉx(s),ˉy(s))]+p12|G1x(ˉx(s),ˉy(s))|2}ds,
    Q2=Et0p|x(s)|p2|x(s)ˉx(s)||H1x(ˉx(s),ˉy(s))|ds,
    Q3=Et0p|x(s)|p2|x(s)ˉx(s)||H2x(ˉx(s),ˉy(s))|ds.

    By Lemma 3.2 and Young inequality aβ+b1ββa+(1β)b, for a,b0 and β(0,1), we have

    Q1(p2)[K21(p1)+2K1+212ρ1]Et0|x(s)|pds+2[K21(p1)+2K1+212ρ1]Et0(1+|ˉψ(s)|2)p2ds(p2)[K21(p1)+2K1+212ρ1]Et0|x(s)|pds+2p[K21(p1)+2K1+212ρ1]Et0(1+|ˉψ(s)|p)ds2p[K21(p1)+2K1+212ρ1]T{1+t0E|ψ(s)|pds+t0E|ˉψ(s)|pds}. (3.10)

    Similarly, we can show that

    Q2C[1+t0(E|ψ(s)|p+E|ˉψ(s)|p)ds]. (3.11)

    Moreover, by the Young inequality and Eq (3.3), we have

    Q3212ρ1Et0{(p2)|x(s)|p+2|x(s)ˉx(s)|p2|ψ(s)|p2}ds212ρ1(p2)Et0|x(s)|pds+232ρ1Et0|x(s)ˉx(s)|p2|ψ(s)|p2ds. (3.12)

    On the other hand, there is a unique nonnegative constant k which satisfies that tkstk+1 for any s[0,T]. By Itô integral and the Hölder, Burkholder-Davis-Gundy (BDG) inequalities, we have

    E|x(s)ˉx(s)|p2=E|stk(H1x(ˉx(u),ˉy(u))+H2x(ˉx(u),ˉy(u)))du+stkG1x(ˉx(u),ˉy(u))dw(u)+stkˉx(u)adu|p23p21E|stk(H1x(ˉx(u),ˉy(u))+H2x(ˉx(u),ˉy(u)))du|p2+3p21E|stkˉx(u)adu|p2+3p21E|stkG1x(ˉx(u),ˉy(u))dw(u)|p23p21Δp21E|stk2p21(Hp21x(ˉx(u),ˉy(u))+Hp22x(ˉx(u),ˉy(u)))du|+3p21E|stkˉx(u)adu|p2+Cp3p21Kp21E|stk(1+|ˉψ(u)|)2du|p43p2Δp22Estk(Hp21x(ˉx(u),ˉy(u))+Hp22x(ˉx(u),ˉy(u)))du+CΔp2+Cp3p21Kp21E|stk2(1+|ˉψ(u)|2)du|p43p2Δp22Kp21Estk(1+|ˉψ(u)|)p2du+CΔp22Estkρp21|ˉψ(u)|p2du+CΔp2+Cp3p1Kp21(Δ+Estk|ˉψ(u)|2du)p4CΔp4+CΔp4E|ˉψ(s)|p2, (3.13)

    where C is a positive constant and Cp is defined by

    Cp={(64p)p4,0<p<4,4,p=4,[(p2)p2+12(p21)p21]p4,p>4.

    Substituting Eq (3.13) into Eq (3.12), we have

    Q3212ρ1Et0(p2)|x(s)|pds+232ρ1Et0|x(s)ˉx(s)|p2|ψ(s)|p2ds212ρ1(p2)Et0|x(s)|pds+232Cρ1Δp4t0(1+E|ˉψ(s)|p2)|ψ(s)|p2ds212ρ1(p2)t0E|x(s)|pds+252Cρ1Δp4t0E|ˉψ(s)|pds. (3.14)

    Now, substituting Eqs (3.10), (3.11) and (3.14) into Eq (3.9), we can derive that

    E|x(t)|p|x0|p+12A2ˇγ2pEt0|x(s)|pds+C[1+Et0|ψ(s)|pds+Et0|ˉψ(s)|pds]C(1+t0sup0υsE|ψ(υ)|pds). (3.15)

    Since the right-hand side of Eq (3.15) is non-decreasing for t[0,T], a conclusion can be draw that

    sup0υtE|x(υ)|pC(1+t0sup0υsE|ψ(υ)|pds). (3.16)

    Similarly, we can get that

    E|y(t)|p|y0|pEt0p|y(s)|p2y(s)a,y(s)ds+Et0p(p1)2|y(s)|p2|G1y(ˉx(s),ˉy(s))|2ds+Et0p|y(s)|p2[|y(s)|(H1y(ˉx(s),ˉy(s))+H2y(ˉx(s),ˉy(s)))]ds.

    Since

    y(s)a,y(s)=A0y(s)day(s)=12[A0β(s,a)y(s,a)da]212A0β2(s,a)daA0y2(s,a)da12Aˇβ2(s,a)|y(s,a)|2,

    where ˇβ is the maximum value of continuous function β(t,a), and

    E|y(t)|p|y0|pEt0p|y(s)|p2{|ˉy(s)|[H1y(ˉx(s),ˉy(s))+H2y(ˉx(s),ˉy(s))]+p12|G1y(ˉx(s),ˉy(s))|2}ds+Et0p|y(s)|p2{|y(s)ˉy(s)|[H1y(ˉx(s),ˉy(s))+H2y(ˉx(s),ˉy(s))]}ds+12Aˇβ2pEt0|y(s)|pds. (3.17)

    Let

    G1=Et0p|y(s)|p2{|ˉy(s)|[H1y(ˉx(s),ˉy(s))+H2y(ˉx(s),ˉy(s))]+p12σ22|ˉy(s)|2}ds,
    G2=Et0p|y(s)|p2|y(s)ˉy(s)||H1y(ˉx(s),ˉy(s))|ds,
    G3=Et0p|y(s)|p2|y(s)ˉy(s)||H2y(ˉx(s),ˉy(s))|ds.

    By Lemma 3.2, we have

    G1pEt0|y(s)|p2[K22(p1)+2K2+212ρ2](1+|ˉψ(s)|2)ds. (3.18)

    By Eq (3.6), we have

    G2pK2Et0|y(s)|p2|y(s)ˉy(s)|(1+|ψ(s)|)dsC[1+t0(E|ψ(s)|p+E|ˉψ(s)|p)ds]. (3.19)

    Moreover, using the Young inequality, we can get

    G3212pρ2Et0|y(s)|p2|y(s)ˉy(s)||ψ(s)|ds212ρ2Et0(p2)|y(s)|p+232ρ2Et0|y(s)ˉy(s)|p2|ψ(s)|p2ds. (3.20)

    Similarly, by the Hölder, BDG inequalities and the Itô integral, we derive that

    E|y(s)ˉy(s)|p2CΔp4+CΔp4E|ˉψ(u)|p2. (3.21)

    Substituting Eq (3.21) into Eq (3.20), we get

    G3212ρ2(p2)Et0|y(s)|pds+232ρ2CΔp4t0(1+E|ˉψ(s)|p2)|ψ(s)|p2ds. (3.22)

    According to Eqs (3.18), (3.19), (3.22), and Eq (3.17), we have

    E|y(t)|pC(1+t0sup0υsE|ψ(υ)|pds). (3.23)

    For any t[0,T], since the right-hand side of Eq (3.23) is non-decreasing function, we obtain

    sup0υtE|y(υ)|pC(1+t0sup0υsE|ψ(υ)|pds). (3.24)

    Using Eqs (3.16), (3.24) and the norm of ψ(t), we have

    sup0υtE|ψ(υ)|pC(1+t0sup0υsE|ψ(υ)|pds).

    Applying the continuous Gronwall inequality, it yields

    sup0tTE|ψ(t)|pC.

    Now, let us present a theorem which shows that ψ(t) and ˉψ(t) are coincide with the discrete solution.

    Theorem 3.2. Under assumptions (A1)–(A5) and p>2, there exists a constant C such that

    E|ψ(t)ˉψ(t)|p<CΔp2.

    Proof. For any t[0,T], there exists a positive integer k such that t[tk,tk+1], and

    E|x(t)ˉx(t)|p=E|ttkx(s)ads+ttk[H1x(ˉx(s),ˉy(s))+H2x(ˉx(s),ˉy(s))]ds+ttkG1x(ˉx(s),ˉy(s))dw(s)|p4p1(E|ttkx(s)ads|p+E|ttkH1x(ˉx(s),ˉy(s))ds|p+E|ttkH2x(ˉx(s),ˉy(s))ds|p+E|ttkG1x(ˉx(s),ˉy(s))dw(s)|p).

    According to the Hölder, BDG, and elementary inequalities, we have

    E|x(t)ˉx(t)|p4p1[Δp1Ettk|x(s)a|pds+Kp1Δp1Ettk(1+|ψ(s)|)pds+2p2ρp1Δp1Ettk|ψ(s)|pds+CpE|ttkG21x(s)ds|p2]4p1[CΔp+2p1Kp1Δp1(Δ+Ettk|ψ(s)|pds)+2p2ρp1Δp1Ettk|ψ(s)|pds+CpE|ttkG21x(s)ds|p2]4p1[(C+2p1Kp1+Cp2p1Kp1)Δp2+(2p1Kp1+2p2ρp1+Cp2p1Kp1)Δp2E|ψ(t)|p]CΔp2(1+E|ψ(t)|p), (3.25)

    where Cp is defined in the proof of Theorem 3.1. Similarly, we have

    E|y(t)ˉy(t)|p=E|ttky(s)ads+ttk[H1y(ˉx(s),ˉy(s))+H2y(ˉx(s),ˉy(s))]ds+ttkG1y(ˉx(s),ˉy(s))dw(s)|p4p1(E|ttky(s)ads|p+E|ttkH1y(ˉx(s),ˉy(s))ds|p+E|ttkH2y(ˉx(s),ˉy(s))ds|p+E|ttkG1y(ˉx(s),ˉy(s))dw(s)|p)4p1[Δp1Ettk|y(s)a|pds+Kp2Δp1Ettk(1+|ψ(s)|)pds+2p2ρp2Δp1Ettk|ψ(s)|pds+CpE|ttkG21y(s)ds|p2]4p1[CΔp+2p1Kp2Δp1(Δ+Ettk|ψ(s)|pds)+2p2ρp2Δp1Ettk|ψ(s)|pds+Cp2p1Kp2Δp2(1+E|ψ(t)|p)]CΔp2(1+E|ψ(t)|p). (3.26)

    Using Eqs (3.25) and (3.26), we deduce that

    E|ψ(t)ˉψ(t)|pCΔp2(1+E|ψ(t)|p).

    By Theorem 3.1, we have

    E|ψ(t)ˉψ(t)|pCΔp2.

    Theorem 3.2 shows that the numerical solution converges to the step process. Now, we can discuss the convergence order rate between the exact and numerical solution of Eq (2.3).

    Theorem 4.1. Under assumptions (A1)–(A5) and p>2, we have

    Esup0tT|Ψ(t)ψ(t)|pCΔp2.

    where C is a constant and Ψ(t):=(X(t),Y(t)) is the true solution of Eq (2.4).

    Proof. For any t[0,T], an application of Itô formula yields

    d|X(t)x(t)|p=p|X(t)x(t)|p2(X(t)x(t),H1x(Ψ(t))H1x(ˉψ(t)))dt+p|X(t)x(t)|p2(X(t)x(t),H2x(Ψ(t))H2x(ˉψ(t)))dt+p|X(t)x(t)|p2(X(t)x(t),[G1x(Ψ(t))G1x(ˉψ(t))]dw(t))+p(p1)2(|X(t)x(t)|p2,|G1x(Ψ(t)G1x(ˉψ(t)))|2)dtp|X(t)x(t)|p2X(t)x(t),(X(t)x(t))adtp(|X(t)x(t)|p1,[G1x(Ψ(t))G1x(ˉψ(t))]dw(t))+p(p1)2L21|X(t)x(t)|p2|Ψ(t)ˉψ(t)|2dtp|X(t)x(t)|p1,(X(t)x(t))adt+pρ1212|X(t)x(t)|p1|Ψ(t)ˉψ(t)|dt+pL1|X(t)x(t)|p1|Ψ(t)ˉψ(t)|dt,

    hence we have

    Esups[0,t]|X(s)x(s)|p(pL1+212pρ1)t0Esupr[0,s]|X(r)x(r)|p1|Ψ(r)ˉψ(r)|ds+pEs0sups[0,t](|X(s)x(s)|p1,[G1x(Ψ(s))G1x(ˉψ(s))]dw(s))+p(p1)2L21t0Esupr[0,s]|X(r)x(r)|p2|Ψ(r)ˉψ(r)|2ds+p2A2γ2(t,a)t0Esupr[0,s]|X(r)x(r)|pds(2pL1+232pρ1+L212p(p1))t0Esupr[0,s]|X(r)x(r)|p2|Ψ(r)ˉψ(r)|2ds+psups[0,t]s0E|X(s)x(s)|p2(X(s)x(s),[G1x(Ψ(s))G1x(ˉψ(s))]dw(s))+(p2A2γ2(t,a)+2pL1+232pρ1)t0Esupr[0,s]|X(r)x(r)|pds. (4.1)

    Using the Young inequality, we derive that

    t0Esupr[0,s]|X(r)x(r)|p2|Ψ(r)ˉψ(r)|2dst0Esupr[0,s]{(12p)|X(r)x(r)|p+2p|Ψ(r)ˉψ(r)|p}ds. (4.2)

    By the BDG inequality, Young inequality, and elementary inequality, we have

    Esups[0,t]s0((X(s)x(s))p1,[G1x(Ψ(s))G1x(ˉψ(s))]dw(s))C1E{t0|((X(s)x(s))p1,G1x(Ψ(s))G1x(ˉψ(s)))|2ds}12C1Esups[0,t]|X(s)x(s)|p1{t0|G1x(Ψ(s))G1x(ˉψ(s))|2ds}12C1Esups[0,t]p1p2|X(s)x(s)|p+C1Lp1p2pEsups[0,t]{t0|Ψ(s)ˉψ(s)|2ds}p2C1(p1)p2Esups[0,t]|X(s)x(s)|p+C1Lp1Tp22p2pEsups[0,t]t0|Ψ(s)ˉψ(s)|pds. (4.3)

    Now, substituting Eqs (4.2) and (4.3) into Eq (4.1), we can derive

    Esups[0,t]|X(s)x(s)|p{p2A2ˇγ2(t,a)+2pL1+232pρ1+(2pL1+232pρ1+L212p(p1))(12p)}t0Esupr[0,s]|X(r)x(r)|pds+{2pp(2pL1+232pρ1+L212p(p1))+C12p1Lp1pp1Tp22}Esups[0,t]t0|Ψ(s)ψ(s)|pds+{2pp(2pL1+232pρ1+L212p(p1))+C12p1Lp1pp1Tp22}Et0|ψ(s)ˉψ(s)|pds+p1pEsups[0,t]|X(s)x(s)|p. (4.4)

    Similarly, we can show

    Esups[0,t]|Y(s)y(s)|p(pL2+212ρ2p)t0Esupr[0,s]|Y(r)y(r)|p1|Ψ(r)ˉψ(r)|ds+pEsups[0,t]s0(|Y(s)y(s)|p1,[G1y(Ψ(s))G1y(ˉψ(s))]dw(s))+p(p1)2L22t0Esupr[0,s]|Y(r)y(r)|p2|Ψ(r)ˉψ(r)|2ds+p2A2ˇβ2(t,a)t0Esupr[0,s]|Y(r)y(r)|pds(2pL2+232pρ2+p(p1)2L22)t0Esupr[0,s]|Y(r)y(r)|p2|Ψ(r)ˉψ(r)|2ds+pEsups[0,t]s0(|Y(s)y(s)|p1,[G1y(Ψ(s))G1y(ˉψ(s))]dw(s))+(p2A2ˇβ2(t,a)+2pL2+232pρ2)t0Esupr[0,s]|Y(r)y(r)|pds. (4.5)

    Therefore, due to Young inequality,

    t0Esupr[0,s]|Y(r)y(r)|p2|Ψ(r)ˉψ(r)|2dst0Esupr[0,s][(12p)|Y(r)y(r)|p+2p|Ψ(r)ˉψ(r)|p]ds. (4.6)

    Moreover, applying the BDG inequality, Young inequality, and Hölder's inequality, it yields

    Esups[0,s]s0(|Y(s)y(s)|p1,[G1y(Ψ(s))G1y(ˉψ(s))]dw(s))p1p2Esups[0,t]|Y(s)y(s)|p+C2Esups[0,t][t0|G1y(Ψ(s))G1y(ˉψ(s))|2ds]p2p1p2Esups[0,t]|Y(s)y(s)|p+C2Lp2Tp22Esups[0,t]t0|Ψ(s)ˉψ(s)|pds. (4.7)

    Now, substituting Eqs (4.6) and (4.7) into Eq (4.5), we obtain

    1pEsups[0,t]|Y(s)y(s)|p[p2A2ˇβ2+2pL2+232pρ2+(2pL2+232pρ2+p(p1)2L22)(12p)]t0Esupr[0,s]|Y(r)y(r)|pds+[2pp(2pL2+232pρ2+p(p1)2L22)+2p1C2pTp22Lp2]Esups[0,t]t0|Ψ(s)ψ(s)|pds+[2pp(2pL2+232pρ2+p(p1)2L22)+2p1C2pTp22Lp2]Et0|ψ(s)ˉψ(s)|pds. (4.8)

    By applying Eqs (4.4) and (4.8), Theorem 3.2, and the Gronwall inequality, we have

    Esups[0,T]|Ψ(s)ψ(s)|pCΔp2.

    Therefore, the proof is complete.

    Corollary 4.2. Assume the assumptions (A1)–(A5) and p>2 hold, the numerical solution of Eq (2.3) will converge to the exact solution in the sense

    limΔt0E(sup0tT|Ψ(t)ψ(t)|p)=0.

    Remark 4.1. Theorem 4.1 and Corollary 4.2 show that the exact solution Ψ(t) and the numerical solution ψ(t) are close to each other, which means that the numerical scheme constructed in this paper is effective.

    In this part, in order to simulate the numerical approximation of the EM scheme, we consider the following stochastic age-dependent cooperative LV system.

    {dtX=Xadt+X[1(1a)2X+cos2aY2]dt12dw(t),(t,a)Q,dtY=Yadt+Y[sin2aXe1aY+3a2]dt12dw(t),(t,a)Q,X(t,0)=1011aX(t,a)da,t[0,1],Y(t,0)=1011aY(t,a)da,t[0,1],X(0,a)=e21a,Y(0,a)=e21a,a[0,1], (5.1)

    where α11(a)=1(1a)2, α12(a)=cos2a, α21(a)=sin2a, α22(a)=e1a, γ(t,a)=β(t,a)=11a, Q=(0,1)×(0,1). All computations are run with T=1,Δ=0.005, and Δa=0.05 for Eq (5.1). In Equation (5.1), we simulate the EM numerical solution of the two species respectively (see Figure 1). Figure 1 shows that both the two species are bounded, which is consistent with Theorem 3.1. From the range of pictures (a) and (b) in Figure 1, we conclude that the fluctuation of species Y(t,a) is larger than that of X(t,a), which indicates that the number of species X(t,a) is relatively stable.

    Figure 1.  The EM numerical solutions of X(t,a) and Y(t,a) of Eq (5.1).

    Now, we apply the exact solution of Eq (2.3) by the method [29] to demonstrate the effectiveness of the EM scheme. In Figure 2, we present the exact solutions of X(t,a) and Y(t,a) with perturbation respectively.

    Figure 2.  The exact solutions of X(t,a) and Y(t,a) with perturbation of Eq (5.1).

    Next, in order to test and verify the convergence of the EM method in the stochastic age-dependent model, we conduct the error analysis of Eq (5.1). Figure 3 shows the errors between the exact solution and numerical solution of species X(t,a). It is easy to observe that the maximum value of the squared error is below 0.2, which means that the numerical solution tends to the exact solution.

    Figure 3.  The error simulation of species X(t,a) of Eq (5.1).

    Similarly, we get the error estimate of species Y(t,a). From Figure 4, it shows that the square error of Y(t,a) is less than 0.05. This confirms our theoretical results in Theorem 4.1.

    Figure 4.  The error simulation of species Y(t,a) of Eq (5.1).

    To study the convergence of EM algorithm in Eq (5.1) more precisely, we calculate the ι-th order error between x(t,a) and X(t,a). From Table 2, we observe that under the same value of order ι, the error (X(t,a)x(t,a))ι becomes smaller as step size Δ gets smaller. In addition, we find that increasing the order ι can reduce the error of X(t,a) and x(t,a) under the same step size. Similarly, Table 3 illustrates the same phenomena.

    Table 2.  Error simulation between X(t,a) and x(t,a).
    (X(t,a)x(t,a))ι Δ=0.01 Δ=0.005 Δ=0.002 Δ=0.001
    Error Time/s Error Time/s Error Time/s Error Time/s
    ι=4 1.2×102 26.881276 1.5×102 63.584526 5×103 373.050756 2×103 984.649714
    ι=6 3×103 22.164523 1.8×103 76.028217 4×104 1040.348042 2×104 1323.191920
    ι=8 1.6×104 23.557376 1.2×104 40.686763 1×104 588.323881 2×105 1570.733100
    ι=10 1.2×104 22.441634 0.5×104 63.584526 1×105 681.486228 1×106 1608.228890

     | Show Table
    DownLoad: CSV
    Table 3.  Error simulation between Y(t,a) and y(t,a).
    (Y(t,a)y(t,a))ι Δ=0.01 Δ=0.005 Δ=0.002 Δ=0.001
    Error Time/s Error Time/s Error Time/s Error Time/s
    ι=4 1×102 41.756866 5×103 81.872743 4×103 646.558116 2.5×103 2773.836731
    ι=6 6×104 42.402488 1.5×104 76.028217 6×105 587.713828 2×105 2851.378369
    ι=8 2×105 42.937773 8×106 80.371846 4×106 680.575583 2×106 1570.733100
    ι=10 1.2×106 42.679908 1×106 82.256243 3×107 992.707197 1×107 1608.228890

     | Show Table
    DownLoad: CSV

    Inspired by [33], in order to examine the stability of the numerical scheme for Eq (5.1), we present 3-dimensional and 2-dimensional trajectory numerical solution respectively. The paths of x(t,a) and y(t,a) for Eq (5.1) are given in Figure 5(a) and Figure 5(c). The two figures show that the limit of x(t,a) and y(t,a) are 0. In order to study the relationship of t and x(t,a), we fix the age a[0,1], then obtain the 2-dimensional trajectory numerical solution of x(t,a). Similarly, we could obtain the trajectory of y(t,a). Figure 5 shows that the numerical scheme is asymptotically stable, (a.s.).

    Figure 5.  The stability of numerical solution for Eq (5.1).

    In this paper, the stochastic age-dependent cooperative LV system has been established. Since it is extremely difficult to obtain the explicit solution for this system, the main contribution of this article is to investigate a numerical scheme to approximate the exact solution of the stochastic age-dependent cooperative LV system. By using the theory of M-matrices, we obtain the sufficient conditions which ensure the stochastic age-dependent cooperative LV Eq (2.3) has a global unique positive solution (see Theorem (2.1)). Since the coefficients of Eq (2.3) do not satisfy the Lipschitz and linear growth conditions, we can not apply the EM scheme directly. However, we apply the EM scheme to Eq (2.3) in a suitable region successfully through studying the equilibrium points. After that, we show that the approximation algorithm constructed in this paper is p-th moments boundedness (see Theorem 3.1) and it converges to the exact solution in the strong sense (see Theorem 4.1). Lastly, we present some numerical experiments which support our theoretical results. In the future, we will construct new algorithms which have a higher convergence order.

    The research was supported in part by the Natural Science Foundation of China (11661064), the Basic Research Project of North Minzu University (2018XYSYK01) and the Natural Science Foundation of Ningxia Province (2020AAC03065).

    The authors declare that they have no competing interests.



    [1] Y. Liu, H. Yang, S. Gong, Y. Liu, X. Xiong, A daily activity feature extraction approach based on time series of sensor events, Math. Biosci. Eng., 17 (2020), 5173–5189. https://doi.org/ 0.3934/mbe.2020280
    [2] H. Li, J. Tong, A novel clustering algorithm for time-series data based on precise correlation coefficient matching in the IoT, Math. Biosci. Eng., 16 (2019), 6654–6671. https://doi.org/10.3934/mbe.2019331 doi: 10.3934/mbe.2019331
    [3] H. M. Srivastava, I. C. Area Carracedo, J. L. Nieto, Power-series solution of compartmental epidemiological models, Math. Biosci. Eng., 18 (2021), 3274–3290. https://doi.org/10.3934/mbe.2021163 doi: 10.3934/mbe.2021163
    [4] M. Li, S. Chen, X. Chen, Y. Zhang, Y. Wang, Q. Tian, Symbiotic graph neural networks for 3d skeleton-based human action recognition and motion prediction, IEEE Trans. Pattern Anal. Mach. Intell., 44 (2021), 3316–3333. https://doi.org/10.1109/TPAMI.2021.3053765 doi: 10.1109/TPAMI.2021.3053765
    [5] M. Gan, Y. Cheng, K. Liu, G. Zhang, Seasonal and trend time series forecasting based on a quasi-linear autoregressive model, Appl. Soft Comput., 24 (2014), 13–18. https://doi.org/10.1016/j.asoc.2014.06.047 doi: 10.1016/j.asoc.2014.06.047
    [6] J. Wang, S. Zhang, An improved deep learning approach based on exponential moving average algorithm for atrial fibrillation signals identification, Neurocomputing, 513 (2013), 127–136. https://doi.org/10.1016/j.neucom.2022.09.079 doi: 10.1016/j.neucom.2022.09.079
    [7] Y. Hu, F. Hao, C. Meng, L Sun, D. Xu, T. Zhang, Spatial general autoregressive model-based image interpolation accommodates arbitrary scale factors, Math. Biosci. Eng., 17 (2020), 6573–6600. https://doi.org/ 10.3934/mbe.2020343 doi: 10.3934/mbe.2020343
    [8] X. Yu, Z. Chen, L. Qi, Comparative study of SARIMA and NARX models in predicting the incidence of schistosomiasis in China, Math. Biosci. Eng., 16 (2019), 2266–2276. https://doi.org/10.3934/mbe.2019112 doi: 10.3934/mbe.2019112
    [9] H. Tong, Non-Linear Time Series: A Dynamical System Approach, Oxford University Press, 1990.
    [10] D. T. Tran, A. Iosifidis, J. Kanniainen, M. Gabbouj, Temporal attention-augmented bilinear network for financial time-series data analysis, IEEE Trans. Neural Networks Learn. Syst., 30 (2018), 1407–1418. https://doi.org/10.1109/TNNLS.2018.2869225 doi: 10.1109/TNNLS.2018.2869225
    [11] D. Li, X. Wang, J. Sun, H. Yang, AI-HydSu: An advanced hybrid approach using support vector regression and particle swarm optimization for dissolved oxygen forecasting, Math. Biosci. Eng., 18 (2021), 3646–3666. https://doi.org/10.3934/mbe.2021182 doi: 10.3934/mbe.2021182
    [12] Y. C. Kuan, C. T. Hong, P. C. Chen, W. T. Liu, C. C. Chung, Logistic regression and artificial neural network-based simple predicting models for obstructive sleep apnea by age, sex, and body mass index, Math. Biosci. Eng., 19 (2022), 11409–11421. https://doi.org/10.3934/mbe.2022532 doi: 10.3934/mbe.2022532
    [13] F. Yang, D. Wang, F. Xu, Z. Huang, K. L. Tsui, Lifespan prediction of lithium-ion batteries based on various extracted features and gradient boosting regression tree model, J. Power Sources, 476 (2020), 228654. https://doi.org/10.1016/j.jpowsour.2020.228654 doi: 10.1016/j.jpowsour.2020.228654
    [14] Y. Liang, S. Zhang, H. Qiao, Y. Cheng, iEnhancer-MFGBDT: Identifying enhancers and their strength by fusing multiple features and gradient boosting decision tree, Math. Biosci. Eng., 18 (2021), 8797–8814. https://doi.org/10.3934/mbe.2021434 doi: 10.3934/mbe.2021434
    [15] H. Wan, S. Guo, K. Yin, X. Liang, Y. Lin, CTS-LSTM: LSTM-based neural networks for correlated time series prediction, Knowl. Based Syst., 191 (2020), 105239. https://doi.org/10.1016/j.knosys.2019.105239 doi: 10.1016/j.knosys.2019.105239
    [16] Y. Rizk, M. Awad, On extreme learning machines in sequential and time series prediction: A non-iterative and approximate training algorithm for recurrent neural networks, Neurocomputing, 325 (2019), 1–19. https://doi.org/10.1016/j.neucom.2018.09.012 doi: 10.1016/j.neucom.2018.09.012
    [17] Y. Liu, C. Gong, L. Yang, Y. Chen, DSTP-RNN: A dual-stage two-phase attention-based recurrent neural network for long-term and multivariate time series prediction, Expert Syst. Appl., 143 (2020), 113082. https://doi.org/10.1016/j.eswa.2019.113082 doi: 10.1016/j.eswa.2019.113082
    [18] Y. Bengio, P. Simard, P. Frasconi, Learning long-term dependencies with gradient descent is difficult, IEEE Trans. Neural Networks, 5 (1994), 157–166. https://doi.org/10.1109/72.279181 doi: 10.1109/72.279181
    [19] S. Hochreiter, J. Schmidhuber, Long short-term memory, Neural Comput., 9 (1997), 1735–1780.https://doi.org/10.1162/neco.1997.9.8.1735 doi: 10.1162/neco.1997.9.8.1735
    [20] V. Eramo, F. G. Lavacca, T. Catena, P. J. P. Salazar, Application of a long short term memory neural predictor with asymmetric loss function for the resource allocation in NFV network architectures, Comput. Networks, 193 (2021), 108104. https://doi.org/10.1016/j.comnet.2021.108104 doi: 10.1016/j.comnet.2021.108104
    [21] V. Eramo, T. Catena, Application of an innovative convolutional/LSTM neural network for computing resource allocation in NFV network architectures, IEEE Trans. Network Service Manage., 19 (2022), 2929–2943. https://doi.org/10.1109/TNSM.2022.3142182 doi: 10.1109/TNSM.2022.3142182
    [22] T. Catena, V. Eramo, M. Panella, A. Rosato, Distributed LSTM-based cloud resource allocation in network function virtualization architectures, Comput. Networks, 213 (2022), 109111. https://doi.org/10.1016/j.comnet.2022.109111 doi: 10.1016/j.comnet.2022.109111
    [23] M. Schuster, K. K. Paliwal, Bidirectional recurrent neural networks, IEEE Trans. Signal Process., 45 (1997), 2673–2681. https://doi.org/10.1109/78.650093 doi: 10.1109/78.650093
    [24] A. A. Ewees, M. A. Al-qaness, L. Abualigah, M. Abd Elaziz, HBO-LSTM: Optimized long short term memory with heap-based optimizer for wind power forecasting, Energy Convers. Manage., 268 (2022), 116022. https://doi.org/10.1016/j.enconman.2022.116022 doi: 10.1016/j.enconman.2022.116022
    [25] J. Liu, X. Huang, Q. Li, Z. Chen, G. Liu, Y. Tai, Hourly stepwise forecasting for solar irradiance using integrated hybrid models CNN-LSTM-MLP combined with error correction and VMD, Energy Convers. Manage., 280 (2023), 116804. https://doi.org/10.1016/j.enconman.2023.116804 doi: 10.1016/j.enconman.2023.116804
    [26] M. Neshat, M. M. Nezhad, N. Y. Sergiienko, S. Mirjalili, G. Piras, D. Astiaso Garcia, Wave power forecasting using an effective decomposition-based convolutional Bi-directional model with equilibrium Nelder-Mead optimizer, Energy, 256 (2022), 124623. https://doi.org/10.1016/j.energy.2022.124623 doi: 10.1016/j.energy.2022.124623
    [27] Y. Li, Z. Zhu, D. Kong, H. Han, Y. Zhao, EA-LSTM: Evolutionary attention-based LSTM for time series prediction, Knowl. Based Syst., 181 (2019), 104785. https://doi.org/10.1016/j.knosys.2019.05.028 doi: 10.1016/j.knosys.2019.05.028
    [28] G. E. Hinton, R. R. Salakhutdinov, Reducing the dimensionality of data with neural networks, Science, 313 (2006), 504–507. https://doi.org/10.1126/science.1127647 doi: 10.1126/science.1127647
    [29] X. Sun, T. Li, Q. Li, Y. Huang, Y. Li, Deep belief echo-state network and its application to time series prediction, Knowl. Based Syst., 130 (2017), 17–29. https://doi.org/10.1016/j.knosys.2017.05.022 doi: 10.1016/j.knosys.2017.05.022
    [30] X. Li, Q. Liu, Y. Wu, Prediction on blockchain virtual currency transaction under long short-term memory model and deep belief network, Appl. Soft Comput., 116 (2022), 108349. https://doi.org/10.1016/j.asoc.2021.108349 doi: 10.1016/j.asoc.2021.108349
    [31] Z. Wu, Q. Li, H. Zhang, Chain-structure echo state network with stochastic optimization: Methodology and application, IEEE Trans. Neural Networks Learn. Syst., 33 (2021), 1974–1985. https://doi.org/10.1109/TNNLS.2021.3098866 doi: 10.1109/TNNLS.2021.3098866
    [32] H. Zhang, B. Hu, X. Wang, J. Xu, L. Wang, Q. Sun, et al., Self-organizing deep belief modular echo state network for time series prediction, Knowl. Based Syst., 222 (2021), 107007. https://doi.org/10.1016/j.knosys.2021.107007 doi: 10.1016/j.knosys.2021.107007
    [33] G. E. Hinton, S. Osindero, Y. W. The, A fast learning algorithm for deep belief nets, Neural Comput., 18 (2006), 1527–1554. https://doi.org/10.1162/neco.2006.18.7.1527 doi: 10.1162/neco.2006.18.7.1527
    [34] T. Tieleman, Training restricted Boltzmann machines using approximations to the likelihood gradient, in Proceedings of the 25th International Conference on Machine Learning, 2008, 1064–1071. https://doi.org/10.1145/1390156.1390290
    [35] A. Saltelli, P. Annoni, I. Azzini, F. Campolongo, M. Ratto, S. Tarantola, Variance based sensitivity analysis of model output, Design and estimator for the total sensitivity index, Comput. Phys. Commun., 181 (2010), 259–270. https://doi.org/10.1016/j.cpc.2009.09.018 doi: 10.1016/j.cpc.2009.09.018
    [36] Y. LeCun, Y. Bengio, G. Hinton, Deep learning, Nature, 521 (2015), 436–444. https://doi.org/10.1038/nature14539 doi: 10.1038/nature14539
    [37] G. Kurnaz, A. S. Demir, Prediction of SO2 and PM10 air pollutants using a deep learning-based recurrent neural network: Case of industrial city Sakarya. Urban Climate, 41 (2022), 101051. https://doi.org/10.1016/j.uclim.2021.101051 doi: 10.1016/j.uclim.2021.101051
    [38] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, preprint, arXiv: 1412.6980, 2014. https://doi.org/10.48550/arXiv.1412.6980
    [39] H. Wang, Q. Wu, J. Xin, J. Wang, H. Zhang, Optimizing deep belief echo state network with a sensitivity analysis input scaling auto-encoder algorithm, Knowl, Based Syst., 191 (2020), 105257. https://doi.org/10.1016/j.knosys.2019.105257 doi: 10.1016/j.knosys.2019.105257
    [40] F. Zamora-Martinez, P. Romeu, P. Botella-Rocamora, J. Pardo, On-line learning of indoor temperature forecasting models towards energy efficiency, Energy Build, , 83 (2014), 162–172. https://doi.org/10.1016/j.enbuild.2014.04.034 doi: 10.1016/j.enbuild.2014.04.034
    [41] T. H. Fanaee, J. Gama, Event labeling combining ensemble detectors and background knowledge, Progress Artif. Intell., 2 (2014), 113–127. https://doi.org/10.1007/s13748-013-0040-3 doi: 10.1007/s13748-013-0040-3
    [42] J. L. Elman, Finding structure in time, Cognit. Sci., 14 (1990), 179–211. https://doi.org/10.1207/s15516709cog1402_1 doi: 10.1207/s15516709cog1402_1
    [43] J. Chung, C. Gulcehre, K. H. Cho, et al. Empirical evaluation of gated recurrent neural networks on sequence modeling, preprint, arXiv: 1412.3555.
    [44] S. Kim, M. Kang, Financial series prediction using attention LSTM, preprint, arXiv: 1902.10877.
    [45] Z. Cui, R. Ke, Z. Pu, Y. Wang, Stacked bidirectional and unidirectional LSTM recurrent neural network for forecasting network-wide traffic state with missing values, Trans. Res. Part C Emerging Technol., 118 (2020), 102674. https://doi.org/10.1016/j.trc.2020.102674 doi: 10.1016/j.trc.2020.102674
    [46] F. Karim, S. Majumdar, H. Darabi, S. Harford, Multivariate LSTM-FCNs for time series classification, Neural Networks, 116 (2019), 237–245. https://doi.org/10.1016/j.neunet.2019.04.014 doi: 10.1016/j.neunet.2019.04.014
  • This article has been cited by:

    1. Chaofeng Zhao, Zhibo Zhai, Qinghui Du, Optimal control of stochastic system with Fractional Brownian Motion, 2021, 18, 1551-0018, 5625, 10.3934/mbe.2021284
    2. Mengqing Zhang, Quanxin Zhu, Jing Tian, Convergence Rates of Partial Truncated Numerical Algorithm for Stochastic Age-Dependent Cooperative Lotka–Volterra System, 2024, 16, 2073-8994, 1659, 10.3390/sym16121659
  • Reader Comments
  • © 2023 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(3309) PDF downloads(248) Cited by(7)

Figures and Tables

Figures(21)  /  Tables(3)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog