Numerical solvers for a poromechanic problem with a moving boundary

1.   Introduction

With the rapid development of industry and agriculture, the environment pollution has caused many serious ecological problems (see ^[1,2,3]), such as the reduction of species diversity and the extinction of some species. Therefore, it motivates many scholars' interest to study dynamic behavior of population in a polluted environment by establishing mathematical models. The population model in a polluted environment was first proposed by Hallam et al.^[4,5]. From then on, more investigations and discussions on the dynamic behavior of the deterministic population model can be found (see ^[6,7,8,9,10]). But in practical problems, population changes are affected not only by environmental noise but also by sudden changes of temperature and climate. Thus, several scholars have introduced random perturbations into population model to study dynamic behavior. For example, Liu and Wang ^[11] established the stochastic population model with impulsive toxicant input and obtained sufficient conditions on extinction, persistence, stability in the mean. Subsequently, Yu et al. ^[12] proved the existence of global positive solution for the stochastic population model with Allee effect under regime switching and established the threshold. In ^[13], Wei et al. proposed a stochastic population model with partial tolerance, discussed the conditions for population the extinction and proved the stationary distribution with ergodicity by constructing the Lyapunov function. Liu et al. ^[14] considered the significance of white noise and color noise on population persistence and extinction and studied stochastic population model with Markov switching. More research results on the persistence, extinction, and stability of random population models and others have been presented (see ^{[15,16,17,18,19]}). However, the above mentioned references didn't consider the invariant measure of population system with diffusion.

In fact, in the real world, the population and toxins in the ecology spread around the medium such as soil and water. In addition, we also know that the existence and uniqueness of invariant measure is one of the important properties for stochastic population model with Markov switching and diffusion. Nevertheless, if we introduce diffusion into stochastic population model, the corresponding Kolmogorov-Fokker-Planck (KFP) equation will become more complicated. Furthermore, the invariant measure of stochastic population model with Markov switching and diffusion is difficult to obtain. Therefore, it is of great significance to choose an effective numerical approximation method. To the best of our knowledge, the explicit Euler-Maruyama (EM) method has the advantages of easy calculation and small calculation amount. Motivated by ^[21,20], in the paper, we first develop a new stochastic population model with Markov switching and diffusion. Under suitable regularity assumptions, we discuss the existence and uniqueness of numerical invariant measure generated by the EM method. Subsequently, we prove that numerical invariant measure converges to the invariant measure of exact solution in the Wasserstein distance sense. In particular, the main contributions of the paper are as follows:

● We establish a novel stochastic population model with diffusion and Markov switching in a polluted environment. By using the Chebyshev's inequality, we obtain the existence and uniqueness of invariant measure for the model.

● Under local Lipschitz conditions, we study the approximation of numerical invariant measure generated by the EM method for the newly developed model.

The structure of this article is as follows: In Section2, we introduce some necessary preliminary knowledge results for the following analysis. In Section3, based on the Perron-Frobenius theorem, we study existence and uniqueness of invariant measure for the exact solution. In Section4, we mainly study the existence and uniqueness of numerical invariant measure for the EM scheme. In addition, we also prove that the numerical invariant measure of the EM scheme converges to invariant measure of exact solution. In Section5, the numerical example is given to verify our theoretical results. In Section6, we give the conclusions of this study.

2.   Preliminaries

In this paper, we introduce Markov switching and spatial diffusion into the model mentioned by Liu and Wang ^[15], and obtain the following model

where $\mathcal{L}: = (0, T)\times \Gamma$, $\Gamma$ is a bounded domain in $R^{3}$ with smooth boundary $\partial \Gamma$, $t\in (0, T)$; $X_{1}(t, x)$ denotes the population density at the location $x$ at time $t$. $X_{2}(t, x)$ is the concentration of toxicant in the organism at time $t$ and in spatial position $x$. The concentration of toxicant in the environment at the location $x$ at time $t$ is described by $X_{3}(t, x)$. $K(\Lambda_{t})$ is the net organismal uptake rate of toxicant from the environment at time $t$. $M(\Lambda_{t})$ is the total loss rate of the toxicant from the environment. $\mu(t, x, X_{2}(t, x), \Lambda_{t})$ denotes the decreasing rate function of the population at time $t$ and in spatial position $x$. $k_{i} > 0, i = 1, 2, 3$ is the diffusion coefficient. $\beta(t, x, X_{2}(t, x), \Lambda_{t})$ describes the intrinsic growth rate function of the population at time $t$ and in spatial position $x$. $u(t, x)$ denotes the exogenous total toxicant input into environment at time $t$ and in spatial position $x$. $l(\Lambda_{t})$ is the net organismal excretion rate of toxicant and $m(\Lambda_{t})$ is depuration rate of toxicant due to metabolic process and other losses.

Throughout the paper, Let $(V, \|\cdot\|)$ and $(H, |\cdot|)$ be two separable Hilbert spaces, with norm denoted by $\|\cdot\|$ and $|\cdot|$, respectively. $V$ is viewed as a subspace of $H$ with a continuous dense embedding. $V\Subset H$ represents the embedding is compact. $V'$ and $H'$ are the dual of $V$, $H$. We set $H_{3}: = H\times H \times H$. Let $(\Omega, \mathcal{F}, \mathbb{P})$ be a complete probability space with $\{\mathcal{F}_{t}\}_{0\leq t\leq T}$ the natural filtration generated by the Brownian motion $W_{t}$, which means $\mathcal{F}_{t} = \sigma\{W_{s}; 0\leq s\leq t\}$ augmented with all $\mathbb{P}$-null sets of $\mathcal{F}_{0}$. To construct such a filtration, we denote by $\mathcal{N}$ the collection of $\mathbb{P}$-null sets, i.e. $\mathcal{N} = \{B\in\mathcal{F}:\mathbb{P}(B) = 0\}$. In the paper, $C > 0$ represents different positive constants. Let $\Lambda_{t},$ $t > 0$, be a right-continuous Markov chain on the probability space taking values in a finite state $S = \{1, 2, \ldots, N\}$ for some positive integer $N < \infty$. The generator of $\{\Lambda_{t}\}_{t > 0}$ is specified by $Q = (q_{ij})_{N\times N}$, such that for a sufficiently small $\Delta$,

where $\Delta > 0$, $o(\Delta)$ satisfies $\lim_{\Delta\rightarrow 0}\frac{o(\Delta)}{\Delta} = 0$. Here $q_{ij}$ is the transition rate from $i$ to $j$ satisfying $q_{ii} = -\sum\limits_{i\neq j}q_{ij}$. We assume that the Markov chain $\{\Lambda_{t}\}$ defined on the probability space above is independent of the standard Brownian motion $\{W_{t}\}_{t\geq 0}$ and the $Q$ matrix is irreducible and conservative. Therefore, the Markov chain $\{\Lambda_{t}\}_{t\geq 0}$ has a unique stationary distribution $\pi: = (\pi_{1}, \ldots, \pi_{N})$ which can be determined by solving the linear equation

Let $\mathcal{P}(H_{3}\times\mathbb{S})$ stand for the family of all probability measures on $H_{3}\times\mathbb{S}$. For $\xi = (\xi_{1}, \xi_{2}, \xi_{3})^{\ast} \in H_{3}$, $\xi\gg \mathbf{0}$ means each component $\xi_{i} > 0$, $i = 1, 2, 3$.

Next, let's give some necessary assumptions:

$(\mathbb{H}1)$ Setting $X_{k, t}^{s_{k}, i}: = X_{k}^{s_{k}, i}(t, x)$, $k = 1, 2, 3$, there exists a positive constant $\rho_{i}$ such that for $i\in \mathbb S$, $(t, x)\in \mathcal{L}$

where $s_{1}$ and $\bar{s}_{1}$ are the different initial values of the first equation for system (2.1).

From $(\mathbb{H}1)$, for each $i\in \mathbb S$ and $X_{1, t}^{s_{1}, i}\in H$, we can obtain that for $(t, x)\in \mathcal{L}$

where $C$ depends on the initial value of the function $g(t, X_{1, t}^{s_{1}, i}, i)$.

$(\mathbb{H}2)$ For each $i\in \mathbb S$, there exist positive constants $\bar{M}$, $\bar{\beta}$ and $\bar{\mu}$ such that

$(\mathbb{H}3)$ $u(t, x)$ is non-negative measurable in $\mathcal{L}$, there exists a positive constant $\bar{u}$ such that

We replace $((X_{1, t}, X_{2, t}, X_{3, t}), \Lambda_{t})$ with $((X_{1, t}^{s_{1}, i}, X_{2, t}^{s_{2}, i}, X_{3, t}^{s_{3}, i}), \Lambda_{t}^{i})$, especially the initial value

For any $p\in (0, 1]$, we set $s: = (s_{1}, s_{2}, s_{3})$ and define a metric on $H_{3}\times\mathbb{S}$ as follows

where $I_{A}$ denotes the indicator function of the set $A$, and $\bar{s}: = (\bar{s}_{1}, \bar{s}_{2}, \bar{s}_{3})$ is different initial value. For $p\in (0, 1]$, we define the Wassertein distance between $\nu\in\mathcal{P}(H_{3}\times\mathbb{S})$ and $\nu'\in\mathcal{P}(H_{3}\times\mathbb{S})$ by

where the infimum is taken over all pairs of random variables $X_{k}$, $X_{k'}$ on $H_{3} \times\mathbb{S}$ with respective laws $\nu$, $\nu'$. Let $\mathbb{P}_{t}((s_{1}, s_{2}, s_{3}), i;\cdot)$ be the transition probability kernel of the pair $((X_{1, t}^{s_{1}, i}, X_{2, t}^{s_{2}, i}, X_{3, t}^{s_{3}, i}), \Lambda_{t}^{i})$, a time homogeneous Markov process (see ^[22]). Recall that $\pi\in \mathcal{P}(H_{3}\times\mathbb{S})$ is called an invariant measure of $((X_{1, t}^{s_{1}, i}, X_{2, t}^{s_{2}, i}, X_{3, t}^{s_{3}, i}), \Lambda_{t}^{i})$ if

holds. For any $p > 0$, let

where $\rho_{i}$ is introduced in the assumptions and $\gamma \in \rm{spec}(Q_{p})$, $\rm{spec}(Q_{p})$ denotes the spectrum of $Q_{p}$(i.e., the multi-set of its eigenvalues). $Re\gamma$ is the real part of $\gamma$ and ${\rm{diag}}(\rho_{1}, \ldots, \rho_{N})$ denotes the diagonal matrix whose diagonal entries are $\rho_{1}, \ldots, \rho_{N}$, respectively.

3.   Existence and uniqueness of invariant measures

In this section, we mainly prove the existence and uniqueness of the invariant measure for the exact solution, under the assumption conditions $(\mathbb{H}1)$–$(\mathbb{H}3)$. Firstly, in order to prove the existence and uniqueness of the underlying invariant measure, we prepare the following lemma.

Lemma 3.1. (see ^[22]) Let $N < \infty$ and assume further that

where $\mu_{i}$ is the stationary distribution of Markov chain $\{\Lambda_{t}\}_{t\geq 0}$, and $\rho_{i}$ is introduced in the assumption $(\mathbb{H}1)$. Then

$(1)$ $\eta_{p} > 0$ if $\max\limits_{i\in \mathbb{S}}\rho_{i}\leq0$;

$(2)$ $\eta_{p} > 0$ for $p < \max\limits_{i\in \mathbb{S}, \rho_{i} > 0}\{-2q_{ii}/\rho_{i}\}$ if $\max\limits_{i\in \mathbb{S}}\rho_{i} > 0$.

Remark 1: The system (2.1) is said to be attractive " in average " if Eq (3.1) holds. The Lemma 3.1 provides great convenience to study the existence and uniqueness of invariant measure for exact solution, i.e., the proof of Theorem 3.1.

Theorem 3.1. Let $N < \infty$ and assume further that $(\mathbb{H}1)$–$(\mathbb{H}3)$ hold with $\max_{i\in S}\rho_{i} > 0$. Then the exact solution of system (2.1) admits a unique invariant measure $\pi\in \mathcal{P}(H_{3}\times\mathbb{S})$.

Proof. The key point of proof is to divide the whole proof into two parts of existence and uniqueness.

$(I)$ Existence of invariant measure. Let $((Y^{s_{1}, i}_{1, t}, Y^{s_{2}, i}_{2, t}, Y^{s_{3}, i}_{3, t}), \Lambda_{t}^{i})$ be the exact solution of system (2.1) with $((s_{1}, s_{2}, s_{3}), i)$ as initial value, where $((s_{1}, s_{2}, s_{3}), i)\in H_{3}\times \mathbb{S}$. A simple application of the Feynman-Kac formula show that let $Q_{p, t} = e^{tQ_{p}}$, where $Q_{p}$ is given in Eq (2.8). Then, the spectral radius Ria$(Q_{p, t})$ (i.e., Ria$(Q_{p, t}) = \sup_{\lambda\in \rm{spec}(Q_{p, t})}|\lambda|$) of $Q_{p, t}$ equals to $e^{-\eta_{p}t}.$ Since all coefficients of $Q_{p, t}$ are positive, by the Perron-Frobenius theorem (see ^[23]) yields that $-\eta_{p}$ is a simple eigenvalue of $Q_{p}$, all other eigenvalues have a strictly smaller real part. Note that the eigenvector of $Q_{p, t}$ corresponding to $e^{-\eta_{p}t}$ is also an eigenvector of $Q_{p}$ corresponding to $-\eta_{p}$. According to Perron-Frobenius theorem, for $Q_{p}$ it can be found that there is a positive eigenvector $\xi^{(p)} = (\xi_{1}^{(p)}, \ldots, \xi_{N}^{(p)})\gg {\mathbf{0}}$ corresponding to the eigenvalue $-\eta_{p}$, and $\xi^{(p)}\gg \mathbf{0}$ means that each component $\xi_{i}^{(p)} > 0$. Let

where $1\wedge \min\limits_{i\in\mathbb{S}, \rho_{i} > 0}\{-2q_{ii}/\rho_{i}\}: = \min\{1, \min\limits_{i\in\mathbb{S}, \rho_{i} > 0}\{-2q_{ii}/\rho_{i}\}\}$. Combined with Lemma 3.1, we can get

In order to investigate the existence and uniqueness of invariant measure for exact solution, we need to prove the boundedness of exact solution for system (2.1). In other words, we need to prove whether the following inequality holds.

First, using the It$\hat{\rm{o}}$'s formula (see ^[24], Theorem 1.45 of p.48), we can have

Using $p(p-2)/2 < 0$, due to $p\in (0, p_{0}),$ and combining with the following inequality,

where $0\leq k_{0}\leq k(t, x) < \infty$ ($k_{0}$ is a constant). Further, we have

Therefore, based on assumption conditions $(\mathbb{H}1)$–$(\mathbb{H}3)$ and the inequality $2ab\leq \varepsilon a^{2}+\frac{1}{\varepsilon}b^{2}$, $\varepsilon > 0$ we can obtain

where $\bar{K}: = \max\limits_{i}\{K(\Lambda_{i})\}$, for all $i\in\mathbb{S}$, $0 < \bar{K} < \infty$. Then, setting $C_{1}: = 2(\bar{\beta}-\mu_{0})+\rho_{0}$, $\rho_{0}: = \max\limits_{i\in \mathbb{S}}|\rho_{\Lambda_{\epsilon}^{i}}|$ $C_{2}: = \max\limits_{i\in \mathbb{S}}2(l(\Lambda_{\epsilon}^{i})+m(\Lambda_{\epsilon}^{i}))+\varepsilon_{1}\bar{K}^{2}$ and $C_{3}: = 2\bar{M}+\frac{1}{\varepsilon_{2}}+\frac{1}{\varepsilon_{1}}$, $C_{4}: = \varepsilon_{2}+c$ are different constants and using the inequality

we can further estimate

Finally, by the Gronwall's lemma, we can get the result

and further estimates can be obtained as follows

For $\forall$ $t > 0$, we can define a probability measure

Then, let $Y_{t}^{s, i}: = (Y_{1, t}^{s_{1}, i}, Y_{2, t}^{s_{2}, i}, Y_{3, t}^{s_{3}, i})$, for any $\varepsilon > 0$, by Eq (3.7) and Chebyshev's inequality, there exists an $r > 0$ sufficiently large such that

Hence, $\chi_{t}$ is tight since the compact embedding $V\Subset H$, then $K_{r} = \{s\in H_{3}; |s|\leq r\}$ is a compact subset of $H_{3}$ (see ^[25], Definition 2, p.27) for each $i\in \mathbb{S}$. Combined with the Fellerian property of transition seimgroup for $\mathbb{P}_{t}(s, i;\cdot)$ and according to Krylov-Bogoliubov theorem (see ^[26]), $((Y_{1, t}^{s_{1}, i}, Y_{2, t}^{s_{2}, i}, Y_{3, t}^{s_{3}, i}), \Lambda_{t}^{i})$ has an invariant measure (see ^[27]). Next, we prove the uniqueness of the invariant measure for $((Y_{1, t}^{s_{1}, i}, Y_{2, t}^{s_{2}, i}, Y_{3, t}^{s_{3}, i}), \Lambda_{t}^{i})$.

$(II)$ Uniqueness of invariant measure. First, let $((Y_{1, t}^{s_{1}, i}, Y_{2, t}^{s_{2}, i}, Y_{3, t}^{s_{3}, i}), \Lambda_{t}^{i})$ and $((Y_{1, t}^{\bar{s}_{1}, i}, Y_{2, t}^{\bar{s}_{2}, i}, Y_{3, t}^{\bar{s}_{3}, i}), \Lambda_{t}^{i})$ be the solutions of the system (2.1) satisfying the initial values $((s_{1}, s_{2}, s_{3}), i)$ and $((\bar{s}_{1}, \bar{s}_{2}, \bar{s}_{3}), i)$, respectively. Under assumption conditions $(\mathbb{H}1)$–$(\mathbb{H}3)$, we take $\forall\varepsilon \in (0, 1)$ and use It$\hat{\rm{o}}$'s formula, combined with Eq (3.4), we have

where $C_{i}$, $i = 1, 2, 3$ have been explained before and $\rho_{\Lambda_{\epsilon}^{i}}$ is introduced in the assumption $(\mathbb{H}1)$. In addition, using the result of Eqs (3.5) and (3.6), we can get

when $\varepsilon \rightarrow 0$, we can get the following result

Define the stopping time

According to the definition of $\mathbb{S}$ and irreducibility of $Q$, there exists $\theta > 0$ such that

Due to $p\in (0, p_{0})$, and choose $q > 1$ such that $0 < pq < p_{0}$, where $p_{0}$ is introduced in Eq (3.2). Using H$\ddot{o}$lder's inequality, we can have

Applying the result of Eq (3.11), we further obtain

where $\sigma: = \frac{(q-1)\theta }{2q}\wedge \frac{\eta_{p}}{2}$, and in the last step, it follows from Eqs (3.7) and (3.10) such that

and

Thus, we also have assertion

Then, according to Eq (3.11), we can get

Next, according to Eqs (3.14) and (3.13) that

where $\sigma^{*}: = \sigma\wedge\theta$. Assume $\pi$, $\nu\in \mathcal{P}(H_{3}\times\mathbb{S})$ are invariant measures of $((Y_{1, t}^{s_{1}, i}, Y_{2, t}^{s_{2}, i}, Y_{3, t}^{s_{3}, i}), \Lambda_{t}^{i})$, it follows from Eq (3.15) that

When $t\rightarrow \infty$, we find $W_{p}(\pi, \nu) \rightarrow 0$. Hence, uniqueness of invariant measure follows immediately. The proof of Theorem 3.1 has been completed.

In the following section, we will investigate existence and uniqueness of numerical invariant measure and prove the convergence of numerical invariant measure.

4.   Numerical invariant measure

In this section, we mainly discuss existence and uniqueness of numerical invariant measure for system (2.1) under the assumption conditions $(\mathbb{H}1)$–$(\mathbb{H}3)$. In order to facilitate the discussion, we consider the numerical solution in the discrete-time for system (2.1). For a given step size $\delta \in (0, 1)$, we define the discrete-time Euler-Maruyama (EM) scheme associated with model (2.1) as follows

where $n\geq 0$ and $\Delta W_{n}\triangleq W_{(n+1)\delta}-W_{n\delta}$ denotes Brownian motion increment, $\Delta \bar{X}_{k, n\delta}^{s_{k}, i}$ is the Laplace of $\bar{X}_{k, n\delta}^{s_{k}, i}$, with the initial data $((\bar{X}_{1}^{0}, \bar{X}_{2}^{0}, \bar{X}_{3}^{0}), \Lambda_{0}) = ((s_{1}, s_{2}, s_{3}), i)\in H_{3}\times\mathbb{S}$ which is introduced before. Equations (4.1) and (4.2) are the discrete-time EM scheme and continuous-time EM scheme of the corresponding system (2.1), respectively. For convenience, we define the corresponding approximate solution to the system (2.1) on continuous time.

where $t > 0, \Lambda_{0}^{i} = i\in \mathbb{S}$, $\forall b\geq 0$, $\lfloor b\rfloor$ is the interger part of $b$. Obviously, by a straightforward calculation, we can have $(X_{1, \lfloor \epsilon/\delta\rfloor\delta}^{s_{1}, i}, X_{2, \lfloor \epsilon/\delta\rfloor\delta}^{s_{2}, i}, X_{3, \lfloor \epsilon/\delta\rfloor\delta}^{s_{3}, i}) = (\bar{X}_{1, \lfloor \epsilon/\delta\rfloor\delta}^{s_{1}, i}, \bar{X}_{2, \lfloor \epsilon/\delta\rfloor\delta}^{s_{2}, i}, \bar{X}_{3, \lfloor \epsilon/\delta\rfloor\delta}^{s_{3}, i}).$

Let $\mathbb{P}_{n\delta}^{\delta}((s_{1}, s_{2}, s_{3}), j; \cdot)$ be the transition probability kernel of $((\bar{X}_{1, n\delta}^{s_{1}, i}, \bar{X}_{2, n\delta}^{s_{2}, i}, \bar{X}_{3, n\delta}^{s_{3}, i}), \Lambda_{n\delta}^{i})$. If $\pi^{\delta}\in\mathcal{P}(H_{3}\times\mathbb{S})$ satisfies the following equation

then we call $\pi^{\delta}\in\mathcal{P}(H_{3}\times\mathbb{S})$ an invariant measure of $((\bar{X}_{1, n\delta}^{s_{1}, i}, \bar{X}_{2, n\delta}^{s_{2}, i}, \bar{X}_{3, n\delta}^{s_{3}, i}), \Lambda_{n\delta}^{i})$ or a numerical invariant measure of $((X_{1, t}^{s_{1}, i}, X_{2, t}^{s_{2}, i}, X_{3, t}^{s_{3}, i}), \Lambda_{t}^{i})$. Let

Our main result in this section is as follows

Lemma 2. Under the conditions of Lemma 3.1 and combining Eq (3.2) with (3.3), it holds that

for any $p\in(0, p_{0})$, $(s, i) = ((s_{1}, s_{2}, s_{3}), i)$, $(\bar{s}, j) = ((\bar{s}_{1}, \bar{s}_{2}, \bar{s}_{3}), j) \in H_{3}\times\mathbb{S}.$ $p_{0}$ is given in Eq (3.2).

Lemma 4.1 shows that numerical solution $(\bar{X}_{1, n\delta}^{s_{1}, i}, \bar{X}_{2, n\delta}^{s_{2}, i}, \bar{X}_{3, n\delta}^{s_{3}, i})$ tends to $(\bar{X}_{1, n\delta}^{\bar{s}_{1}, j}, \bar{X}_{2, n\delta}^{\bar{s}_{2}, j}, \bar{X}_{3, n\delta}^{\bar{s}_{3}, j})$ when $n\rightarrow \infty$ and $\delta \rightarrow 0$ under different initial values and states. This lemma provides a great convenience for the proof of Theorem 4.1. Applying a method similar to Theorem 3.1 can prove the conclusion of Lemma 4.1, so it is omitted.

Theorem 4.1. Under the conditions of Theorem 3.1, there exists a sufficiently small $\delta^{*}$ such that for any $\delta\in (0, \delta^{*})$, the solutions of the EM method (4.2) converge to a unique invariant measure $\pi^{\delta}\in \mathcal{P}(H_{3}\times \mathbb{S})$ with some exponential rate $\bar{\gamma} > 0$ in the Wassertein distance.

Proof. In fact, for any the initial data $(s_{1}, s_{2}, s_{3})$, by Eq (4.2) and the Chebyshev's inequality, we derive that $\{\delta_{(s_{1}, s_{2}, s_{3})}\mathbb{P}^{\delta}_{n\delta}\}$ is tight. Therefore, there exists an exact subsequence which converges weakly to an invariant measure denoted by $\pi^{\delta}\in\mathcal{P}(H_{3}\times \mathbb{S})$. According to the Eq (3.14), we have the following result

For any $n > 0,$ combining with Eq (4.4), it is not difficult to get

where $\bar{\gamma}: = \varrho\wedge \theta$, and using the Kolmogorov-Chapman equation and Eq (4.6), for any $n, m > 0$, we have

where $H_{1} = \mathbb{P}_{m\delta}^{\delta}((s_{1}, s_{2}, s_{2}), i;d(\bar{s}_{1}, \bar{s}_{2}, \bar{s}_{3}), j)$, then taking $m\rightarrow \infty$ such that

in other words, $\pi^{\delta}$ is the unique invariant measure of $\{\delta_{(s_{1}, s_{2}, s_{3})}\mathbb{P}^{\delta}_{n\delta}\}$. $\forall \pi^{\delta}, \nu^{\delta}\in \mathcal{P}(H_{3}\times\mathbb{S})$ are invariant measures of $((\bar{X}_{1, n\delta}^{s_{1}, i}, \bar{X}_{2, n\delta}^{s_{2}, i}, \bar{X}_{3, n\delta}^{s_{3}, i}), \Lambda_{n\delta}^{i})$ and $((\bar{X}_{1, n\delta}^{\bar{s}_{1}, j}, \bar{X}_{2, n\delta}^{\bar{s}_{2}, j}, \bar{X}_{3, n\delta}^{\bar{s}_{3}, j}), \Lambda_{n\delta}^{j})$, respectively. Further, we have

The uniqueness for the numerical invariant measure have been completed. Therefore, the proof of Theorem 4.1 is complete. To show that the numerical invariant measure $\pi^{\delta}$ converges to the invariant measure of the corresponding exact solution under the Wasserstein distance, the following theorem is given.

Theorem 4.2. Under the assumptions of Theorem 4.1 and Eq (4.8), for $\delta\in(0, 1)$ there exists $C > 0$ such that

where $p_{0} > 0$ is defined in Eq (3.2).

Proof. For $p\in (0, p_{0})$, due to

and

Then based on the assumption conditions of $(\mathbb{H}1)$–$(\mathbb{H}3)$ and Eq (4.8), there exists a sufficiently small $\delta^{*}$ such that for any $\delta\in (0, \delta^{*})$, there is $n > 0$ sufficiently large such that

For fixed $n > 0$ and using the triangle inequality, and by the similar way of ^[22], we can obtain $\lim\limits_{\delta\rightarrow 0}W_{p}(\delta_{((s_{1}, s_{2}, s_{2}), i)}\mathbb{P}_{n\delta}, \delta_{((s_{1}, s_{2}, s_{2}), i)}\mathbb{P}_{n\delta}^{\delta}) = 0$. In other words, there exists a positive constant $\bar{\nu}$ such that $W_{p}(\delta_{((s_{1}, s_{2}, s_{2}), i)}\mathbb{P}_{n\delta}, \delta_{((s_{1}, s_{2}, s_{2}), i)}\mathbb{P}_{n\delta}^{\delta})\leq Ce^{\bar{\nu} \delta n}\delta^{\frac{p}{2}}$. According to Theorem 3.1 and Eq (4.8), we can get the following result

where $\gamma^{*}: = \sigma^{*}\wedge\bar{\gamma}.$ Let $\bar{C}$ be the integer part of constant $-p\ln\delta/[2(\bar{\nu}+\gamma^{*})\delta]$, obviously, $\bar{C}\rightarrow 0$ as $\delta \rightarrow 0$. On the other hand, we have $e^{\bar{\nu} \bar{C}\delta}\delta^{\frac{p}{2}}\leq \delta^{\frac{p\sigma^{*}}{2(\bar{\nu}+\gamma^{*})}}\leq\delta^{\frac{p}{2}}$, $e^{-\sigma^{*}\bar{C}\delta}\leq e^{\gamma^{*}\delta^{*}}\delta^{\frac{p}{2}}$. Therefore, $W_{p}(\pi, \pi^{\delta})\leq C\delta^{\frac{p}{2}}$ holds.

Further, to illustrate the validity of our theory which are discussed in the previous section, we will give a numerical example.

5.   Numerical examples

Let $\Lambda_{t}$ be a Markov chain with the state space $\mathbb{S} = \{1, 2\}$, and the generator

It is easy to show that its unique stationary distribution $\pi = (\pi_{1}, \pi_{2})$ is given by $\pi_{1} = 1/2$, $\pi_{1} = 1/2$. On the other hand, we give the following setting: $V(\Lambda_{t}): = l(\Lambda_{t})+m(\Lambda_{t})$, when $\Lambda_{t} = 1$, we choose $M(1) = \frac{1}{2}\exp(\frac{2}{1+2t})$, $K(1) = 0.01\sin(\frac{1}{(3+0.2t)^2})$ and $V(1) = l(1)+m(1) = 1.99$; when $\Lambda_{t} = 2$, we choose $M(2) = \frac{9}{10}(\frac{1}{1+t})$, $K(2) = 0.05\sin(\frac{1}{(3+0.2t)^2})$ and $V(2) = l(2)+m(2) = 1.6$. In the state $1$ and $2$, setting $T = 1$, $t\in (0, 1)$, $\beta: = \beta(t, x, X_{2}(t, x), \Lambda_{t}) = \frac{1}{2}(1-\frac{0.5X_{2}(t, x)}{0.5(1+X_{2}(t, x))})(1-\frac{x}{5+x})$, $\mu: = \mu(t, x, X_{2}(t, x), \Lambda_{t}) = \frac{3}{10}(0.5-\frac{0.8X_{2}(t, x)}{1+0.5X_{2}(t, x)})(1-\frac{x}{0.5+x})$, $g: = g(t, x, X_{1}(t, x), \Lambda_{t}) = 0.05+0.3X_{1}(t, x)$, and taking $k_{1} = 0.005$, $k_{2} = k_{3} = 0.05$, $s_{1}(x) = s_{2}(x) = \frac{0.2}{(1+x)^2}$, $s_{3}(x) = \frac{0.2}{(1+1.5x)^2}$, the system (2.1) is described as follows

First, for the system (5.1), we use the discrete-time EM method for numerical simulation. Figure 1 is a simulation of Markov chain which describes switching between different states.

Then, taking $T = 1$, $N = 100$, $|W_{i+1}-W_{i}| = \sqrt{\delta }$ and $t\in(0, 1)$, step sizes $\delta = 0.005$. Among them, the values of $X_{3}(t, x)$ and $X_{2}(t, x)$ do not exceed 0.4. This satisfies the practical significance, i.e., $0\leq X_{2}(t, x)\leq 1$, $0\leq X_{3}(t, x)\leq1$.

As far as we know, the exact solution for system (5.1) is difficult to find. Inspired by ^[30] and based on the method of ^[28], we can take the "explicit solution" $Y_{1}(t, x) = \exp(\frac{1}{2}-\frac{1}{1-x}-\frac{t^2}{2})(1+\Delta W)$, $Y_{2}(t, x) = K\int_{0}^{t}Y_{3}(t, x)\exp\{(l+m)(s-t)\}ds+C_{Y_{2}}\exp\{-(l+m)t\}$ and $Y_{3}(t, x) = \int_{0}^{t}u(t, x)\exp\{h(s-t)\}ds+C_{Y_{3}}\exp\{-ht\}$ replace exact solution, where $C_{Y_{2}}, C_{Y_{3}}$ are initial values of $Y_{2}$ and $Y_{3}$, respectively. Setting $C_{Y_{3}} = \frac{0.2}{(1+1.5x)^2}$, $C_{Y_{2}} = \frac{0.2}{(1+x)^2}$, $K = 0.05$, $l+m = 1.9$ and $h = 0.5$, $u(t, x) = \frac{1}{5}(\frac{4}{(1+2x)^2}-\frac{1}{2})(1-t)^3$. Then, The simulation results are presented separately in Figure 2(a), Figure 4(a) and Figure 6(a). In Figure 6(b) and Figure 4(b) reflect the numerical simulation of $X_{3}(t, x)$ and $X_{2}(t, x)$ with Markov switching when the step size is 0.005 under the state $"1$" and "$2$" switching.

In addition, Figure 7, Figure 5 and Figure 3 show mean-square error between "explicit solutions $Y_{3}$, $Y_{2}$ and $Y_{1}$" and the corresponding numerical solutions $X_{3}$, $X_{2}$ and $X_{1}$(Figure 6(b), Figure 4(b) and Figure 2(b)) of stochastic population with diffusion and Markov switching in a polluted environment system (5.1), when we take step sizes $\delta = 0.005, 0.0001$. Obviously, when the step size $\delta$ changes from 0.005 to 0.0001, the error values decreases from 0.14, 0.4 and 0.04 to 0.012, 0.025 and 0.02, respectively. Combining Figure 7, Figure 5 and Figure 3, we have the assertion that the smaller the step size, the smaller the error. Hence, it is not difficult to conclude that when $\delta \rightarrow 0$, the numerical solution $X_{3}(t, x)$, $X_{2}(t, x)$, $X_{1}(t, x)$ under discrete-time EM method converges to the explicit solution $Y_{3}(t, x)$, $Y_{2}(t, x)$, $Y_{1}(t, x)$, respectively.

6.   Concluding remarks

In this paper, we establish a new stochastic population model with Markov chain and diffusion in a polluted environment. Based on the Perron-Frobenius theorem, when the diffusion coefficient satisfies the local Lipschitz, the criterion on the existence and uniqueness of invariant measure for the exact solution is given. Moreover, we also discuss the existence and uniqueness of numerical invariance measure for model (2.1) under the discrete-time Euler-Maruyama scheme, and prove that numerical invariance measure converges to invariance measure of the corresponding exact solution in the Wasserstein distance sense. At the end of this paper, the accuracy of the theoretical results is verified by numerical simulation.

Acknowledgements

The authors are very grateful to the anonymous reviewers for their insightful comments and helpful suggestions. The research was supported by the Natural Science Foundation of China (Grant numbers 11661064). This research was funded by the "Major Innovation Projects for Building First-class Universities in China's Western Region" (ZKZD2017009).

Conflict of interest

The authors declare there is no conflict of interest.