The asymptotic stability of numerical analysis for stochastic age-dependent cooperative Lotka-Volterra system

1.   Introduction

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 $R_{0}$ 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 $p$th-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 $p$th-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.

2.   Preliminaries

2.1.   Model formulation

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

where $x(t), \ y(t)$ represent the density of the two cooperative species, $r_1, r_2$ are the intrinsic growth rates of the two species, $\alpha_{11}, \ \alpha_{22}$ denote the intraspecific competition rates, $\alpha_{12}, \alpha_{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:

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

Now, we consider the effect of random environment noise in Eq (2.2) and assume $r_{i}(a)\rightarrow r_{i}(a)+\sigma_i \dot{\omega}(t), \ (i = 1, 2)$, where $w(t)$ is a standard Brownian motion. So, we have the stochastic age-dependent cooperative LV system.

where $d_tX$ and $d_tY$ are the differential of $X(t, a)$ and $Y(t, a)$ relative to $t$, i.e., $d_tX = \frac{\partial X}{\partial t}dt,$ and $d_tY = \frac{\partial Y}{\partial t}dt$. Features of the parameters in Eq (2.3) are showed in Table 1.

The integral form of Eq (2.3) is given by

where $X(t): = X(t, a), \ Y(t): = Y(t, a), \ X_0: = X(0, a),$ and $Y_0: = 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 $r_1, r_2$ 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, $r_i(a)\rightarrow r_i(a) + \sigma_i\dot{\omega}(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) = \int_0^A \gamma(t, a)X(t, a)da, Y(t, 0) = \int_0^A \beta(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$.

2.2.   Basic concepts

Throughout the paper, let $(\Omega, \mathfrak{F}, \{\mathfrak{F}_t\}_{t\geq 0}, \mathbb{P})$ be a complete probability space. Here, the filtration $\{\mathfrak{F}_t\}_{t\geq 0}$ satisfies the usual conditions (that is, it is increasing and right continuous with $\mathfrak{F}_0$ containing all $\mathbb{P}$-null sets), $\mathbb{E}$ denotes the expectation corresponding to $\mathbb{P}$. For a set $\mathbb{A}$, its indicator function is denoted by $\bf{1}_\mathbb{A} = \begin{cases} 1, &x\in \mathbb{A}, \\ 0, &x\notin \mathbb{A}. \end{cases}$ We also denote by $\mathbb{R}^n_+$ the positive cone in $\mathbb{R}^n$, that is $\mathbb{R}^n_+ = \{x\in \mathbb{R}^n: x_i > 0$ for all $1\leq i \leq n\}$. If $x\in \mathbb{R}^n$, its norms are denoted as $|x|_{\iota} = \begin{cases} |x_1|+|x_2|+\cdots+|x_n|, & \iota = 1, \\ (\sum^n_{i = 1}x^2_i)^{\frac{1}{2}}, & \iota = 2. \end{cases}$

Let $V\equiv\big\{\varphi| \varphi\in L^p([0, A]), \frac{\partial \varphi}{\partial a}\in L^p([0, A]),$ where $\frac{\partial \varphi}{\partial a}$ is generalized partial derivatives with respect to age a$\big\}$. $H = L^2([0, A])$ satisfy that $V\hookrightarrow H \equiv H'\hookrightarrow V'$, where $V'$ is the dual space of $V$, $H'$ is the dual space of $H$. The norm in $H$ is denoted as $|\cdot|$. The duality product between $V, V'$ is written as $\langle \cdot, \cdot \rangle$, the scalar product in $H$ is denoted by $(\cdot, \cdot)$. $\mathbb{C} = \mathbb{C}([0, T]; H)$ is the space of all continuous functions from $[0, T]$ into $H$. $I^p([0, T]; V)$ denotes the space of all $V$-valued processes $(P_t)_{t\in [0, T]}$, $L^p_V = L^p([0, T]; V)$. $W: = (I^p([0, T]; V)\bigcap L^2(\Omega; \mathbb{C}([0, T]; H)))\times (I^p([0, T]; V)\bigcap L^2(\Omega; \mathbb{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_{\alpha\rightarrow A^-}\int_{\alpha}^A \gamma(t, a)X(t, a)da = \lambda_x > 0$; $\lim_{\alpha\rightarrow A^-}\int_{\alpha}^A \beta(t, a)Y(t, a)da = \lambda_y > 0$.

(A2) The fertility rate of females $\gamma(t, a), \beta(t, a)\in \mathbb{C}([0, T]\times [0, A]; H)$.

(A3) $\alpha_{ij}(a)\in \mathbb{C}([0, A]; \mathbb{R}_{+}), i, j \in \{1, 2\}$ and satisfy $\sup_{a\in [0, A]} \alpha_{12}(a)\alpha_{21}(a) < \inf_{a\in [0, A]} \alpha_{11}(a)\alpha_{22}(a)$.

(A4) $r_i(a)\in \mathbb{C}\big([0, A]; \mathbb{R}\big)$ and $r_i(a)$ is nondecreasing with $-\infty < r_i(0) < 0 < r_i(A) < \infty$, for $i = 1, 2.$

(A5) $\sup_{a\in [0, A]}\alpha_{12}(a)\alpha_{21}(a) < 1$ and $r_1(0)(1-\alpha_{12} \alpha_{21})+\alpha_{12}(r_2(A)+\alpha_{21} r_1(A)) < 0$; $r_2(0)(1-\alpha_{12} \alpha_{21})+\alpha_{21}(r_1(A)+\alpha_{12} r_2(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 $\lambda_x$, which is a positive constant. This means that the species $X$ at age A has not lost its reproductive capacity.

(A2) $\gamma(t, a)$ and $\beta(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) $r_i(0)$ represents the intrinsic growth rate at age 0. "$r_i(0) < 0$" means that the birth rate is less than the death rate. "$0 < r_i(A) < \infty$" means that the birth rate is greater than the death rate at age $A$. Combined with the actual biological background, "$0 < r_i(A) < \infty$" also indicates that the species will not overreproduce at age A.

(A5) The value of $\alpha_{12}$ and $\alpha_{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 $-\frac{\alpha_{12}(r_2(A)+\alpha_{21}r_1(A))}{1-\alpha_{12}\alpha_{21}}$.

Theorem 2.1. Under the assumptions (A1)–(A3), for any given initial value $(X_0(a_0), Y_0(a_0))\in \mathbb{R}_+^2$, 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 $(X_0(a_0), Y_0(a_0))\in \mathbb{R}^2_+$, there is a unique local solution $(X(t, a), Y(t, a))\in W$ when $t\in [0, \tau_e)$, where $\tau_e$ is the explosion time of Eq (2.3). Now we need to proof that $\tau_e = \infty$, a.s. Choosing a positive constant $k_0 > 1$ such that for any initial value $(X_0(a_0), Y_0(a_0))$, we have $(X_0(a_0), Y_0(a_0))\in [\frac{1}{k_0}, k_0]\times[\frac{1}{k_0}, k_0]$. Then, for each $k\geq k_0$, $a \in [0, A]$, we define the stopping time $\tau_k$ as $\tau_k = \inf\{t \in [0, \tau_e): (X, Y)\notin (\frac{1}{k}, k)\times(\frac{1}{k}, k)\}.$

We can deduce that $\tau_k$ is increasing when $k\rightarrow \infty$ and $\lim_{k\rightarrow \infty}\tau_k = \tau_{\infty}$ with $\tau_{\infty}\leq \tau_e$. Assuming that there are two constants $T > 0$ and $\varepsilon\in (0, 1)$ such that $\mathbb{P}\{\tau_{\infty}\leq T\} > \varepsilon$, then, for the constant $k_0$, there is an integer $k_1$ such that

Now, for any positive constants $c_x, c_y$, we define a function $V: W\rightarrow \mathbb{R}^+$ as

we can check that $V(X, Y)\geq0$. By the Itô formula ^[15],

where $\Psi = \begin{pmatrix} X(t, a)\\ Y(t, a) \end{pmatrix},$ $\bar{R} = \begin{pmatrix} r_1(a)\\ r_2(a) \end{pmatrix}$, $\mbox{diag}(c_x, c_y) = \begin{pmatrix} c_x & 0\\ 0 & c_y \end{pmatrix}$, $A' = \begin{pmatrix} -\alpha_{11}(a) & \alpha_{12}(a)\\ \alpha_{21}(a) & -\alpha_{22}(a) \end{pmatrix}$, $\bar{C}^{\top} = \begin{pmatrix} c_x\\ c_y \end{pmatrix}$.

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

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

Now, let $M_2 = \frac{c_x\vee c_y}{c_x \wedge c_y}, M_3 = c_x \wedge c_y$, then

Since $|\Psi| = (X^2+Y^2)^{\frac{1}{2}}\leq X+Y$, we have

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

where $M_4 = M_1(5\vee \frac{2}{M_3})$. Then, we have

where $M_5 = V(X_0, Y_0)+M_4T$.

Applying the Gronwall inequality, we obtain

Since $\mathbb{P}(\Omega_k) = \mathbb{P}\{\tau_k\leq T\}\geq \varepsilon$ for $k\geq k_1$, for every $\omega \in \Omega_k, (t, a) \in Q$, $X(\tau_k, \omega)$ and $Y(\tau_k, \omega)$ equal either $k$ or $\frac{1}{k}$, hence

and

Taking $k\rightarrow \infty$, we have the contradiction that

The proof is completed.

3.   The approximate EM methods

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).

Equation (2.3) has four equilibria: $E_0(0, 0), \ E_1(\frac{r_1(A)}{\alpha_{11}}, 0), \ E_2(0, \frac{r_2(A)}{\alpha_{22}}),$ and $E_*(x^*, y^*)$, where

Now, let us first present some lemmas.

Remark 3.1. We note $H_{ix}(\Psi_j): = H_{ix}(X_j, Y_j), H_{iy}(\Psi_j): = H_{iy}(X_j, Y_j), H_{ix}(\psi_j): = H_{ix}(x_j, y_j), H_{iy}(\psi_j): = H_{iy}(x_j, y_j)$, $(i = 1, 2), (j = 1, 2)$ as follows.

Lemma 3.1. Under assumptions (A4) and (A5), for any $\psi_1, \psi_2\in W$, we have

and

where $L_1, \ L_2$ are constants and $\rho_1 = 2x^*-r_1(0)+\alpha_{12}(x^*+y^*), \ \rho_2 = 2y^*-r_2(0)+\alpha_{21}(x^*+y^*)$, $\psi_i^{\top}(t, a) = (x_i(t, a), y_i(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 $K_1$ and $K_2$ such that

and

where $\psi(t, a) = (x(t, a), y(t, a))^{\top}$ is the numerical solution corresponding to the exact solution $\Psi$ of Equation (2.3), $H_{1x}(\psi) = H_{1x}(x, y), H_{1x}(\psi) = H_{1x}(x, y)$

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

and

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

by the definition of $|\psi|$ and the positive of $|\psi|^2-|\psi|+1$, we derive that

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

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 $\Delta \in (0, 1)$, let $t_k = k\Delta$ for $k = 0, 1, 2, \cdots, [\frac{T}{\Delta}]$, where $[\frac{T}{\Delta}]$ denotes the integer part of $\frac{T}{\Delta}$. So, the discrete time EM approximate solution to Eq (2.3) is defined as

where $\Delta w_k = w_{t_{k+1}}-w_{t_k}$, $X(0) = X_0(a)$, and $Y(0) = Y_0(a)$.

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

where

are called the step processes. We can easily see that $\bar{x}(t) = x_k$ and $\bar{y}(t) = y_k$ for $t\in [t_k, t_{k+1})$ when $k = 0, 1, 2, \cdots, [\frac{T}{\Delta}]$. Without loss of generality, we denote $\psi(t) = (x(t), y(t))^{\top}$ and $\bar{\psi}(t) = (\bar{x}(t), \bar{y}(t))^{\top}$.

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

where $C$ only dependents on $T$ and $p$.

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

Since

where ${\check{\gamma}}$ means the maximum value of continuous function $\gamma(t, a)$. Then we have

Let

By Lemma 3.2 and Young inequality $a^\beta+b^{1-\beta}\leq\beta a+(1-\beta)b$, for $a, b \geq0$ and $\beta\in (0, 1)$, we have

Similarly, we can show that

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

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

where $C$ is a positive constant and $C_p$ is defined by

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

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

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

Similarly, we can get that

Since

where $\check{\beta}$ is the maximum value of continuous function $\beta(t, a)$, and

Let

By Lemma 3.2, we have

By Eq (3.6), we have

Moreover, using the Young inequality, we can get

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

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

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

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

Using Eqs (3.16), (3.24) and the norm of $\psi(t)$, we have

Applying the continuous Gronwall inequality, it yields

Now, let us present a theorem which shows that $\psi(t)$ and $\bar\psi(t)$ are coincide with the discrete solution.

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

Proof. For any $t\in [0, T]$, there exists a positive integer $k$ such that $t\in [t_k, t_{k+1}]$, and

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

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

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

By Theorem 3.1, we have

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).

4.   Convergence rate over the time interval $[0, T]$

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

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

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

hence we have

Using the Young inequality, we derive that

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

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

Similarly, we can show

Therefore, due to Young inequality,

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

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

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

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

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

5.   Numerical experiments

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

where $\alpha_{11}(a) = \frac{1}{(1-a)^2}, \ \alpha_{12}(a) = \cos^2a, \ \alpha_{21}(a) = \sin^2a, \ \alpha_{22}(a) = e^{-\frac{1}{a}}, \ \gamma(t, a) = \beta(t, a) = \frac{1}{1-a}, \ Q = (0, 1)\times(0, 1)$. All computations are run with $T = 1, \Delta = 0.005,$ and $\Delta 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.

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.

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.

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.

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

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 \in [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.).

6.   Conclusions and future research

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.

Acknowledgements

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).

Conflict of interest

The authors declare that they have no competing interests.