Electronic Journal of Differential Equations, Vol. 2023 (2023), No. 02, pp. 1–18. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu or https://ejde.math.unt.edu ASYMPTOTIC STABILITY OF A STOCHASTIC AGE-STRUCTURED COOPERATIVE LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS MENGQING ZHANG, JING TIAN, KEYUE ZOU Abstract. In this article, we study a stochastic age-structured cooperative Lotka-Volterra system with Poisson jumps. Applying the M-matrix theory, we prove the existence and uniqueness of a global solution for the system. Then we use an optimized Euler-Maruyama numerical scheme to approximate the solution. We obtain second-moment boundedness and convergence rate of the numerical solutions. The numerical solutions illustrate the theoretical results. 1. Introduction The determined cooperative Lotka-Volterra system [13, 20] is described as dx(t) dt = x(t)(−α11x(t) + α12y(t) + r1), dy(t) dt = y(t)(α21x(t)− α22y(t) + r2), (1.1) where x(t), y(t) are the densities of the two cooperative species, r1, r2 present the intrinsic growth rates of the two species, α11, α22 are the intraspecific competition rates, and α12, α21 denote the interspecific cooperation rates. This system has been widely studied and has many applications [8, 24]. In particular, it can be used to describe the cooperative relationship between multiple species, such as bees and flowers, algal-fungal associations of lichens, etc. Besides the studies on the deterministic properties of cooperative Lotka-Volterra systems [4, 16, 21], there are important studies on stochastic Lotka-Volterra systems [1, 9, 17, 18]. Especially, noises with jumps are necessary to describe phenomena that burst out in nature. For example, the impact of sudden pesticide spraying or tsunami on the population can not be ignored. There are some works on the stochastic differential equations with Poisson jumps[3, 7, 22] to consider the discontinuous random effects. One example is system dx(t) = x(−α11(a)x+ α12(a)y + r1(a))dt+ σ1xdw(t) + h1xdN(t), dy(t) = y(α21(a)x− α22(a)y + r2(a))dt+ σ2ydw(t) + h2ydN(t). (1.2) 2020 Mathematics Subject Classification. 35Q30, 35B41, 76D05. Key words and phrases. Cooperative Lotka-Volterra system; age-structured; rate of convergence; Poisson jumps; Euler-Maruyama scheme. ©2023. This work is licensed under a CC BY 4.0 license. Submitted September 10, 2022. Published January 6, 2023. 1 2 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 Considering the actual biological background of system (1.2), population reproduc- tion is determined by the fertility of female individuals which is related to their age. According to [2, 5], the birth rate is sensitive to age. Therefore, we need to include an age-structured fertility rate into system (1.2). The age-structured theory was first introduced by Lotka [13] to describe the fertility of species and latter been applied in many systems [11, 19, 26, 29, 30, 31]. For example, Zhang and Liu [29] focused on the age-structured predator-prey model and discussed the non-trivial periodic oscillation phenomenon. Solis and Carrillo [19] investigated the type of theoretical predation of the predator-prey Lotka-Volterra model with age-structured. However, very few studies are on the cooperative Lotka-Volterra system. To fill this gap, we establish the stochastic age-structured cooperative Lotka-Volterra system with Poisson jumps, dtX = −∂X ∂a dt+X(−α11(a)X + α12(a)Y + r1(a))dt+ σ1Xdw(t) + h1XdN(t), dtY = −∂Y ∂a dt+ Y (α21(a)X − α22(a)Y + r2(a))dt+ σ2Y dw(t) + h2Y dN(t), X(t, 0) = ∫ A 0 γ(t, a)X(t, a)da, t ∈ [0, T ], (1.3) Y (t, 0) = ∫ A 0 β(t, a)Y (t, a)da, t ∈ [0, T ], X(0, a) = Xa(0), Y (0, a) = Ya(0), a ∈ [0, A], where X(t, a) and Y (t, a) denote the densities of the two species at time t, age a. dtX and dtY denote the differential of X(t, a) and Y (t, a) relative to t. w(t) is a standard Brownian motion, N(t) is a scalar Poisson process independent of w(t). The compensated Poisson process is Ñ(t) = N(t)− λ1t, where the parameter λ1 is the jump intensity. In the first two equations (t, a) ∈ Q and Q = (0, T ) × (0, A). Features of the parameters in system (1.3) are showed in [27]. It is extremely hard to find the explicit form of the exact solution to (1.3), so we target on finding the numerical approximation the exact solution. So we can describe the changes of number of species and predict the species size of the system. However, the Euler-Maruyama (EM) method cannot discretize the age-structured cooperative Lotka-Volterra system directly because the system does not satisfy the Lipschitz continuous condition and the linear growth condition on the drift coeffi- cients. There are many numerical approximate methods to deal with the different algebraic structure of the model, such as the positive preserving method, the trun- cated EM method, the tamed EM method and Runge-Kutta-Fehlberg method, etc. [6, 12, 14, 28]. Zhang et al. [28] constructed a positive preserving numerical method for the stochastic R&D model and proved that the convergence order is 1 2 (1 − 1 p ). Yang et al. [23] constructed the truncated EM method for the stochastic differen- tial equations to deal with the blasting phenomenon, while, the method is hard to guarantee the positive numerical solution. Liu et al. [12] established the tamed EM approximation for McKean-Vlasov stochastic differential equations with super- linear drift and Hölder diffusion coefficients. The convergence rate was obtained by replacing the one-sided local Lipschitz condition with the global condition. How- ever, most of the above mentioned methods focused on the stochastic differential equations but not involving the age-structured. In this study, we use an optimized Euler-Maruyama scheme to find a reasonable and effective numerical solution for EJDE-2023/02 LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS 3 system (1.3). Using the theory of equilibrium point, we have successfully avoided the restrictions on the coefficients. Two main contributions in this are the following: • The existence and uniqueness of the positive global solution of the stochastic age- structured cooperative Lotka-Volterra system with Poisson jumps (this is proved for the first time). • The finite-time convergence and boundedness of the numerical approximation to solution. This article is arranged as follows. In section 2, we prove the existence and uniqueness of the global solution for the stochastic age-structured cooperative Lotka-Volterra system with Poisson jumps. In section 3, we introduce the EM approximation method to the system and obtain the 2nd-moment estimations of this algorithm. The convergence of the algorithm is presented in section 4. We find a strong convergence order. The numerical simulation is reported in section 5. 2. Existence and uniqueness of a solution 2.1. Preliminaries. Throughout the paper, let (Ω,F, {Ft}t≥0,P) be a complete probability space. Here, the filtration {Ft}t≥0 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. We denote by Rn+ the positive cone in Rn, that is Rn+ = {x ∈ Rn : xi > 0 for all 1 ≤ i ≤ n}. γ̌(t, a) and γ̂(t, a) are the maximum and minimum value of the continuous function γ(t, a), respectively. For a pair of real numbers, a and b, we have a ∨ b = max{a, b} and a ∧ b = min{a, b}. For a set A, its indicator function is denoted by 1A = { 1, x ∈ A, 0, x /∈ A. For x ∈ Rn, its norms is denoted by |x|ι = { |x1|+ |x2|+ · · ·+ |xn|, ι = 1,(∑n i=1 x 2 i )1/2 , ι = 2. Let V = { ϕ ∈ L2([0, A]) : ∂ϕ ∂a ∈ L 2([0, A]), where ∂ϕ ∂a is the generalized partial derivatives with respect to age a } . H = L2([0, A]) and V ↪→ H ≡ H ′ ↪→ V ′, where V ′ is the dual space of V and H ′ is the dual space of H. The norms in V,H, V ′ are denoted as ‖‖, | · |, and ‖‖∗, respectively. The duality product be- tween 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. I2([0, T ];V ) denotes the space of all V -valued processes (Pt)t∈[0,T ], L 2 V = L2([0, T ];V ). W := (I2([0, T ];V ) ⋂ L2(Ω;C([0, T ];H)))× (I2([0, T ];V ) ⋂ L2(Ω;C([0, T ];H))). To discuss the order of convergence of the numerical methods, we have the following assumptions. (A1) limα→A− ∫ A α γ(t, a)X(t, a)da = λx > 0; limα→A− ∫ A α β(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). 4 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 (A4) ri(a) ∈ C ( [0, A];R ) and ri(a) is nondecreasing with −∞ < ri(0) < 0 < ri(A) < ∞, for i = 1, 2, and hi(t) is the continuous function on [0, T] 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. Theorem 2.1. Under Assumptions (A1)–(A5), for each initial value (Xa0(0), Ya0(0)) ∈ R2 +, there exists at most one solution (X(t), Y (t)) of system (1.3) on W . Proof. Assume Ψ1(t) := (X1(t), Y1(t)) and Ψ2(t) := (X2(t), Y2(t)) are two solutions of (1.3). Applying the Itô formula to |X1(t)−X2(t)|2, we have |X1(t)−X2(t)|2 = 2 ∫ t 0 〈−∂X1(s) ∂a + ∂X2(s) ∂a ,X1(s)−X2(s)〉ds + 2 ∫ t 0 (H1x(Ψ1(s))−H1x(Ψ2(s)), X1(s)−X2(s))ds + 2 ∫ t 0 (H2x(Ψ1(s))−H2x(Ψ2(s)), X1(s)−X2(s))ds + λ1 ∫ t 0 |J1x(Ψ1(s))− J1x(Ψ2(s))|2ds + ∫ t 0 |G1x(Ψ1(s))−G1x(Ψ2(s))|2ds + 2 ∫ t 0 (X1(s)−X2(s), (J1x(Ψ1(s))− J1x(Ψ2(s)))dN(s)) + 2 ∫ t 0 (X1(s)−X2(s), (G1x(Ψ1(s))−G1x(Ψ2(s)))dw(s)). Since 〈∂(X1(s)−X2(s)) ∂a ,X1(s)−X2(s)〉 ≤ 1 2 A2γ̌2|X1(s)−X2(s)|2, and by Assumptions (A1-A5), we have |X1(t)−X2(t)|2 ≤ (A2γ̌2 + λ1 + 2) ∫ t 0 |X1(s)−X2(s)|2ds + 2(L2 1 + 2ρ2 1 + λ1L 2 1) ∫ t 0 |Ψ1(s)−Ψ2(s)|2ds + 2 ∫ t 0 (X1(s)−X2(s), (G1x(Ψ1(s))−G1x(Ψ2(s)))dw(s)) + 2 ∫ t 0 (X1(s)−X2(s), (J1x(Ψ1(s))− J1x(Ψ2(s)))dÑ(s)). EJDE-2023/02 LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS 5 Therefore, E sup 0≤s≤t |X1(s)−X2(s)|2 ≤ (A2γ̌2 + λ1 + 2) ∫ t 0 E|X1(s)−X2(s)|2ds + 2(L2 1 + 2ρ2 1 + λ1L 2 1) ∫ t 0 E|Ψ1(s)−Ψ2(s)|2ds + 2E sup 0≤s≤t ∫ s 0 (X1(r)−X2(r), (G1x(Ψ1(r))−G1x(Ψ2(r)))dw(r)) + 2E sup 0≤s≤t ∫ s 0 (X1(r)−X2(r), (J1x(Ψ1(r))− J1x(Ψ2(r)))dÑ(r)). (2.1) Applying the Burkholder-Davis-Gundy (BDG) inequality, we have E sup 0≤s≤t ∫ s 0 (X1(r)−X2(r), (G1x(Ψ1(r))−G1x(Ψ2(r)))dw(r)) ≤ 1 8 E[ sup 0≤s≤t |X1(s)−X2(s)|2] + k1L 2 1 ∫ t 0 E|Ψ1(s)−Ψ2(s)|2ds, (2.2) and E sup 0≤s≤t ∫ s 0 (X1(r)−X2(r), (J1x(Ψ1(r))− J1x(Ψ2(r)))dÑ(r)) ≤ 1 8 E[ sup 0≤s≤t |X1(s)−X2(s)|2] + k2L 2 1 ∫ t 0 E|Ψ1(s)−Ψ2(s)|2ds, (2.3) where k1 and k2 are determined by the BDG inequality. Substituting (2.2) and (2.3) into (2.1), we have E sup 0≤s≤t |X1(s)−X2(s)|2 ≤ 1 2 E[ sup 0≤s≤t |X1(s)−X2(s)|2] + (A2γ̌2 + λ1 + 2) ∫ t 0 E sup 0≤r≤s |X1(s)−X2(s)|2ds + 2(L2 1 + 2ρ2 1 + λ1L 2 1 + k1L 2 1 + k2L 2 1) ∫ t 0 E sup 0≤r≤s |X1(s)−X2(s)|2ds. Then E sup 0≤s≤t |X1(s)−X2(s)|2 ≤ 4(L2 1 + 2ρ2 1 + λ1L 2 1 + k1L 2 1 + k2L 2 1) ∫ t 0 E sup 0≤r≤s |Ψ1(s)−Ψ2(s)|2ds + 2(A2γ̌2 + λ1 + 2) ∫ t 0 E sup 0≤r≤s |X1(s)−X2(s)|2ds. (2.4) 6 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 Similarly, we have E sup 0≤s≤t |Y1(s)− Y2(s)|2 ≤ 2(A2β̌2 + λ1 + 1 + 2ρ2 2)) ∫ t 0 E sup 0≤r≤s |Y1(r)− Y2(r)|2dr + 2(3L2 2 + 2ρ2 2 + λ1L 2 2 + 2k1L 2 2 + 2k2L 2 2) ∫ t 0 E sup 0≤r≤s |Y1(r)− Y2(r)|2ds. (2.5) By (2.4) and (2.5), we deduce that E sup 0≤s≤t |Ψ1(s)−Ψ2(s)|2 ≤ 2C ∫ t 0 E sup 0≤r≤s |Ψ1(r)−Ψ2(r)|2ds,∀t ∈ [0, T ], where C > 0 is a genetic constant whose values may vary for its different appear- ance. L1 is defined in Lemma 3.1. Now, by the Gronwall’s lemma, we obtain the uniqueness of the solution. � Theorem 2.2. Under Assumptions (A1)–(A5), for each given initial value (Xa0 (0), Ya0 (0)) ∈ R2 +, there exists a nonnegative global solution (X(t), Y (t)) of system (1.3) on W . Proof. It is easy to see that for any given initial value (Xa0 (0), Ya0 (0)) ∈ R2 +, there exists a local solution (X(t, a), Y (t, a)) ∈ W when t ∈ [0, τe), where τe is the explosion time of system (1.3). Next, we need to show that τe =∞, a.s. We choose a positive constant k0 > 1 such that for any initial value (Xa0 (0), Ya0 (0)) ∈ [ 1 k0 , k0]× [ 1 k0 , k0]. Then, for each k ≥ k0, the stopping time τk is defined by τk = inf{t ∧ a ∈ [0, τe) : (X,Y ) /∈ ( 1 k , k)× ( 1 k , k)}. Next, we should prove that τk is an increasing function and limk→∞ τk = τ∞, where τ∞ ≤ τ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 satisfying P{τk ≤ T} ≥ ε,∀k ≥ k1. (2.6) Now, for any constants cx, cy > 0, we define a function V : W → R+ by V (X,Y ) = cx(X − logX − 1) + cy(Y − log Y − 1). Using the Itô formula [15], we have LV (X,Y ) = cx 〈 1− 1 X ,−∂X ∂a +X(−α11X + α12Y + r1) 〉 + 1 2 (cxσ 2 1 + cyσ 2 2) + cy 〈 1− 1 Y ,−∂Y ∂a + Y (α21X − α22Y + r2) 〉 + 1 2 λ1(cxh 2 1 + cyh 2 2) ≤ Ψ> diag(cx, cy)R̄− C̄(ĀΨ + R̄) + 1 2 Ψ>(diag(cx, cy)Ā+ Ā> diag(cx, cy))Ψ + cx ∫ A 0 1 X daX + cy ∫ A 0 1 Y daY + 1 2 (cxσ 2 1 + cyσ 2 2) + 1 2 λ1(cxh 2 1 + cyh 2 2) EJDE-2023/02 LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS 7 ≤ cx log λx + cy log λy − [ cx log | ∫ A 0 γ(t, a)X(t, a)da| + cy log | ∫ A 0 β(t, a)Y (t, a)da| ] + 1 2 (cxσ 2 1 + cyσ 2 2) + Ψ> diag(cx, cy)R̄ − C̄(ĀΨ + R̄) + 1 2 λ1(cxh 2 1 + cyh 2 2) + 1 2 Ψ>(diag(cx, cy)Ā+ Ā> diag(cx, cy))Ψ, where Ψ = ( X(t, a) Y (t, a) ) , R̄ = ( r1(a) r2(a) ) , C̄> = ( cx cy ) , Ā = ( −α11(a) α12(a) α21(a) −α22(a) ) , diag(cx, cy) = ( cx 0 0 cy ) . By Assumptions (A1)–(A5), since −Ā is a nonsingular M-matrix, we have LV (X,Y ) ≤ cx log λx + cy log λy − [ cx log ∣∣ ∫ A 0 γ(t, a)X(t, a)da ∣∣ + cy log ∣∣ ∫ A 0 β(t, a)Y (t, a)da ∣∣]+ 1 2 (cxσ 2 1 + cyσ 2 2) + Ψ> diag(cx, cy)R̄− C̄(ĀΨ + R̄) + 1 2 λ1(cxh 2 1 + cyh 2 2) ≤ cx log λx + cy log λy − [ cx log ∣∣ ∫ A 0 γ(t, a)X(t, a)da ∣∣ + cy log ∣∣ ∫ A 0 β(t, a)Y (t, a)da ∣∣]+ {R̄> diag(cx, cy)− C̄Ā}Ψ + 1 2 (cxσ 2 1 + cyσ 2 2)− r1cx − r2cy + 1 2 λ1(cxh 2 1 + cyh 2 2), and LV (x, y) ≤ ∣∣R̄> diag(cx, cy)− C̄Ā ∣∣ · |Ψ|+ C ≤M1(1 + |Ψ|). (2.7) Letting M2 = cx∨cy cx∧cy and M3 = cx ∧ cy, we obtain M2V (X,Y ) = cx ∨ cy cx ∧ cy [ cx(X − 1− logX) + cy(Y − 1− log Y ) ] ≥ V (X,Y ). Since |Ψ| = (X2 + Y 2)1/2 ≤ X + Y , it yields |Ψ| ≤ 2(X − 1− logX) + 2(Y − 1− log Y ) + 4 ≤ 4 + 2 M3 V (X,Y ). (2.8) By (2.7) and (2.8), we obtain LV (X,Y ) ≤M4(1 + V (X,Y )), where M4 = M1(5 ∨ 2/M3). Then, we have EV (X(t ∧ τk), Y (t ∧ τk)) ≤ V (X0, Y0) + E ∫ t∧τk 0 M4(1 + V (X(s), Y (s)))ds ≤M5 +M4 ∫ t 0 EV (X(s ∧ τk), Y (s ∧ τk))ds, where M5 = V (X0, Y0) +M4T . 8 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 By the Gronwall inequality, we have EV (X(τk ∧ T ), Y (τk ∧ T )) ≤M5e M4T . According to (2.6) and for any k ≥ k1, ω ∈ Ωk, (t, a) ∈ Q, X(τk, ω) and Y (τk, ω) equal either k or 1 k , hence we have V (X(τk, ω), Y (τk, ω)) ≥ [ cx(k − 1− log k) + cy(k − 1− log k) ] ∧ [ cx( 1 k − 1 + log k) + cy( 1 k − 1 + log k) ] , and M5e M4T ≥ ε [ cx(k−1−log k)+cy(k−1−log k) ] ∧ [ cx( 1 k −1+log k)+cy( 1 k −1+log k) ] . Taking k →∞, we obtain a contradictiion, ∞ > M5e M4T ≥ ∞. � 3. Euler Maruyama (EM) method In this section, we apply the EM approximate solution to system (1.3), and discuss the boundedness moments of the scheme. First, we introduce notation. H1x(Ψ) := H1x(X,Y ) := r1(a)X(t, a), H1y(Ψ) := H1y(X,Y ) := r2(a)Y (t, a), G1x(Ψ) := G1x(X,Y ) := σ1X(t, a) := σ1X(t), G1y(Ψ) := G1y(X,Y ) := σ2Y (t, a) := σ2Y (t), H2x(Ψ) := H2x(X,Y ) := X(t, a)[−α11(a)X(t, a) + α12(a)Y (t, a)], H2y(Ψ) := H2y(X,Y ) := Y (t, a)[α21(a)X(t, a)− α22(a)Y (t, a)], J1x(Ψ) := J1x(X,Y ) := h1(t)X(t, a) := h1(t)X(t), J1y(Ψ) := J1y(X,Y ) := h2(t)Y (t, a) := h2(t)Y (t). The equilibria of (1.3) are: 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) . Lemma 3.1. Under Assumptions (A1)–(A5), for each ψ1, ψ2 ∈W , we have |H2x(ψ1)−H2x(ψ2)| ≤ ρ1|ψ1 − ψ2|1 ≤ 21/2ρ1|ψ1 − ψ2|2, |H2y(ψ1)−H2y(ψ2)| ≤ ρ2|ψ1 − ψ2|1 ≤ 21/2ρ2|ψ1 − ψ2|2, |H1x(ψ1)−H1x(ψ2)| ∨ |G1x(ψ1)−G1x(ψ2)| ∨ |J1x(ψ1)− J1x(ψ2)| ≤ L1|ψ1 − ψ2|, |H1y(ψ1)−H1y(ψ2)| ∨ |G1y(ψ1)−G1y(ψ2)| ∨ |J1y(ψ1)− J1y(ψ2)| ≤ L2|ψ1 − ψ2|, where L1, L2 are constants, ρ1 = 2x∗ − r1(0) + α12(x∗ + y∗), ρ2 = 2y∗ − r2(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. [25], we omit it here. Corollary 3.2. From Lemma 3.1, we can conclude that there exist constants K1 and K2 such that |H1x(ψ)| ∨ |G1x(ψ)| ∨ |J1x(ψ)| ≤ K1(1 + |ψ|), |H1y(ψ)| ∨ |G1y(ψ)| ∨ |J1y(ψ)| ≤ K2(1 + |ψ|), EJDE-2023/02 LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS 9 and x ·H2x(ψ) ≤ 21/2ρ1(1 + |ψ|2), y ·H2y(ψ) ≤ 21/2ρ2(1 + |ψ|2), where ψ(t, a) is the numerical solution of (1.3) and ψ(t, a) = (x(t, a), y(t, a))>. From now on, we fix T > 0 and the time step ∆ ∈ (0, 1). Let tk = k∆ for k = 0, 1, 2, . . . , [T/∆], where [T/∆] denotes the integer part of T/∆. Then, we form the discrete time EM solutions by setting X(0) = Xa(0), Y (0) = Ya(0) and computing xk+1 = xk + { − ∂xk ∂a + r1(a)xk − α11(a)xkxk + α12(a)xkyk } ∆ + σ1xk∆wk + h1xk∆Nk, yk+1 = yk + { − ∂yk ∂a + r2(a)yk + α21(a)xkyk − α22(a)ykyk } ∆ + σ2yk∆wk + h2yk∆Nk, where ∆wk = wtk+1 − wtk , ∆Nk = Ntk+1 −Ntk . The continuous time process is defined as x(t) = x0 − ∫ t 0 ∂x̄(s) ∂a ds+ ∫ t 0 {r1(a)x̄(s)− α11(a)x̄(s)x̄(s) + α12(a)x̄(s)ȳ(s)}ds + ∫ t 0 σ1x̄(s)dw(s) + ∫ t 0 h1x̄(s)dN(s), y(t) = y0 − ∫ t 0 ∂ȳ(s) ∂a ds+ ∫ t 0 {r2(a)ȳ(s) + α21(a)x̄(s)ȳ(s)− α22(a)ȳ(s)ȳ(s)}ds + ∫ t 0 σ2ȳ(s)dw(s) + ∫ t 0 h2ȳ(s)dN(s). (3.1) where x̄(t) = [ T ∆ ]∑ k=0 xk1[tk,tk+1)(t), ȳ(t) = [ T ∆ ]∑ k=0 yk1[tk,tk+1)(t) are the step processes. We can easily see that x̄(t) = xk and ȳ(t) = yk for t ∈ [tk, tk+1) when k = 0, 1, 2, . . . , [ T∆ ]. In the following proof, we use ψ(t) = (x(t), y(t)) to represent the numerical solution of the system (1.3). Theorem 3.3. Under Assumptions (A1)–(A5), there exists a constant C > 0 such that sup 0≤t≤T E|ψ(t)|2 ≤ C, where C dependents only on T . Proof. Applying the Itô formula to the first equation of (3.1), we have |x(t)|2 = |x0|2 + 2 ∫ t 0 〈−∂x(s) ∂a +H1x +H2x, x(s)〉ds+ ∫ t 0 ||G1x||2ds + 2 ∫ t 0 (x(s), G1xdw(s)) + 2 ∫ t 0 (x(s), J1xdN(s)) + λ1 ∫ t 0 |J1x|2ds. Since 〈 − ∂x(s) ∂a , x(s) 〉 = − ∫ A 0 x(s)dax(s) = 1 2 [ ∫ A 0 γ(s, a)x(s, a)da ]2 10 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 ≤ 1 2 ∫ A 0 γ2(s, a)da ∫ A 0 x2(s, a)da ≤ 1 2 A2γ̌2(s, a)|x(s, a)|2, by the assumptions we have |x(t)|2 ≤ |x0|2 +A2γ̌2(s, a) ∫ t 0 |x(s)|2ds+ 2 ∫ t 0 (H1x +H2x, x(s))ds+ 2 ∫ t 0 (x(s), G1xdw(s)) + ∫ t 0 ||G1x||2ds+ 2 ∫ t 0 (J1xdN(s), x(s)) + λ1 ∫ t 0 |J1x|2ds ≤ |x0|2 +A2γ̌2(s, a) ∫ t 0 |x(s)|2ds+ 2 ∫ t 0 |H1x|2ds+ 2 ∫ t 0 |H2x|2ds+ ∫ t 0 |x(s)|2ds + 2 ∫ t 0 (x(s), G1xdw(s)) + ∫ t 0 ||G1x||2ds+ 2 ∫ t 0 (J1xdN(s), x(s)) + λ1 ∫ t 0 |J1x|2ds ≤ |x0|2 + (A2 · γ̌2 + 1) ∫ t 0 |x(s)|2ds+ [(3 + λ1)L2 1 + 4ρ2 1] ∫ t 0 |ψ(s)|2ds + 2 ∫ t 0 (x(s), G1xdw(s)) + 2 ∫ t 0 (x(s), J1xdÑ(s)). Then E sup 0≤s≤t |x(s)|2 ≤ (A2 · γ̌2 + 1) ∫ t 0 E sup 0≤s≤t |x(s)|2ds + [(3 + λ1)L2 1 + 4ρ2 1] ∫ t 0 E sup 0≤s≤t |ψ(s)|2ds + E|x0|2 + 2E sup 0≤s≤t ∫ s 0 (x(r), G1xdw(r)) + 2E sup 0≤s≤t ∫ s 0 (x(r), J1xdÑ(r)). (3.2) By the BDG inequality, we obtain E sup 0≤s≤t ∫ s 0 (x(r), G1xdw(r)) ≤ 1 8 E[ sup 0≤s≤t |x(s)|2] + C1L 2 1 ∫ t 0 E|ψ(s)|2ds, (3.3) E sup 0≤s≤t ∫ s 0 (x(r), J1xdÑ(r)) ≤ 1 8 E[ sup 0≤s≤t |x(s)|2] + C2 ∫ t 0 E|ψ(s)|2ds, (3.4) where C1 and C2 are determined by the BDG inequality. According to (3.2), (3.3), and (3.4), we have E sup 0≤s≤t |x(s)|2 ≤ [(3 + λ1 + 2C1 + 2C2)L2 1 + 4ρ2 1] ∫ t 0 E sup 0≤s≤t |ψ(s)|2ds + 2E|x0|2 + (A2 · γ̌2 + 1) ∫ t 0 E sup 0≤r≤s |x(r)|2ds. (3.5) EJDE-2023/02 LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS 11 Similarly, E sup 0≤s≤t |y(s)|2 ≤ 2E|y0|2 + (A2β̌2 + 2λ1 + 1) ∫ t 0 E sup 0≤s≤t |y(s)|2ds + [(3 + 2λ1 + 2C1 + 2C2)L2 2 + 4ρ2 2] ∫ t 0 E sup 0≤s≤t |ψ(s)|2ds. (3.6) By (3.5) and (3.6), we have E sup 0≤s≤t |ψ(t)|2 ≤ 2E|ψ0|2 + [2A2(β̌2 ∨ γ̌2) + 2 + 4λ1] ∫ t 0 E sup 0≤r≤s |ψ(r)|2ds + [(6 + 2λ1 + 2C1 + 2C2)(L2 1 ∨ L2 2) + 2(ρ2 1 ∨ ρ2 2)] ∫ t 0 E sup 0≤s≤t |ψ(s)|2ds. By Gronwall’s inequality, E supt∈[0,T ] |ψ(t)|2 ≤ C. � Based on Theorem 3.3, we can prove the convergence of ψ(t) and ψ̄(t). Theorem 3.4. Under Assumptions (A1)–(A5), we have E sup 0≤t≤T |ψ(t)− ψ̄(t)|2 ≤ C∆, where C only dependents on T . Proof. First, we have x(t)− x̄(t) = − ∫ t tk ∂x(s) ∂a ds+ ∫ t tk [H1x +H2x]ds+ ∫ t tk G1xdw(s) + ∫ t tk J1xdN(s), thus |x(t)− x̄(t)|2 ≤ 4| ∫ t tk ∂x(s) ∂a ds|2 + 4| ∫ t tk [H1x +H2x]ds|2 + 4| ∫ t tk G1xdw(s)|2 + 4| ∫ t tk J1xdN(s)|2 ≤ 4∆ ∫ t tk |∂x(s) ∂a |2ds+ 4∆ ∫ t tk |H1x +H2x|2ds+ 4| ∫ t tk G1xdw(s)|2 + 8|λ1 ∫ t tk J1xds|2 ≤ 4∆ ∫ t tk |∂x(s) ∂a |2ds+ 8∆[ ∫ t tk (|H2 1x|+ |H2 2x|)ds] + 4| ∫ t tk G1xdw(s)|2 + 8| ∫ t tk J1xdN̄(s)|2 + 8|λ1 ∫ t tk J1xds|2. by Lemma 3.1, E sup t∈[0,T ] |x(t)− x̄(t)|2 ≤ 5E sup t∈[0,T ] max k=0,1,...,N−1 ∣∣ ∫ t tk G1x(ψ)dw(s) ∣∣2 + 8E sup t∈[0,T ] max k=0,1,...,N−1 ∣∣ ∫ t tk J1x(ψ)dÑ(s) ∣∣2 + 5∆ ∫ t tk |∂x(s) ∂a |2ds+ 8TC[L2 1 + 2ρ1 + 8λ2 1L 2 1]∆. 12 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 According to the Doob inequality, we obtain E sup t∈[0,T ] |x(t)− x̄(t)|2 ≤ 5∆ ∫ t tk |∂x(s) ∂a |2ds+ 8TC[L2 1 + 2ρ1 + 8λ2 1L 2 1]∆ + 5 max k=0,1,...,N−1 ∫ (k+1)∆ k∆ E|G1x(x, y)|2ds + 8λ1 max k=0,1,...,N−1 ∫ (k+1)∆ k∆ E|J1x(x, y)|2ds ≤ 5∆ ∫ t tk |∂x(s) ∂a |2ds+ 8TC[L2 1 + 2ρ1 + 8λ2 1L 2 1]∆ + 5L2 1C∆ + 8λ1L 2 1C∆. (3.7) Since y(t)− ȳ(t) = ∫ t tk ∂y(s) ∂a ds+ ∫ t tk [H1y +H2y]ds+ ∫ t tk G1ydw(s) + ∫ t tk J1ydN(s), we have |y(t)− ȳ(t)|2 ≤ 4| ∫ t tk ∂y(s) ∂a ds|2 + 4| ∫ t tk [H1y +H2y]ds|2 + 4| ∫ t tk G1ydw(s)|2 + 4| ∫ t tk J1ydN(s)|2 ≤ 4∆ ∫ t tk |∂y(s) ∂a |2ds+ 4∆ ∫ t tk |H1y +H2y|2ds+ 4| ∫ t tk G1ydw(s)|2 + 8| ∫ t tk J1ydN̄(s)|2 + 8|λ1 ∫ t tk J1yds|2 ≤ 4∆ ∫ t tk |∂y(s) ∂a |2ds+ 8∆ ∂y(s) ∂a [|H2 1y + |H2y|2]ds+ 4| ∫ t tk G1ydw(s)|2 + 8| ∫ t tk J1ydN̄(s)|2 + 8|λ1J1yds|2. Using Lemma 3.1, we obtain E sup t∈[0,T ] |y(t)− ȳ(t)|2 ≤ 4∆ ∫ t tk |∂y(s) ∂a |2ds+ 8TC∆[L2 2 + 2ρ2] + 8λ2 1L 2 2TC∆ + 4E sup t∈[0,T ] max 0≤k≤N−1 ∣∣ ∫ t k∆ G1y(ψ)dw(s) ∣∣2 + 4E sup t∈[0,T ] ∣∣ ∫ t k∆ J1y(ψ)dÑ(s) ∣∣2 ≤ 4∆ ∫ t tk |∂y(s) ∂a |2ds+ C[8T (L2 2 + 2ρ2 + λ2 1L 2 2) + 4L2 2 + 8λ1L 2 2]∆. (3.8) By (3.7) and (3.8), we deduce that E sup t∈[0,T ] |ψ(t)− ψ̄(t)|2 ≤ 4∆ ∫ t tk |∂ψ(s) ∂a |2ds+ C ′∆. EJDE-2023/02 LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS 13 Then, we obtain the result. � Theorem 3.4 shows that the numerical solution converges to the step process. Now, we can discuss the convergence relation between the true and numerical so- lution of the system (1.3). 4. Convergence rates over the time interval [0, T ] Theorem 4.1. Under Assumptions (A1)–(A5), there exists a constant C such that E|Ψ(t)− ψ(t)|2 < C∆. Proof. For any t ∈ [0, T ], there exists a positive integer k such that t ∈ [tk, tk+1], and X(t)− x(t) = − ∫ t 0 ∂X(s)− x(s) ∂a ds + ∫ t 0 [H1x(Ψ(s)) +H2x(Ψ(s))−H1x(ψ̄(s))−H2x(ψ̄(s))]ds + ∫ t 0 (G1x(Ψ(s))−G1x(ψ(s)))dw(s) + ∫ t 0 (J1x(Ψ(s))− J1x(ψ̄(s)))dN(s). According to the Itô formula, we have d|X(t)− x(t)|2 = −2〈X(t)− x(t), ∂(X(t)− x(t)) ∂a 〉dt+ 2(X(t)− x(t), H1x(Ψ(s)−H1x(ψ̄(s))))dt + 2(X(t)− x(t), H2x(Ψ(s)−H2x(ψ̄(s))))dt + 2(X(t)− x(t), (G1x(Ψ(t))−G1x(ψ̄(t)))dw(t)) + |G1x(Ψ(t))−G1x(ψ̄(t))|2dt+ λ1|J1x(Ψ(t))− J1x(ψ̄(t))|2dt + 2(X(t)− x(t), (J1x(Ψ(t))− J1x( ¯ψ(t)))dN(t)) ≤ A2γ̌2|X(t)− x(t)|2dt+ L1|X(t)− x(t)|2dt+ L1|Ψ(t)− ψ̄(t)|2dt + |X(t)− x(t)|2dt+ 2ρ2 1|Ψ(t)− ψ̄(t)|2dt + 2(X(t)− x(t), (G1x(Ψ(t))−G1x(ψ̄(t)))dw(t)) + L2 1|Ψ(t)− ψ(t)|2dt+ λ1L 2 1|Ψ(t)− ψ̄(t)|2dt + 2(X(t)− x(t), (J1x(Ψ(t))− J1x(ψ̄(t)))dÑ(t)) + λ1|X(t)− x(t)|2 + λ1|J1x(Ψ(t))− J1x(ψ̄(t))|2dt. 14 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 Moreover, E sup s∈[0,t] |X(t)− x(t)|2 ≤ ((A2)γ̌2 + L1 + λ1 + 1) ∫ t 0 E sup s∈[0,t] |X(r)− x(r)|2ds + ((2λ1 + 1)L2 1 + 2ρ2 1 + 1)E ∫ t 0 |Ψ(s)− ψ̄(s)|2ds + E sup s∈[0,t] ∫ t 0 2(X(s)− x(s), (G1x(Ψ(s))−G1x(ψ̄(s)))dw(s)) + 2E sup s∈[0,t] ∫ t 0 (X(s)− x(s), (J1x(Ψ(s))− J1x(ψ̄(s)))dÑ(s)). (4.1) By the BDG inequality, we have E sup s∈[0,t] ∫ s 0 (X(s)− x(s), (G1x(Ψ(s))−G1x(ψ̄(s)))dw(s)) ≤ 1 8 E[ sup 0≤s≤t |X(s)− x(s)|2] + k1 ∫ t 0 E|Ψ(s)− ψ̄(s)|2ds, (4.2) and E sup s∈[0,t] ∫ s 0 (X(s)− x(s), (J1x(Ψ(s))− J1x(ψ̄(s)))dÑ(s)) ≤ 1 8 E[ sup 0≤s≤t |X(s)− x(s)|2] + k2 ∫ t 0 E|Ψ(s)− ψ̄(s)|2ds, (4.3) where k1 and k2 are two positive constants. Substituting (4.2) and (4.3) into (4.1), we have E sup s∈[0,t] |X(s)− x(s)|2 ≤ (A2γ̌2 + L1 + λ1 + 1) ∫ t 0 E sup r∈[0,s] |X(s)− x(s)|2ds + 2(L2 1 + 2λ1L 2 1 + 2ρ2 1 + L1) ∫ t 0 E sup r∈[0,s] |Ψ(s)− ψ(s)|2ds + 2(L2 1 + 2λ1L 2 1 + 2ρ2 1 + L1)E sup r∈[0,s] |ψ(s)− ψ̄(s)|2ds + 1 4 E[ sup 0≤s≤t |X(s)− x(s)|2] + 2k1 ∫ t 0 E|Ψ(s)− ψ̄(s)|2ds + 1 4 E[ sup 0≤s≤t |X(s)− x(s)|2] + 2k2 ∫ t 0 E|Ψ(s)− ψ̄(s)|2ds = 1 2 E[ sup 0≤s≤t |X(s)− x(s)|2] + (A2γ̌2 + L1 + λ1 + 1) ∫ t 0 E sup r∈[0,s] |X(s)− x(s)|2ds + 2(L2 1 + 2λ1L 2 1 + 2ρ2 1 + L1 + 2k1 + 2k2) ∫ t 0 E sup r∈[0,s] |Ψ(s)− ψ(s)|2ds + 2(L2 1 + 2λ1L 2 1 + 2ρ2 1 + L1 + 2k2) ∫ t 0 E sup s∈[0,t] |ψ(s)− ψ̄(s)|2ds. We apply the same method to |Y (t)− y(t)|2 and complete the proof. � EJDE-2023/02 LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS 15 Corollary 4.2. Under Assumptions (A1)–(A5), the numerical solution of system (1.3) will converge to the true solution in the sense that lim ∆t→0 E ( sup 0≤t≤T |Ψ(t)− ψ(t)|2 ) = 0. From Theorem 4.1 and Corollary 4.2, we can conclude that the numerical solution ψ(t) and the true solution Ψ(t) are close to each other, which shows that the numerical algorithm constructed in this paper is effective. 5. Numerical experiments To illustrate the theorems in this paper, we consider the stochastic age-structured cooperative LV system for numerical simulations, dtX = −∂X ∂a dt+X[− 1 (1− a)2 X + cos2 aY − 2]dt− 1 2 dw(t) + 1 3 x(t−)dN(t), dtY = −∂Y ∂a dt+ Y [sin2 aX − e 1 aY + 3a− 2]dt− 1 2 dw(t) + 1 3 y(t−)dN(t), X(t, 0) = ∫ 1 0 1 1− a X(t, a)da, t ∈ [0, 1], (5.1) Y (t, 0) = ∫ 1 0 1 1− a Y (t, a)da, t ∈ [0, 1], X(0, a) = e −2 1−a , Y (0, a) = e −2 1−a , a ∈ [0, 1], where α11(a) = 1 (1−a)2 , α12(a) = cos2 a, α21(a) = sin2 a, α22(a) = e− 1 a , and γ(t, a) = β(t, a) = 1 1−a . In the first two equations, we have (t, a) ∈ Q We simulate the EM numerical approximate solution of the two species respectively (see Figure 1). Approximate solution X(t, a) Approximate solution Y (t, a) Figure 1. EM numerical approximate solutions of system (5.1). Next, we simulate the true solution of the system (5.1) by using the method in literature [10]. The Figure 2 shows that the true solutions of X(t, a) and Y (t, a) with perturbation, respectively. Figure 3 and Figure 4 are the error between the true solution and the numerical approximate solution of system (5.1). We can observe that the maximum value of 16 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 True solution X(t, a) True solution Y (t, a) Figure 2. True solutions with the perturbation of system (5.1). the squared error of X(t, a) and Y (t, a) are below 0.3 and 0.15, respectively, which means that the numerical approximate solution convergent to the true solution effectively. Error beqtween X(t, a) and x(t, a) Error square beqtween X(t, a) and x(t, a) Figure 3. Simulation error in species X(t, a) of system (5.1). Acknowledgements This research was supported by the Research Start-Up Project of North Minzu University (2021KYQD33), by the Basic Research Project of North Minzu Univer- sity (2021XYSYK01). References [1] Y. T. Cai, C. C. Wang, D. Fan; Stability and bifurcation in a delayed predator-prey model with Holling-type IV response function and age structure, Electron. J. Differ. Equ., vol. 2022 no. 42 (2021) pp. 1-16. [2] M. Delgado, A. Suárez; Age-dependent diffusive Lotka-Volterra type systems, Math. Comput. Model., vol. 45 (2007) pp. 668-680. [3] S.N. Deng, W. Y. Fei, W. Liu, X. R. Mao; The truncated EM method for stochastic differential equations with Poisson jumps, J. Comput. Appl. Math., vol. 355 (2019) pp. 232-257. EJDE-2023/02 LOTKA-VOLTERRA SYSTEM WITH POISSON JUMPS 17 Error between Y (t, a) and y(t, a) Error square between Y (t, a) and y(t, a) Figure 4. Simulation error in species Y (t, a) of system (5.1). [4] J. H. Huang, X. F. Zou; Traveling wavefronts in diffusive and cooperative Lotka-Volterra system with delays, J. Math. Anal. Appl., vol. 271 (2002) pp. 455-466. [5] N. Hritonenko, Y. Yatsenko; The structure of optimal time- and age-dependent harvesting in the Lotka-McKendrik population model, Math. Biosci., vol. 208 (2007) pp. 48-62. [6] R. A. Jahdali, L. Dalcin, R. Boukharfane, I. R. Nolasco, D. E. Keyes, M. Parsani; Optimized explicit Runge-Kutta schemes for high-order collocated discontinuous Galerkin methods for compressible fluid dynamics. Comput. Math. Appl., 118 (2022) pp. 1-17. [7] H. D. Li, Q. X. Zhu; The pth moment exponential stability and almost surely exponential stability of stochastic differential delay equations with Poisson jump, J. Math. Anal. Appl., vol. 471 (2019) pp. 197-210. [8] S. B. Li, Y. Xiao, Y. Dong; Diffusive predator-prey models with fear effect in spatially het- erogeneous environment, Electron. J. Differ. Equ., 2021 (2021), no. 70, pp. 1-31. [9] S. Li, S. Guo; Dynamics of stochastic Lotka-Volterra predator-prey models driven by three independent Brownian motions, Electron. J. Differ. Equ., 2022 no. 32 (2022) pp. 1-28. [10] X. Y. Li; Variational iteration method for nonlinear age-structured population models, Com- put. Math. Appl., vol. 58 (2009) pp. 2177-2181. [11] Y. Li, M. Ye, Q. M. Zhang; Strong convergence of the partially truncated Euler-Maruyama scheme for a stochastic age-structured SIR epidemic model, Appl. Math. Comput., vol. 362 (2019) pp. 1-22. [12] H. G. Liu, B. B. Shi, F. K. Wu; Tamed Euler-Maruyama approximation of McKean-Vlasov stochastic differential equations with super-linear drift and Hölder diffusion coefficients. Appl. Numer. Math., 183 (2023) pp. 56-85. [13] A. J. Lotka; Relation between birth rates and death rates, Science, vol. 26 (1907) pp. 121-130. [14] X. R. Mao, F. Y. Wei, T. Wiriyakraikul; Positivity preserving truncated Euler-Maruyama Method for stochastic Lotka-Volterra competition model, J. Comput. Appl. Math., 394 (2021) pp. 113566. [15] X. R. Mao; Stochastic Differential Equations and Applications, second ed., Horwood, UK, 2007. [16] Y. Nakata, Y. Muroya; Permanence for nonautonomous Lotka-Volterra cooperative systems with delays, Nonlinear. Anal-Real., vol. 11 (2010) pp. 528-534. [17] N. N. Nguyen, G. Yin; Stochastic Lotka-Volterra competitive reaction-diffusion systems per- turbed by space-time white noise: Modeling and analysis, J. Differ. Equations, vol. 282 (2021) pp. 184-232. [18] R. Rudnicki; Long-time behaviour of a stochastic prey-predator model, Stoch. Process. their Appl., vol. 108 (2003), 93?107. [19] F. J. Solis, R. A. Ku-Carrillo; Generic predation in age structure predator-prey models, Appl. Math. Comput., vol. 231 (2014) pp. 205-213. [20] V. Volterra; Variazioni e fluttuazioni del numero d’individui in specie animali conviventi. Societa anonima tipografica Leonardo da Vinci, 1926. 18 M. ZHANG, J. TIAN, K. ZOU EJDE-2023/02 [21] L. Xu, J. Y. Liu, G. Zhang; Pattern formation and parameter inversion for a discrete Lotka- Volterra cooperative system, Chaos. Soliton. Fract., vol. 110 (2018) pp. 226-231. [22] Z. G. Yan, M. Zhang, G. Z. Chang, H. Lv, J. H. Park; Finite-time annular domain stability and stabilization of Itô stochastic systems with Wiener noise and Poisson jumps-differential Gronwall inequality approach, Appl. Math. Comput., vol. 412 (2022) pp. 126589. [23] H. Yang, F. Wu, P. Kloeden, X. Mao; The truncated Euler-Maruyama method for stochastic differential equations with Hölder diffusion coefficients, J. Comput. Appl. Math. vol. 366 (2020) pp. 112379. [24] K. Yang, H. L. Smith; Convergence in Lotka?Volterra-type delay systems without instanta- neous feedbacks, Proc. R. Soc. Edinb. A: Math., 123 (1993) pp. 45-58. [25] Y. Yang, C. F. Wu, Z. X. Li; Forced waves and their asymptotics in a Lotka-Volterra coop- erative model under climate change, Appl. Math. Comput. vol. 353 (2019) pp. 254-264. [26] Q. M. Zhang, W. A. Liu, Z. K. Nie; Existence, uniqueness and exponential stability for stochastic age-dependent population, Appl. Math. Comput., vol. 154 (2004) pp. 183-201. [27] M. Q. Zhang, Q. M. Zhang, J. Tian, X. N. Li; The asymtotic stability of numerical analysis for stochastic age-dependent cooperative Lotka-Volterra system, Math. Biosci. Eng., vol. 18 (2021) pp. 1425-1449. [28] M. Q. Zhang, Q. M. Zhang; A positivity preserving numerical method for stochastic R&D model, Appl. Math. Comput., vol. 351 (2019) pp. 193-203. [29] X. M. Zhang, Z. H. Liu; Periodic oscillations in age-structured ratio-dependent predator-prey model with Michaelis-Menten type functional response, Physica D. vol. 389 (2019) pp. 51-63. [30] R. H. Li, P. K. Leung, W. K. Pang; Convergence of numerical solutions to stochastic age- dependent population equations with Markovian switching, J. Comput. Appl. Math. vol. 233 (2009) pp. 1046-1055. [31] Y. Zhao, S. L. Yuan, Q.M. Zhang; The effect of Levy noise on the survival of a stochastic competitive model in an impulsive polluted environment, Appl. Math. Model., vol. 40 (2016) pp. 7583-7600. Mengqing Zhang School of Mathematics and Information Science, North Minzu University, Yinchuan 750021, China Email address: mengqingzhang@126.com Jing Tian (corresponding author) Department of Mathematics, Towson University, Towson, MD 21252, USA Email address: jtian@towson.edu Keyue Zou School of Finance, Southwestern University of Finance and Economics, Chengdu 611130, China Email address: zoukeyue@outlook.com 1. Introduction 2. Existence and uniqueness of a solution 2.1. Preliminaries 3. Euler Maruyama (EM) method 4. Convergence rates over the time interval [0, T] 5. Numerical experiments Acknowledgements References