An ordinary differential equation system is called a Poisson system^[1], if it can be written in the following form

$ \mathrm{d} \boldsymbol{y}(t)=\boldsymbol{B}(\boldsymbol{y}(t)) \nabla H(\boldsymbol{y}(t)) \mathrm{d} t $

(1)

where $\boldsymbol{y} \in \mathbb{R}^n, \boldsymbol{B}(\boldsymbol{y})=\left(b_{i j}(\boldsymbol{y})\right)_{n \times n}$ is a smooth skew-symmetric matrix-valued function satisfying

$\sum\limits_{l=1}^n\left(\frac{\partial b_{i j}(\boldsymbol{y})}{\partial y_l} b_{l k}(\boldsymbol{y})+\frac{\partial b_{j k}(\boldsymbol{y})}{\partial y_l} b_{l i}(\boldsymbol{y})+\frac{\partial b_{k i}(\boldsymbol{y})}{\partial y_l} b_{l j}(\boldsymbol{y})\right)=0 $

for all i, j, k=1, …, n, and H is a smooth function. A function $C(\boldsymbol{y})$ is called a Casimir function of the Poisson system (1), if $\forall \boldsymbol{y}, \nabla C(\boldsymbol{y})^{\mathrm{T}} \boldsymbol{B}(\boldsymbol{y})=\boldsymbol{0}$.

The Lotka-Volterra (L-V) model of systems with n interacting components is given by

$ \dot{y}_i=y_i\left(b_i+\sum\limits_{j=1}^n a_{i j} y_j\right), \quad i=1, 2, \cdots, n, $

(2)

where a_ij, b_i(i, j=1, 2, …, n) are real parameters. In Ref. [2], the Poisson structure of a class of Lotka-Volterra systems was analyzed, which can be written in the form of (1) with

$ \begin{aligned} &H(\boldsymbol{y})=\sum\limits_{i=1}^n \beta_i y_i-p_i \ln y_i, \\ &\boldsymbol{B}(\boldsymbol{y})=\left(b_{i j} y_i y_j\right)_{n \times n}= \\ &\operatorname{diag}\left(y_1, \cdots, y_n\right) \boldsymbol{B} \operatorname{diag}\left(y_1, \cdots, y_n\right), \end{aligned} $

(3)

where $\boldsymbol{y}(t)=\left(y_1(t), \cdots, y_n(t)\right)^{\mathrm{T}}, \boldsymbol{B}=\left(b_{i j}\right)_{n \times n}$ is a skew-symmetric constant matrix, and β_i≠0 (i=1, …, n). It is not difficult to check that (1) with (3) can be of the form (2).

As is well known, however, stochastic perturbations are unavoidable and universal in real world. Accordingly, there have been literature exploring the Lotka-Volterra system (2) under Gaussian white noise perturbations, which were written in stochastic differential equations (SDEs) of Itô form^[3-6]. For instance, Mao et al.^[4] showed that the environmental noise suppresses the explosion in the Lotka-Volterra system, i.e., with probability 1 the solution of the stochastic L-V system exists (no explosion in finite time) uniquely under certain non-global Lipschitz conditions; Rudnicki and Pichór^[6] analyzed the influence of various stochastic perturbations on prey-predator systems. To the best of our knowledge, however, there is few research on Lotka-Volterra systems under stochastic perturbations of Stratonovich sense, except for the work by Khasminskii and Klebaner^[7], which investigated the long term behavior of the solution of the Lotka-Volterra systems under small random perturbations of Stratonovich sense on the birth and death rate.

Poisson systems under certain Stratonovich white noises perturbations, namely the stochastic Poisson systems, got attention in recent years, see e.g. Refs. [8-11], where in Ref. [9], the general form of stochastic Poisson systems was given as

$ \begin{gathered} \mathrm{d} \boldsymbol{y}(t)=\boldsymbol{B}(\boldsymbol{y}(t))\left(\nabla H_0(\boldsymbol{y}(t)) \mathrm{d} t+\right. \\ \left.\sum\limits_{r=1}^m \nabla H_r(\boldsymbol{y}(t)) \circ \mathrm{d} W_r(t)\right), \end{gathered} $

(4)

where $\boldsymbol{B}(\boldsymbol{y}), H_0(\boldsymbol{y})$ are defined the same way as for $\boldsymbol{B}(\boldsymbol{y})$ and $H(\boldsymbol{y})$, respectively, for the deterministic Poisson systems (1), and $H_r(\boldsymbol{y})(r=1, \cdots, m)$ are smooth functions. $\left(W_1(t), \cdots, W_m(t)\right)$ is an m-dimensional standard Wiener process defined on a complete filtered probability space, and the circle '$\circ$' in front of $\mathrm{d} W_r(t)$ denotes Stratonovich stochastic differential equations.

In this paper we consider the Lotka-Volterra systems (1) with (3) under Stratonovich white noise perturbation, of the following form:

$ \mathrm{d} \boldsymbol{y}(t)=\boldsymbol{B}(\boldsymbol{y}(t)) \nabla H(\boldsymbol{y}(t))(\mathrm{d} t+c \circ \mathrm{d} W(t)), $

(5)

where c>0 denotes the intensity of random noise, and W(t) is a one-dimensional standard Wiener process. Obviously it is a stochastic Poisson system according to (4). The stochastic version of the Lotka-Volterra system studied as an example in Ref. [8] is of this form. However, it can be seen that the coefficients of the system do not satisfy the global Lipschitz and linear growth conditions guaranteeing existence and uniqueness of the solution^[12-13]. When transfered into its equivalent Itô form, the drift coefficient is cubic, which is different from the equations with non-global Lipschitz and non-linearly growing coefficients considered in the Refs. [3-7], wherefore the results regarding well-posedness of the solutions in these literature are not applicable for this system either. It then arises the question, whether the system (5) with (3) has a unique solution, which does not explode in a finite time, with probability 1. Moreover, the definiton of $H(\boldsymbol{y})$ in (3) intrinsically needs y_i>0 (i=1, …, n). Then it is also to prove that the solution of (5) with (3) is positive with probability 1. Up to now, to our knowledge, no answers are available for these questions.

1 A class of invariants of the system

We first introduce some notations. To write a vector $\boldsymbol{\alpha}>0$ is to mean that all its components are greater than 0. $ \mathbb{R}_{+}^n:=\left\{\boldsymbol{y} \in \mathbb{R}^n: \boldsymbol{y}>0\right\}$, and we denote the kernal of the constant matrix $\boldsymbol{B}$ by ${Ker} \boldsymbol{B}:=\left\{\boldsymbol{x} \in \mathbb{R}^n: \boldsymbol{B} \boldsymbol{x}=\boldsymbol{0}\right\} . a \wedge b:=\min \{a, b\}$, and we let $\left(\varOmega, \mathscr{F}, \left\{\mathscr{F}_t\right\}_{t \geqslant 0}, P\right)$ be a complete probability space with the filtration $\left\{\mathscr{F}_t\right\}_{t \geqslant 0}$ satisying the usual conditions, where the one-dimensional Wiener process W(t) is defined.

It is not difficult to verify that

$ C(\boldsymbol{y})=\alpha_1 \ln y_1+\cdots+\alpha_n \ln y_n $

for $\boldsymbol{y} \in \mathbb{R}_{+}^n$ is a Carsimir function of system (5) with (3), i.e., $\nabla C(\boldsymbol{y})^{\mathrm{T}} \boldsymbol{B}(\boldsymbol{y})=\boldsymbol{0}\left(\forall \boldsymbol{y} \in \mathbb{R}_{+}^{\mathrm{n}}\right)$, whenever $\boldsymbol{\alpha}=\left(\alpha_1, \cdots, \alpha_n\right)^{\mathrm{T}} \in \operatorname{Ker} \boldsymbol{B}$. Note that, we do not exclude the case α= 0.

Proposition 1.1 Let $T>0, f\left(x_1, x_2\right)$ be a binary function defined on $\mathbb{R} \times \mathbb{R}$, and $f \in C^1(\mathbb{R} \times$ $\mathbb{R}$. Suppose the solution $\boldsymbol{y}(t)$ of system (5) with (3) is positive on [0, T]. Then $f(H(\boldsymbol{y}(t))$, $C(\boldsymbol{y}(t)))$ is an invariant of the system (5) with (3) on [0, T], where $C(\boldsymbol{y})$ is the Casimir function mentioned above.

Proof By the Stratonovich chain rule, it holds on [0, T]:

$ \begin{aligned} &\mathrm{d} H(\boldsymbol{y}(t))=\nabla H(\boldsymbol{y}(t))^{\mathrm{T}} \circ \mathrm{d} \boldsymbol{y}(t) \\ &=\nabla H(\boldsymbol{y}(t))^{\mathrm{T}} \boldsymbol{B}(\boldsymbol{y}(t)) \nabla H(\boldsymbol{y}(t)) \\ &(\mathrm{d} t+c \circ \mathrm{d} W(t))=0 \\ & \end{aligned} $

where the last equality is due to skew-symmetry of $\boldsymbol{B}(\boldsymbol{y})$, and

$ \begin{aligned} &\mathrm{d} C(\boldsymbol{y}(t)) \\ &=\nabla C(\boldsymbol{y}(t))^{\mathrm{T}} \circ \mathrm{d} \boldsymbol{y}(t) \\ &=\nabla C(\boldsymbol{y}(t))^{\mathrm{T}} \boldsymbol{B}(\boldsymbol{y}(t)) \nabla H(\boldsymbol{y}(t))(\mathrm{d} t+c \circ \mathrm{d} W(t)) \\ &=0 . \end{aligned} $

Thus

$ \begin{aligned} &\mathrm{d} f(H(\boldsymbol{y}(t)), C(\boldsymbol{y}(t)))= \\ &\partial_1 f(H(\boldsymbol{y}), C(\boldsymbol{y})) \circ \mathrm{d} H(\boldsymbol{y}(t))+\partial_2 f(H(\boldsymbol{y}) \\ &C(\boldsymbol{y})) \circ \mathrm{d} C(\boldsymbol{y}(t))=0 . \end{aligned} $

2 Non-explosion and positiveness of the solution

In the following, we will prove that the solution of the system (5) with (3) is globally non-explosive and positive almost surely. To this end, we make the following assumptions.

Hypothesis 2.1 Assume that for the parameters $\boldsymbol{\beta}=\left(\beta_1, \cdots, \beta_n\right)^{\mathrm{T}}, \boldsymbol{p}=\left(p_1, \cdots, p_n\right)^{\mathrm{T}}$ of the system (5) with (3), there exist a real number $s \in \mathbb{R}$ and a vector $\boldsymbol{\alpha} \in \operatorname{Ker} \boldsymbol{B}$ such that

$ \left\{\begin{array}{l} s \boldsymbol{\beta}>0, \\ -s \boldsymbol{p}+\boldsymbol{\alpha} <0 . \end{array}\right. $

Theorem 2.1 Under Hypothesis 2.1, for any given initial value $\boldsymbol{y}(0) \in \mathbb{R}_{+}^n$, there is a unique solution $\boldsymbol{y}(t)$ to the system (5) with (3) on t≥0, and the solution $\boldsymbol{y}(t)$ will remain in $\mathbb{R}_{+}^n$, with probability 1.

Proof The equivalent It form of the system (5) with (3) is

$ \begin{aligned} &\mathrm{d} \boldsymbol{y}(t)=\boldsymbol{B}(\boldsymbol{y}) \nabla H(\boldsymbol{y})(\mathrm{d} t+c \mathrm{~d} W(t))+ \\ &\frac{c^2}{2} \frac{\partial}{\partial \boldsymbol{y}}(\boldsymbol{B}(\boldsymbol{y}) \nabla H(\boldsymbol{y})) \boldsymbol{B}(\boldsymbol{y}) \nabla H(\boldsymbol{y}) \mathrm{d} t . \end{aligned} $

(6)

Taking the concrete expressions of $\boldsymbol{B}(\boldsymbol{y})$ and $H(\boldsymbol{y})$ in (3) into account, it is not difficult to see that its coefficients are locally Lipschitz continuous, and then it has a unique local solution $\boldsymbol{y}(t)$ on [0, τ_e), where τ_e is the explosion time^[12-13]. Next, we show this solution is global, i.e., $\tau_e=+\infty$ almost surely. Choose an integer k₀>0 such that every component of y(0) belongs to [1/k₀, k₀]. For each integer k≥k₀, define the stopping time

$ \begin{gathered} \tau_k:=\inf \left\{t \in\left[0, \tau_e\right): y_i(t) \notin(1 / k, k)\right. \\ \text { for some } i=1, \cdots, n\} \end{gathered} $

on the probability space $(\varOmega, \mathscr{F}, P)$. We set $\emptyset=+\infty$, which corresponds to the case when for certain $k^*, y_i(t) \in\left(1 / k^*, k^*\right)$, for all $i \in\{1$, $\cdots, n\}$ and $t \in\left[0, \tau_e\right)$. This can only happen when $\tau_e=+\infty$, due to continuity and construction of $\boldsymbol{y}(t)$ on $\left[0, \tau_e\right)$^[13].

Clearly, $\left\{\tau_k, k=0, 1, 2, \cdots\right\}$ is an increasing random sequence. Set $\tau_{\infty}=\lim _{k \rightarrow \infty} \tau_k$, then $\tau_{\infty} \leqslant \tau_e$ a.s.. Hence, to show $\tau_e=+\infty$ a.s. and $\boldsymbol{y}(t) \in \mathbb{R}_{+}^n$ a.s. for all t≥0, we only need to prove τ_∞=+∞ a.s..

If this statement is not true, then there exist real numbers T>0 and ε∈(0, 1) such that

$ P\left(\left\{\tau_{\infty} \leqslant T\right\}\right)>\varepsilon, $

which implies

$ P\left(\left\{\tau_k \leqslant T\right\}\right)>\varepsilon \text { for all } k \geqslant k_0, $

since $\left\{\tau_k, k=0, 1, 2, \cdots\right\}$ is an increasing random sequence. By Hypothesis 2.1, let

$ \begin{gathered} \left(a_1, a_2, \cdots, a_n\right):=s \boldsymbol{\beta}>0, \\ \left(d_1, d_2, \cdots, d_n\right):=-s \boldsymbol{p}+\boldsymbol{\alpha} <0 . \end{gathered} $

Now, for y>0 we construct the function

$\begin{aligned} &G(\boldsymbol{y}):=\sum\limits_{j=1}^n a_j y_j+d_j \ln y_j+d_j-d_j \ln \left(-\frac{d_j}{a_j}\right), \\ &G_j\left(y_j\right):=a_j y_j+d_j \ln y_j+d_j-d_j \ln \left(-\frac{d_j}{a_j}\right), \\ &j=1, \cdots, n . \end{aligned} $

(7)

We see that $G_j\left(y_j\right)(j=1, \cdots, n)$ are convex functions and have minimum value 0 on (0, +∞), and $G_j\left(y_j\right) \rightarrow+\infty$ as $y_j \rightarrow 0$ or +∞. Therefore

$ \begin{aligned} G\left(\boldsymbol{y}\left(\tau_k\right)\right) &=\sum\limits_{j=1}^n G_j\left(y_j\left(\tau_k\right)\right) \geqslant \\ G_i\left(y_i\left(\tau_k\right)\right) & \geqslant G_i(1 / k) \wedge G_i(k), \end{aligned} $

where $y_i$ denotes the element of y that runs beyond (1/k, k) at the time τ_k. Since $\boldsymbol{y}(t)$ is positive on [0, τ_k∧T], according to Proposition 1.1,

$ \begin{gathered} G(\boldsymbol{y}(t)) \equiv s H(\boldsymbol{y}(t))+C(\boldsymbol{y}(t))+ \\ \sum\limits_{j=1}^n\left(d_j-d_j \ln \left(-\frac{d_j}{a_j}\right)\right) \end{gathered} $

is an invarint of the system on the time interval [0, τ_k∧T]. Then set $\varOmega_k:=\left\{\tau_k \leqslant T\right\}$, we have

$ \begin{aligned} G(\boldsymbol{y}(0))=& E\left[G\left(\boldsymbol{y}\left(\tau_k \wedge T\right)\right)\right] \\ & \geqslant E\left[1_{\varOmega_{\mathrm{k}}} \mathrm{G}\left(\boldsymbol{y}\left(\tau_k \wedge T\right)\right)\right] \\ &=E\left[1_{\varOmega_{\mathrm{k}}} \mathrm{G}\left(\boldsymbol{y}\left(\tau_k\right)\right)\right] \\ &>\varepsilon\left[G_i(1 / k) \wedge G_i(k)\right] . \end{aligned} $

Let k→∞ in the above inequality, we then draw the contradiction

$ G(\boldsymbol{y}(0))>+\infty . $

Thus it holds τ_∞=+∞ almost surely.

Remark 2.1 When the constant matrix $\boldsymbol{B}$ is non-singular, i.e., Ker$\boldsymbol{B}$ contains only 0, the vector α must be equal to 0, in this case, the conditions in Hypothesis 2.1 can be simply expressed as β>0, p>0 or β < 0, p < 0. Of course, this case only occurs when the system is of even dimension, since an odd-dimensional skew-symmetric matrix $\boldsymbol{B}$ must be singular.

Remark 2.2 Under Hypothesis 2.1, the solution y(t) is positive on [0, +∞) almost surely, then f(H(y(t)), C(y(t))) given in Proposition 1.1 is a class of invariants of the system (5)with (3) on [0, +∞).

3 Boundedness of the positive solution

Based on Theorem 2.1, we can further obtain the boundedness of the positive solution y(t).

Proposition 3.1 Under Hypothesis 2.1, for any given initial value $\boldsymbol{y}(0) \in \mathbb{R}_{+}^n$, the solution y(t) of the system (5) with (3) is almost surely bounded.

Proof According to Remark 2.2, almost surely, the constructed function

$ \begin{aligned} G(\boldsymbol{y}(t))=& \sum\limits_{j=1}^n a_j y_j(t)+d_j \ln \left(y_j(t)\right)+\\ & d_j-d_j \ln \left(-\frac{d_j}{a_j}\right) \end{aligned} $

in (7) is an invariant of the system (5) with (3) on [0, +∞), i.e., G(y(t))≡G(y(0)). Moreover, the function G_j(y_j) has the minimum value $G_j\left(-\frac{d_j}{a_j}\right)$ which is equal to 0 on (0, +∞). Therefore

$ \begin{gathered} G_i\left(y_i(t)\right)=G(\boldsymbol{y}(0))- \\ \sum\limits_{j \neq i}^n G_j\left(y_j(t)\right) \leqslant G(\boldsymbol{y}(0)), \text { i. e. }, \\ a_i y_i(t)+d_i \ln \left(y_i(t)\right)+d_i- \\ d_i \ln \left(-\frac{d_j}{a_j}\right) \leqslant G(\boldsymbol{y}(0)) . \end{gathered} $

Hence, $y_i(t)$ locates in the bounded compact set $\left\{y_i \in \mathbb{R}_{+}: a_i y_i+d_i \ln y_i \leqslant G(\boldsymbol{y}(0))-d_i+\right.$ $\left.d_i \ln \left(-d_i / a_i\right)\right\}$ for all i=1, …, n, almost surely. Further, the equation of the tangent line of $G_i$ at the point $\left(-\frac{2 d_i}{a_i}, G_i\left(-\frac{2 d_i}{a_i}\right)\right)$ is

$z=\frac{a_i}{2}\left(y_i+\frac{2 d_i}{a_i}\right)+G_i\left(-\frac{2 d_i}{a_i}\right) . $

Then by the convexity of $G_i$(y_i) we obtain

$ \begin{gathered} \frac{a_i}{2}\left(y_i(t)+\frac{2 d_i}{a_i}\right)+G_i\left(-\frac{2 d_i}{a_i}\right) \leqslant \\ G_i\left(y_i(t)\right) \leqslant G(\boldsymbol{y}(0)), \end{gathered} $

which implies

$ 0 <y_i(t) \leqslant \frac{2\left[G(\boldsymbol{y}(0))-d_i \ln 2\right]}{a_i} $

for all t≥0 and i=1, …, n.

4 Numerical validations

In this section, we simulate the solutions of two concrete models of the form (5) with (3), by the numerical method proposed in Ref. [8] for stochastic Poisson systems of the form (5), which was proved to be of root mean-square convergence order 1, and reads

$ \begin{aligned} &y_{n+1}=y_n+\boldsymbol{B}\left(\frac{y_n+y_{n+1}}{2}\right) \\ &\int_0^1 \nabla H\left(y_n+\tau\left(y_{n+1}-y_n\right)\right) \mathrm{d} \tau\left(h+c \Delta \hat{W}_n\right) \end{aligned} $

(8)

where h is the time step, $\Delta \hat{W}_n:=\sqrt{h} \zeta_h$ is the truncated Wiener process increment on the time interval $\left[t_n, t_{n+1}\right]$, and

$ \zeta_h= \begin{cases}\xi, & \text { if } \quad \xi \leqslant\left|A_h\right|, \\ A_h, & \text { if } \quad \xi>A_h, \\ -A_h, & \text { if } \quad \xi <-A_h, \end{cases} $

where $A_h=\sqrt{4|\ln (h)|}$, and $\xi$ is a random variable subject to standard normal distribution. In this way we observe the behavior of the solutions, and validate the theoretical results on the solutions.

4.1 A three-dimensional model

Consider the three-dimensional Lotka-Volterra system with Stratonovich white noise perturbation^[8]

$ \begin{gathered} \mathrm{d} \boldsymbol{y}(t)=\boldsymbol{B}(\boldsymbol{y}(t)) \nabla H(\boldsymbol{y}(t))(\mathrm{d} t+c \circ \mathrm{d} W(t)), \\ \boldsymbol{B}(\boldsymbol{y})=\left(\begin{array}{ccc} 0 & v y_1 y_2 & b v y_1 y_3 \\ -v y_1 y_2 & 0 & -y_2 y_3 \\ -b v y_1 y_3 & y_2 y_3 & 0 \end{array}\right), \\ H(y)=a b y_1+y_2+\gamma \ln y_2-a y_3-\mu \ln y_3 . \end{gathered} $

(9)

By simple calculation one can varify that $\boldsymbol{B}(\boldsymbol{y})=$ $\operatorname{diag}\left(y_1, y_2, y_3\right) \boldsymbol{B} \operatorname{diag}\left(y_1, y_2, y_3\right)$ with $\boldsymbol{B}=$ $\left(\begin{array}{ccc}0 & v & b v \\ -v & 0 & -1 \\ -b v & 1 & 0\end{array}\right)$, and $\operatorname{Ker} \boldsymbol{B}=\left\{s_2 \boldsymbol{\alpha}: s_2 \in \mathbb{R}\right\}$, where $\boldsymbol{\alpha}=(-1 / v, -b, 1)^{\mathrm{T}}$ is a basis. Then a Casimir function of this system $C(\boldsymbol{y})=-1 / v \ln y_1-$ $b \ln y_2+\ln y_3$, for $\boldsymbol{y} \in \mathbb{R}_{+}^3$. In this system, $\boldsymbol{\beta}=$ $(a b, 1, -a)^{\mathrm{T}}, \boldsymbol{p}=(0, -\gamma, \mu)^{\mathrm{T}}$, we use the numerical scheme (8) to simulate this system, with two different groups of parameters.

First set $a=-2, b=-1, v=-0.5, \gamma=1, \mu=$ 2, and choose $s_1=1, s_2=-3$, such that $s_1(a b, 1$, $-a)=(2, 1, 2)>0$ and $-s_1(0, -\gamma, \mu)+$ $s_2(-1 / v, -b, 1)=(-6, -2, -5) <0$. Then the Hypothesis 2.1 is satisfied and the function $G(\boldsymbol{y})$ $=2 y_1-6 \ln y_1+y_2-2 \ln y_2+2 y_3-5 \ln y_3-13+$ $\ln \frac{2278125}{8}$. In Fig. 1(a), we show one sample path of $y_1, y_2$, and $y_3$, respectively, it can be seen that the solution is non-explosive, positive and bounded.

Then we set another group of parameters $a=$ $-0.2, b=-1, v=0.5, \gamma=-1, \mu=2$, and choose $s_1=1$, $s_2=0.5$, such that $s_1(a b, 1, -a)=(0.2, 1, 0.2)>0$ and $-s_1(0, -\gamma, \mu)+s_2(-1 / v, -b, 1)=(-1, -0.5, -1.5) <0$. Hence these parameters also meet the Hypothesis 2.1, and the function $G(\boldsymbol{y})=0.2 y_1-\ln y_1+y_2-$ $0.5 \ln y_2+0.2 y_3-1.5 \ln y_3-13+\ln \frac{5 \sqrt{3375}}{4}$. The sample paths of $y_1, y_2$ and $y_3$ under this set of parameters are shown in Fig. 1(b), which are also non-explosive, positive and bounded.

In Fig. 1, we take the parameter $c=0.5$, and the step size $h=10^{-3}$, initial value $\boldsymbol{y}(0)=(1.0$, $1.9, 0.5)^{\mathrm{T}}$ for Fig. 1 (a) and $\boldsymbol{y}(0)=(1.0, 1.5$, $0.5)^{\mathrm{T}}$ for Fig. 1(b), respectively.

4.2 A two-dimensional model

We consider a prey-predator model^[1] with random pertubation

$ \left(\begin{array}{l} \mathrm{d} u \\ \mathrm{~d} v \end{array}\right)=\left(\begin{array}{l} u(v-2) \\ v(1-u) \end{array}\right)(\mathrm{d} t+c \circ \mathrm{d} W(t))= \\ \left(\begin{array}{cc} 0 & u v \\ -u v & 0 \end{array}\right) \nabla H(u, v)(\mathrm{d} t+c \circ \mathrm{d} W(t)), $

(10)

where $H(u, v)=u-\ln u+v-2 \ln v$. In this case, $ \boldsymbol{B}(u, v)=\operatorname{diag}(u, v) \boldsymbol{B} \operatorname{diag}(u, v)$ with $\boldsymbol{B}= \left(\begin{array}{cc}0 & 1 \\ -1 & 0\end{array}\right) $, Ker$\boldsymbol{B}$ contains only 0, and $\boldsymbol{\beta}=(1, 1)^{\mathrm{T}} >0, \boldsymbol{p}=(1, 2)^{\mathrm{T}}>0$. Then the parameters satisfy Hypothesis 2.1, and the function $G(u, v)=u-\ln u+v-2 \ln v-3+2 \ln 2$. We use the numerical scheme (8) with $\boldsymbol{y}=: (u, v)^{\mathrm{T}}$ to simulate the solution of the system (10). As shown in Fig. 2(a) and Fig. 2(b), the solution is non-explosive, positive and bounded.

Here we take c=0.5, and the initial value (u(0), v(0))=(1.5, 2.5), the step size h=10^-3.

5 Conclusion

We prove the almost sure existence (global non-explosion), uniqueness and positiveness of the solution of a class of stochastic Poisson systems, under certain hypothesis, via constructing a function G(y) which is a special class of invariants of the systems. Almost sure boundedness of the solution is also verified. Numerical simulations give support to the theoretical results.