Invariants for a Dynamical System with Strong Random Perturbations

Abstract

In this chapter we consider the invariant method for stochastic system with strong perturbations, and its application to many different tasks related to dynamical systems with invariants. This theory allows constructing the mathematical model (deterministic and stochastic) of actual process if it has invariant functions. These models have a kind of jump-diffusion equations system (stochastic differential Itô equations with a Wiener and a Poisson paths). We show that an invariant function (with probability 1) for stochastic dynamical system under strong perturbations exists. We consider a programmed control with Prob. 1 for stochastic dynamical systems – PSP1. We study the construction of stochastic models with invariant function based on deterministic model with invariant one and show the results of numerical simulation. The concept of a first integral for stochastic differential equation Itô introduce by V. Doobko, and the generalized Itô – Wentzell formula for jump-diffusion function proved us, play the key role for this research.

1. Introduction

Models for actual dynamical processes are based on some restrictions. These restrictions are represented as a conservation law.

The conservation law states that a particular measurable property of an isolated dynamical system does not change as the system evolves over time.

Actual dynamical systems are open, and they are subject to strong external disturbances that violate the laws of conservation for the given system.

Conventionally, deterministic dynamical systems have an invariant function. Doobko¹ V. in [1] proved that stochastic dynamical systems have an invariant function as well. For dynamical system which are described using a system of stochastic differential Itô equations, a first integral – or an invariant function, exists with probability 1 [2, 3, 4, 5, 6, 7, 8, 9, 10].

When we know only a conservation law for a dynamical system, and equations which describing this system are unknown, the invariant functions are a good tool for determination of these equations.

Our method differs for other (see, for example, [11]) preliminary in the fact that we construct a system of differential equation with the given first integral under arbitrary initial conditions. Besides, this algorithm is realized as software and it allows us to choose a set of functions for simulation. Moreover, we can construct both a system of stochastic differential equations and a system of deterministic ones.

The goal of this chapter is representation of modern approach to describe of dynamical systems having a set of invariant functions.

This chapter is structured as follows. Firstly, we show that the invariant functions for stochastic systems exist. Then, the generalized Itô – Wentzell formula is represented. It is a differentiated rule for Jump-diffusion function under variables which solves the Jump-diffusion equations system. This rule is basic for the necessary and sufficient conditions for the stochastic first integral (or invariant function with probability 1) for the Jump-diffusion equations system. The next step is the construction of the differential equations system using the given invariant functions. It can be applied for stochastic and nonstochastic cases. The concept of PCP1 (Programmed control with Prob. 1) for stochastic dynamical systems is introduced. Finally, we show an application of the stochastic invariant theory for a transit from deterministic model with invariant to the same stochastic model. Several examples of application of this theory are given and confirmed by results of numerical calculations.

2. Notation and preliminaries

Now we introduce the main concepts which we will use below.

Let $w\left(t\right)$, $t\in \left[0,\infty \right)$ be a Wiener process or a (standard) Brownian motion, i. e.

• $w\left(0\right)=0$,

• it has stationary, independent increments,

• for every $t>0$, $w\left(t\right)$ has a normal $\mathcal{N}\left(0,t\right)$ distribution,

• it has continuous sample paths,

• every trajectory of $w\left(t\right)$ is not differentiated for all $t\ge 0$.

A $\nu \left(t,A\right)$ is called a Poisson random measure or standard Poisson measure (PM) if it is non-negative integer random variable with the Poisson distribution $\nu \left(t,A\right)\sim \mathit{Poi}\left(t\mathrm{\Pi }\left(A\right)\right)$, and it has the properties of measure:

• $\nu \left(t,A\right)$ is a random variable for every $t\in \left[0,T\right]$, $A\in {\mathbb{R}}^{n^{\prime }}$,

• $\nu \left(t,A\right)\in \mathbb{N}\cup \left\{0\right\}$, $\nu \left(t,\varnothing \right)=0$,

• if $A\cap B=\varnothing$, then $\nu \left(t,A\cup B\right)=\nu \left(t,A\right)+\nu \left(t,B\right)$,

• $\mathbf{E}\left[\nu \left(t,\mathcal{A}\right)\right]=t\mathrm{\Pi }\left(A\right)$,

• if $\#A$ is a number of random events from set $A$ during $t$, then

$${\mathbf{P}}_t\left(\#A=k\right)=\frac{{\left(t\mathrm{\Pi }\left(A\right)\right)}^k}{k!}exp\left\{-t\mathrm{\Pi }\left(A\right)\right\}.$$

$\tilde{\nu }\left(t,A\right)=\nu \left(t,A\right)-\mathbf{E}\left[\nu \left(t,A\right)\right]$ is called a centered Poisson measure (CPM).

Let $\mathbf{w}\left(t\right)={\left(w_1\left(t\right),\dots ,w_m\left(t\right)\right)}^{\ast }$ be an $m$-dimensional Wiener process, such that the one-dimensional Wiener processes $w_k\left(t\right)$ for $k=1,\dots ,m$ is mutually independent.

Take a vector $\gamma \in \mathrm{\Theta }$ with values in ${\mathbb{R}}^{n^{\prime }}$. Denote by $\nu \left(\mathrm{\Delta }t,\mathrm{\Delta }\gamma \right)$ the PM on $\left[0,T\right]\times {\mathbb{R}}^{n^{\prime }}$ modeling independent random variables on disjoint intervals and sets. The Wiener processes $w_k\left(t\right)$, $k=1,\dots ,m$, and the Poisson measure $\nu \left(\left[0,T\right],A\right)$ defined on the specified space are ${\mathcal{F}}_t$-measurable and independent of one another.

Consider a random process $x\left(t\right)$ with values in ${\mathbf{R}}^n$, $n\ge 2$, defined by the Equation [12]:

$$\mathit{dx}\left(t\right)=A\left(t\right)\mathit{dt}+B\left(t\right)d\mathbf{w}\left(t\right)+{\int }_{R_{\gamma }}g\left(t,\gamma ,\mathbf{x}\right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),$$E1

where $A\left(t\right)={\left\{a_1\left(t\right),\dots ,a_n\left(t\right)\right\}}^{\ast }$, $B\left(t\right)=\left(b_{j,k}\left(t\right)\right)$ is $\left(n\times k\right)$ - matrix, and $g\left(t,\gamma \right)={\left\{g_1\left(t,\gamma \right),\dots ,g_n\left(t,\gamma \right)\right\}}^{\ast }\in {\mathbb{R}}^n$, and $\gamma \in {\mathbb{R}}^{n^{\prime }}≕R_{\gamma }$, while $\mathbf{w}\left(t\right)$ is an $m$-dimensional Wiener process. In general the coefficients $A\left(t\right)$, $B\left(t\right)$, and $g\left(t,\gamma \right)$ are random functions depending also on $\mathbf{x}\left(t\right)$. Since the restrictions on these coefficients relate explicitly only to the variables $t$ and $\gamma$, we use precisely this notation for the coefficients of (1) instead of writing $A\left(t,\mathbf{x}\left(t\right)\right)$, $\left(t,\mathbf{x}\left(t\right)\right)$, and $g\left(t,\mathbf{x}\left(t\right),\gamma \right)$.

A system (1) is the stochastic differential Itô equation with Wiener and Poisson perturbations, which named below as a Jump-diffusion Itô equations system (GSDES).

We will consider the dynamical system described using ordinary deterministic differential equations (ODE) system and ordinary stochastic differential Itô equations (SDE) system of different types, taking into account the fact that $\mathbf{x}\in {\mathbf{R}}^n$, $n\ge 2$.

3. An existence of an invariant function (with Prob.1) for stochastic dynamical system under strong perturbations

Consider the diffusion Itô equation in ${\mathbf{R}}^3$ with orthogonal random action with respect to the vector of the solution

$$d\mathbf{v}\left(t\right)=-\mu \mathbf{v}\left(t\right)\mathit{dt}+\frac{b}{\mid \mathbf{v}\left(t\right)\mid }\left[\mathbf{v}\left(t\right)\times d\mathbf{w}\left(t\right)\right],$$E2

where $\mathbf{v}\in {\mathbf{R}}^3$, $\mathbf{w}\in {\mathbf{R}}^3$, and $w_i\left(t\right)$, $i=\mathrm{1,2,3}$ are independent Wiener processes. This equation is a specific form of the Langevin equation.

V. Doobko in [1] showed that the system (2) have an invariant function called a first integral of this system:

$$u\left(t,\mathbf{v}\right)=exp\left\{2\mathit{\mu t}\right\}\left({\left|\mathbf{v}\left(0\right)\right|}^2-\frac{b^2}{\mu }\right).$$

This, in particular, implies that

$$\underset{t\to \infty }{lim}{\left|\mathbf{v}\left(t\right)\right|}^2=\frac{b^2}{\mu },$$

i.e. process $\mid \mathbf{v}\left(t\right)\mid$ is a nonrandom function and the random process $\mathbf{v}\left(t\right)$ itself is generated in a sphere of constant radius $\frac{b}{\sqrt{\mu }}$.

In [4, 5, 10] it is shown that invariant function exists for other stochastic equations of Langevin type. To obtain this result, it is necessary to use the Itô’s formula.

4. The generalized Itô – Wentzell formula for jump-diffusion function

The rules for constructing stochastic differentials, e.g., the change rule, are very important in the theory of stochastic random processes. These are Itô’s formula [13, 14] for the differential of a nonrandom function of a random process and the Itô – Wentzell² formula [15] enabling us to construct the differential of a function which per se is a solution to a stochastic equation. Many articles address the derivation of these formulas for various classes of processes by extending Itô’s formula and the Itô – Wentzell formula to a larger class of functions.

The next level is to obtain a new formula for the generalized Itô Equation [14] which involves Wiener and Poisson components. In 2002, V. Doobko presented [7] a generalization of stochastic differentials of random functions satisfying GSDES with CPM based on expressions for the kernels of integral invariants (only the ideas of a possible proof) were sketched in [7]. The result is called” the generalized Itô – Wentzell formula”.

In contrast to [7], the generalized Itô – Wentzell formula for the noncentered Poisson measure was represented in [9, 16, 17]. The proof [9] of the generalized Itô – Wentzell formula uses the method of stochastic integral invariants and equations for their kernels. In this case the requirement on the character of the Poisson distribution is only a general restriction, as the knowledge of its explicit form is unnecessary. Other proofs in [16, 17] are based on traditional stochastic analysis and the use of approximations to random functions related to stochastic differential equations by averaging their values at each point.

The generalized Itô – Wentzell formula relying on the kernels of integral invariants [9] requires stricter conditions on the coefficients of all equations under consideration: the existence of second derivatives. The reason is that the kernels of invariants for differential equations exist under certain restrictions on the coefficients.

Since the random function $F\left(t,\mathbf{x}\left(t\right)\right)$ has representation as stochastic diffusion Itô equation with jumps, we can use the generalized Itô – Wentzell formula, proved by us by several methods in accordance with different conditions for the equations coefficients. Now we consider only one case.

We will use the following notation: ${\mathcal{C}}_y^s$ is the space of functions having continuous derivatives of order $s$ with respect to $y$, ${\mathcal{C}}_0^s\left(y\right)$ is the space of bounded functions having bounded continuous derivatives of order $s$ with respect to $y$.

Theorem 1.1 (generalized Itô – Wentzell formula). Consider the real function $F\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,2}$, $\left(t,\mathbf{x}\right)\in \left[0,T\right]\times {\mathbb{R}}^n$ with generalized stochastic differential of the form

$$d_{t}F\left(t,\mathbf{x}\right)=Q\left(t,\mathbf{x}\right)\mathit{dt}+\sum _{k=1}^m{D}_k\left(t,\mathbf{x}\right)d{\mathrm{w}}_k\left(t\right)+{\int }_{R_{\gamma }}G\left(t,\mathbf{x},\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right)$$E3

whose coefficients satisfy the conditions:

$$Q\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,2},D_k\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,2},G\left(t,\mathbf{x},\gamma \right)\in {\mathcal{C}}_{t,x,\gamma }^{\mathrm{1,2,1}}.$$

If a random process $\mathbf{x}\left(t\right)$ obeys (1) and its coefficients satisfy the conditions

$$a_i\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,1},b_{\mathit{ij}}\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,2},g_i\left(t,\mathbf{x},\gamma \right)\in {\mathcal{C}}_{t,x,\gamma }^{\mathrm{1,2,1}}.$$E4

then the stochastic differential exists and

$$\begin{array}{c}{d}_{t}F\left(t,\mathbf{x}\left(t\right)\right)=Q\left(t,\mathbf{x}\left(t\right)\right)\mathit{dt}+\sum _{k=1}^m{D}_k\left(t,\mathbf{x}\left(t\right)\right)d{\mathrm{w}}_k+\\ +\left[\sum _{i=1}^n{a}_i\left(t\right)\frac{\partial F\left(t,\mathbf{x}\right)}{\partial x_i}{\left|{}_{\mathbf{x}=\mathbf{x}\left(t\right)}+\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^m{b}_{i,k}\left(t\right)b_{j,k}\left(t\right)\frac{{\partial }^{2}F\left(t,\mathbf{x}\right)}{\partial x_i\partial x_j}\right|}_{\mathbf{x}=\mathbf{x}\left(t\right)}+\right.\\ {\left.+\left.\sum _{i=1}^n{b}_{i,k}\left(t\right){\frac{\partial D_k\left(t,\mathbf{x}\right)}{\partial x_i}}_{\mathbf{x}=\mathbf{x}\left(t\right)}\right]\mathit{dt}+\sum _{i=1}^n\sum _{k=1}^m{b}_{i,k}\left(t\right)\frac{\partial F\left(t,\mathbf{x}\right)}{\partial x_i}\right|}_{\mathbf{x}=\mathbf{x}\left(t\right)}d{\mathrm{w}}_k+\\ +{\int }_{R_{\gamma }}\left[\left(F\right(t,\mathbf{x}\left(t\right)+g(t,\gamma \left)\right)-F(t,\mathbf{x}\left(t\right))\right]\nu \left(\mathit{dt},\mathit{d\gamma }\right)+\\ +{\int }_{R_{\gamma }}G\left(t,\mathbf{x}\left(t\right)+g\left(t,\gamma \right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right).\end{array}$$E5

By analogy with the terminology proposed earlier, let us call formula (5) “the generalized Itô – Wentzell formula for the GSDES with PM” (GIWF).

By analogy with the classical Itô and Itô – Wentzell formulas, the generalized Itô – Wentzell formula is promising for various applications. In particular, it helped to obtain equations for the first and stochastic first integrals of the stochastic Itô system [9], equations for the density of stochastic dynamical invariants, Kolmogorov equations for the density of transition probabilities of random processes described by the generalized stochastic Itô differential Equation [8], as well as the construction of program controls with probability 1 for stochastic systems [18, 19].

5. A first integral for GSDES

In the theory of ODE, there are constructed equations to find deterministic functions, first integrals which preserve a constant value with any solutions to the equation. The concept of a first integral plays an important role in theoretical mechanics, for example, to solve inverse problems of mechanics or in constructing controls of dynamical systems.

It turned out that the first integral exists in the theory of stochastic differential equations (SDE) as well. However, there appears an additional classification connected with different interpretations. This gives a first integral for a system of SDE (see [1]), a first direct integral, and a first inverse integral for a system of Itô SDE (see [20]).

Definition 1.1 [1, 3]. Let $\mathbf{x}\left(t\right)$ be an $n$-dimensional random process satisfying a system of Itô SDE

$${\left.{\mathit{dx}}_i\left(t\right)=a_i(t,\mathbf{x}\left(t\right))\mathit{dt}+\sum _{k=1}^m{b}_{\mathit{ik}}(t,\mathbf{x}\left(t\right))d{\mathrm{w}}_k\left(t\right),\mathbf{x}(t,\mathbf{x}\left(0\right))\right|}_{t=0}=\mathbf{x}\left(0\right),$$E6

whose coefficients satisfy the conditions of the existence and uniqueness of a solution [12]. A nonrandom function $u\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,2}$ is called a first integral of the system of SDE if it takes a constant value depending only on $\mathbf{x}\left(0\right)$ on any trajectory solution to (6) with probability 1:

$$u(t,\mathbf{x}\left(t,\mathbf{x}\left(0\right)\right)=u\left(0,\mathbf{x}\left(0\right)\right)\text{almost}\text{surely},$$

or, in other words, its stochastic differential is equal to zero: $d_{t}u\left(t,\mathbf{x}\left(t\right)\right)=0.$

Another important notion in the theory of deterministic dynamical systems is given by the notion of an integral invariant introduced by Poincaré [21].

As it turned out, there also exist integral invariants for stochastic dynamical systems [2, 3]. In [7] V. Doobko give the concept of a kernel (=density) of a stochastic integral invariant and, based on it, formulate the notion of a stochastic first integral and a first integral as a deterministic function for GSDES with the centered Poisson measure, which makes it possible to compose a list of first integrals for stochastic differential equations.

Consider a random process $\mathbf{x}\left(t\right)$, $\mathbf{x}\in {\mathbb{R}}^n$, which is a solution to GSDES

$$\begin{array}{c}{\mathit{dx}}_i\left(t\right)=a_i\left(t,\mathbf{x}\left(t\right)\right)\mathit{dt}+b_{\mathit{ik}}\left(t,\mathbf{x}\left(t\right)\right)d{\mathrm{w}}_k\left(t\right)+{\int }_{R_{\gamma }}{g}_i\left(t,\mathbf{x}\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\left.\mathbf{x}\left(t\right)=x(t,\mathbf{x}\left(0\right),\omega )\right|}_{t=0}=\mathbf{x}\left(0\right),i=1,\dots ,n,t\ge 0,\end{array}$$E7

whose coefficients (in general, random functions) satisfy the conditions of the existence and uniqueness of a solution [12] and the following smoothness conditions:

$$a_i\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,1},b_{\mathit{ij}}\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,2},g_i\left(t,\mathbf{x},\gamma \right)\in {\mathcal{C}}_{t,x,\gamma }^{\mathrm{1,2,1}}.$$E8

Suppose that $\rho \left(t,\mathbf{x},\omega \right)$ is a random function connected with any deterministic function $f\left(t,\mathbf{x}\right)\in \mathfrak{S}\subset {\mathcal{C}}_0^{1,2}\left(t,\mathbf{x}\right)$ by the relations

$${\int }_{{\mathbb{R}}^n}\rho \left(t,\mathbf{x},\omega \right)f\left(t,\mathbf{x}\right)d\widehat{\mathrm{\Gamma }}\left(x\right)={\int }_{{\mathbb{R}}^n}\rho \left(0,\mathbf{y}\right)f\left(t,\mathbf{x}\left(t,\mathbf{y}\right)\right)d\widehat{\mathrm{\Gamma }}\left(\mathbf{y}\right)$$E9

$${\int }_{{\mathbb{R}}^n}\rho \left(0,\mathbf{x}\right)d\widehat{\mathrm{\Gamma }}\left(\mathbf{x}\right)=1,$$E10

$$\underset{\mid \mathbf{x}\mid \to \infty }{lim}\rho \left(0,\mathbf{x},\omega \right)=\underset{\mid \mathbf{x}\mid \to \infty }{lim}\rho \left(0,\mathbf{x}\right)=0,d\widehat{\mathrm{\Gamma }}\left(\mathbf{x}\right)=\prod _{i=1}^n{\mathit{dx}}_i,$$E11

where $\mathbf{y}≔\mathbf{x}\left(0\right)$, and $\mathbf{x}\left(t,\mathbf{y}\right)$ is a solution to (7), and $\omega$ is a random event.

In the particular case when $f\left(t,\mathbf{x}\right)=1$, conditions (9) and (10) imply that

$${\int }_{{\mathbb{R}}^n}\rho \left(t,\mathbf{x},\omega \right)d\widehat{\mathrm{\Gamma }}\left(\mathbf{x}\right)={\int }_{{\mathbb{R}}^n}\rho \left(0,\mathbf{y}\right)d\widehat{\mathrm{\Gamma }}\left(\mathbf{y}\right)=1,$$E12

i.e., for the random function $\rho \left(t,\mathbf{x},\omega \right)$, there exists a nonrandom functional preserving a constant value:

$${\int }_{{\mathbb{R}}^n}\rho \left(t,\mathbf{x},\omega \right)d\widehat{\mathrm{\Gamma }}\left(\mathbf{x}\right)=1.$$E13

Then, with conditions (10) and (11), Eq. (9) can be regarded as a stochastic integral invariant, and the function (t, x) can be viewed as its density.

Definition 1.2 [3]. A nonnegative random function $\rho \left(t,\mathbf{x},\omega \right)$ is referred to as a stochastic kernel or the stochastic density of a stochastic integral invariant (of $n$th order) if conditions (9), (10), and (11) are held.

Note that a substantial difference which made it possible to consider the invariance of the random volume on the basis of a kernel of an integral operator in [3, 7], is that (9) contains a functional factor. Thus, the notion of a kernel of an integral invariant [3] for a system of ordinary differential equations can be regarded as a particular case by taking $f\left(t,\mathbf{x}\right)=1$ and excluding from (7) the randomness determined by the Wiener and Poisson processes.

Using the GIWF (5), we obtain equation for the stochastic kernel function [9].

$$\begin{array}{c}{d}_t\rho \left(t,\mathbf{x},\omega \right)=-\frac{\partial \rho \left(t,\mathbf{x},\omega \right)b_{ik}\left(t,\mathbf{x}\right)}{\partial x_i}{\mathit{dw}}_k\left(t\right)+(-\frac{\partial \left(\rho \left(t,\mathbf{x},\omega \right)a_i(t,\mathbf{x})\right)}{\partial x_i}+\\ +\frac{1}{2}\frac{{\partial }^2\left(\rho \left(t,\mathbf{x},\omega \right)b_{ik}(t,\mathbf{x}\left)b_{jk}\right(t,\mathbf{x})\right)}{\partial x_i\partial x_j})\mathit{dt}+\\ +{\int }_{R_{\gamma }}\left[\rho \left(t,\mathbf{x}-g(t,{\mathbf{x}}^{-1}\left(t,\mathbf{x},\gamma \right),\gamma ,\omega )\right)\cdot \mathcal{J}\left({\mathbf{x}}^{-1}\left(t,\mathbf{x},\gamma \right)\right)-\rho \left(t,\mathbf{x},\omega \right)\right]\nu \left(\mathit{dt},\mathit{d\gamma }\right),\end{array}$$E14

under restrictions

$$\begin{array}{c}{\left.\rho \left(t,\mathbf{x},\omega \right)\right|}_{t=0}=\rho \left(0,\mathbf{x},\omega \right)=\rho \left(0,\mathbf{x}\right)\in {\mathcal{C}}_0^2\left(\mathbf{x}\right),\\ \underset{\mid x\mid \to \infty }{lim}\rho \left(0,\mathbf{x},\omega \right)=\underset{\mid \mathbf{x}\mid \to \infty }{lim}\rho \left(0,\mathbf{x}\right)=0,\underset{\mid \mathbf{x}\mid \to \infty }{lim}\frac{\partial \rho \left(0,\mathbf{x},\omega \right)}{\partial x_i}=\underset{\mid \mathbf{x}\mid \to \infty }{lim}\frac{\partial \rho \left(0,\mathbf{x}\right)}{\partial x_i}=0.\end{array}$$

This result plays a major role in obtaining of equation for the stochastic first integral.

6. Necessary and sufficient conditions for the stochastic first integral

Lemma 1.1. If $\rho \left(t,\mathbf{x},\omega \right)$ is a stochastic kernel of an integral invariant of $n$ th order of a stochastic process $\mathbf{x}\left(t\right)$ starting from a point $\mathbf{x}\left(0\right)$ then, for every $t$, it satisfies the equality

$$\rho \left(t,\mathbf{x}\left(t,\mathbf{x}\left(0\right)\right),\omega \right)\mathcal{J}\left(t,\mathbf{x}\left(0\right),\omega \right)=\rho \left(0,\mathbf{x}\left(0\right)\right),$$

where $\mathcal{J}\left(t,\mathbf{x}\left(0\right),\omega \right)$ is the Jacobian of transition from $\mathbf{x}\left(t\right)$ to $\mathbf{x}\left(0\right)$.

Definition 1.3 A set of kernels of integral invariants of $n$th order is called complete if any other function that is the kernel of this integral invariant can be presented as a function of the elements of this set.

In [9] it is shown that a system of GSDE (7) whose coefficients satisfy the conditions in (8), has a complete set of kernels consisting of $\left(n+1\right)$ functions.

Suppose that ${\rho }_l\left(t,\mathbf{x},\omega \right)\ne 0$, $l=1,\dots ,m$, $m\le n+1$ are kernels of the integral invariant (9). Lemma 1.1 implies that, for any $l\ne n+1$, the ratio $\frac{{\rho }_l\left(t,\mathbf{x}\left(t,\mathbf{y}\right),\omega \right)}{{\rho }_{n+1}\left(t,\mathbf{x}\left(t,\mathbf{y}\right),\omega \right)}$ is a constant depending only on the initial condition $\mathbf{x}\left(0\right)=\mathbf{y}$ for every solution $\mathbf{x}\left(t\right)$ to the GSDE (7) because

$$\frac{{\rho }_s\left(t,\mathbf{x}\left(t,\mathbf{y}\right),\omega \right)}{{\rho }_{n+1}\left(t,\mathbf{x}\left(t,\mathbf{y}\right),\omega \right)}=\frac{{\rho }_s\left(0,\mathbf{y}\right)}{{\rho }_{n+1}\left(0,\mathbf{y}\right)}.$$E15

Since for some realization ${\omega }_1$ we have

$$u\left(t,\mathbf{x}\left(t,\mathbf{x}\left(0\right)\right)\right)\equiv \frac{{\rho }_l\left(t,\mathbf{x}\left(t,\mathbf{x}\left(0\right)\right),{\omega }_1\right)}{{\rho }_s\left(t,\mathbf{x}\left(t,\mathbf{x}\left(0\right)\right),{\omega }_1\right)}=\frac{{\rho }_l\left(0,\mathbf{x}\left(0\right)\right)}{{\rho }_s\left(0,\mathbf{x}\left(0\right)\right)}\equiv u\left(0,\mathbf{x}\left(0\right)\right),$$

and it means, that $d_{t}u\left(t,\mathbf{x}\left(t\right)\right)=0$.

Definition 1.4 A random function $u\left(t,\mathbf{x},\omega \right)$ defined on the same probability space as a solution to (7) is referred to as a stochastic first integral of the system (7) of Itô ˆ GSDE with NCM if the following condition holds with probability 1:

$$u\left(t,\mathbf{x}\left(t,\mathbf{x}\left(0\right),\omega \right)\right)=u\left(0,\mathbf{x}\left(0\right)\right)\text{almost}\text{surely}$$

for every solution $\mathbf{x}\left(t,\mathbf{x}\left(0\right),\omega \right)$ to (7).

For practical purposes, for example, to construct program controls for a dynamical system under strong random perturbations, the presence of a concrete realization is important, i.e., the parameter $\omega$ is absent in what follows. In this connection, we introduce one more notion.

Definition 1.5 A nonrandom function $u\left(t,\mathbf{x}\right)$ is called a first integral of the system of GSDE (7) if it preserves a constant value with probability 1 for every realization of a random process $\mathbf{x}\left(t\right)$ that is a solution to this system:

$$u\left(t,\mathbf{x}\left(t,\mathbf{x}\left(0\right)\right)\right)=u\left(0,\mathbf{x}\left(0\right)\right)\text{almost}\text{surely}.$$

Thus, a stochastic first integral includes all trajectories (or realizations) of the random process while the first integral is related to one realization.

Construct an equation for $u\left(t,\mathbf{x},\omega \right)$ using the relation

$$ln{u}_s\left(t,\mathbf{x},\omega \right)=ln{\rho }_s\left(t,\mathbf{x},\omega \right)-ln{\rho }_l\left(t,\mathbf{x},\omega \right),$$E16

as a result of assertion (15). Let us differentiate $ln\rho \left(t,\mathbf{x}\right)$ (omit $\omega$) using generalized Itô – Wentzell formula:

$$\begin{array}{c}{d}_{t}ln\rho \left(t,\mathbf{x}\right)=\frac{1}{\rho \left(t,\mathbf{x}\right)}{\tilde{d}}_t\rho \left(t,\mathbf{x}\right)-\frac{1}{2{\rho }^2\left(t,\mathbf{x}\right)}{\left(-\frac{\partial \left(\rho \left(t,\mathbf{x}\right)b_{ik}\left(t,\mathbf{x}\right)\right)}{\partial x_i}\right)}^2\mathit{dt}+\\ +{\int }_{R_{\gamma }}\left[ln\left\{{\rho }_s\left(t,\mathbf{x}-g(t,\mathbf{x}\left(t,\mathbf{y}\right),\gamma ),\gamma \right)\mathcal{J}\left({\mathbf{x}}^{-1}\left(t,\mathbf{x},\gamma \right)\right)\right\}-ln{\rho }_s\left(t,\mathbf{x}\right)\right]\nu \left(\mathit{dt},\mathit{d\gamma }\right),\end{array}$$E17

where ${\tilde{d}}_t\rho \left(t,\mathbf{x}\right)$ is the right side of Eq.(14) without the integral expression. Having written down the equations for $ln{\rho }_s\left(t,\mathbf{x}\right)$ and $ln{\rho }_l\left(t,\mathbf{x}\right)$, and taking into account this result and Eq.(16), we obtain:

$$\begin{array}{c}{d}_{t}u\left(t,\mathbf{x},\omega \right)=\left[-a_i\left(t,x\right)\frac{\partial u\left(t,x,\omega \right)}{\partial x_i}+\frac{1}{2}{b}_{ik}(t,x\left)b_{jk}\right(t,x)\frac{{\partial }^{2}u\left(t,x,\omega \right)}{\partial x_i\partial x_j}-\right.\\ \left.-b_{ik}\left(t,\mathbf{x}\right)\frac{\partial }{\partial x_i}\left(b_{jk}\left(t,\mathbf{x}\right)\frac{\partial u\left(t,\mathbf{x},\omega \right)}{\partial x_j}\right)\right]\mathit{dt}-b_{ik}\left(t,\mathbf{x}\right)\frac{\partial u\left(t,\mathbf{x},\omega \right)}{\partial x_i}d{\mathrm{w}}_k\left(t\right)+\\ +{\int }_{R_{\gamma }}\left[u\left(t,\mathbf{x}-g(t,{\mathbf{x}}^{-1}\left(t,\mathbf{x},\gamma \right),\gamma ),\omega \right)-u\left(t,\mathbf{x},\omega \right)\right]\nu \left(\mathit{dt},\mathit{d\gamma }\right),\end{array}$$E18

which means that a stochastic first integral $u\left(t,\mathbf{x},\omega \right)$ of the Itô generalized Eq. (7) is a solution to the GSDE (18).

For a first integral which is a nonrandom function of one realization, the differential is also defined by an equation of the form of (18).

Theorem 1.2 Let $\mathbf{x}\left(t\right)$ be a solution to the GSDES (7) with conditions (8). A nonrandom function $u\left(t,\mathbf{x}\right)\in {\mathcal{C}}_{t,x}^{1,2}$ is a first integral of system (7) if and only if it satisfies the conditions:

1. $\frac{\partial u\left(t,\mathbf{x}\right)}{\partial t}+\frac{\partial u\left(t,\mathbf{x}\right)}{\partial x_i}\left[a_i\left(t,\mathbf{x}\right)-\frac{1}{2}{b}_{jk}(t,\mathbf{x})\frac{\partial b_{ik}\left(t,\mathbf{x}\right)}{\partial x_j}\right]=0$,

2. $b_{ik}\left(t,\mathbf{x}\right)\frac{\partial u\left(t,\mathbf{x}\right)}{\partial x_i}=0$, for all $k=\overline{1,m}$,

3. $u\left(t,\mathbf{x}\right)-u\left(t,\mathbf{x}+g\left(t,\mathbf{x},\gamma \right)\right)=0$ for any $\gamma \in R_{\gamma }$ in the entire domain of definition of the process.

Theorem (6) allows us to obtain a method for construction of differential equations systems on the basis of the given set of invariant functions.

7. Construction of the differential equations system using the given invariant functions

The concept of a first integral for a system of stochastic differential equations plays a key role in our theory. In this section, we will use a set of first integrals for the construction of a system of differential equations.

Let us write Eq. (7) in matrix form:

$$\begin{array}{c}\mathit{dX}\left(t\right)=A\left(t,X\left(t\right)\right)\mathit{dt}+B\left(t,X\left(t\right)\right)\mathit{dw}\left(t\right)+{\int }_{R_{\gamma }}\mathrm{\Theta }\left(t,X\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right)\\ X\left(0\right)=x_0,t\ge 0.\end{array}$$E19

Theorem 1.3 [22]. Let $X\left(t\right)$ be a solution of the Eq. (19) and let a nonrandom function $s\left(t,x\right)$ be continuous together with its first-order partial derivatives with respect to all its variables. Assume the set $\left\{{\overrightarrow{e}}_o,{\overrightarrow{e}}_1,\dots ,{\overrightarrow{e}}_n\right\}$ defines an orthogonal basis in ${\mathbf{R}}_{+}\times {\mathbf{R}}^n$. If function $s\left(t,x\right)$ is a first integral for the system (19), then the coefficients of Eq. (19) and the function $s\left(t,x\right)$ together are related by the conditions:

1. Functions $B_k\left(t,x\right)={\sum }_{i=1}^n{b}_{\mathit{ik}}\left(t,x\right){\overrightarrow{e}}_i$$\left(k=\left\{1,\dots ,m\right\}\right)$, which determine columns of the matrix $B\left(t,x\right)$, belong to a set

$$B_k\left(t,x\right)\in \left\{\begin{array}{c}{q}_{\mathit{oo}}\left(t,x\right)\cdot det\left[\begin{array}{ccc}{\overrightarrow{e}}_1& \dots & {\overrightarrow{e}}_n\\ \frac{\partial s\left(t,x\right)}{\partial x_1}& \dots & \frac{\partial s\left(t,x\right)}{\partial x_n}\\ f_{31}& \dots & f_{3n}\\ \dots & \dots & \dots \\ f_{n1}& \dots & f_{\mathit{nn}}\end{array}\right]\end{array}\right\},$$E20

where $q_{\mathit{oo}}\left(t,x\right)$ is an arbitrary nonvanishing function,

2. Coefficient $A\left(t,x\right)$ belongs to a set of functions defined by

$$A\left(t,x\right)\in \left\{R\left(t,x\right)+\frac{1}{2}\sum _{k=1}^n\left[\frac{\partial B_k\left(t,x\right)}{\partial x}\right]\cdot B_k(t,x)\right\},$$E21

where a column matrix $R\left(t,x\right)$ with components $r_i\left(t,x\right)$, $i=\left\{1,\dots ,n\right\}$, is defined as follows:

$$C^{-1}\left(t,x\right)\cdot detH\left(t,x\right)={\overrightarrow{e}}_o+\sum _{i=1}^n{r}_i\left(t,x\right){\overrightarrow{e}}_i,$$

$C\left(t,x\right)$ is an algebraic adjunct of the element ${\overrightarrow{e}}_o$ of a matrix $H\left(t,x\right)$ and $detC\left(t,x\right)\ne 0$, a matrix $H\left(t,x\right)$ is defined as

$$H\left(t,x\right)=\left[\begin{array}{cccc}{\overrightarrow{e}}_o& {\overrightarrow{e}}_1& \dots & {\overrightarrow{e}}_n\\ \frac{\partial s\left(t,x\right)}{\partial t}& \frac{\partial s\left(t,x\right)}{\partial x_1}& \dots & \frac{\partial s\left(t,x\right)}{\partial x_n}\\ h_{30}& h_{31}& \dots & h_{3n}\\ \dots & \dots & \dots & \dots \\ h_{n+1,0}& h_{n+1,1}& \dots & h_{n+1,n}\end{array}\right],$$E22

and $\left[\frac{\partial B_k\left(t,x\right)}{\partial x}\right]$ is a Jacobi matrix for function $B_k\left(t,x\right)$,

3. Coefficient $\mathrm{\Theta }\left(t,X,\gamma \right)={\sum }_{i=1}^n{\gamma }_i\left(t,x,\gamma \right){\overrightarrow{e}}_i$, related to Poisson measure, is defined by the representation $\mathrm{\Theta }\left(t,x,\gamma \right)=y\left(t,x,\gamma \right)-x$, where $y\left(t,x,\gamma \right)$ is a solution of the differential equations system

$$\begin{array}{c}\frac{\partial y\left(\cdot ,\gamma \right)}{\partial \gamma }=det\left[\begin{array}{cccc}{\overrightarrow{e}}_1& {\overrightarrow{e}}_2& \cdots & {\overrightarrow{e}}_n\\ \frac{\partial s\left(t,y\left(\cdot ,\gamma \right)\right)}{\partial y_1}& \frac{\partial s\left(t,y\left(\cdot ,\gamma \right)\right)}{\partial y_2}& \cdots & \frac{\partial s\left(t,y\left(\cdot ,\gamma \right)\right)}{\partial y_n}\\ {\phi }_{31}\left(t,y\left(\cdot ,\gamma \right)\right)& {\phi }_{32}\left(t,y\left(\cdot ,\gamma \right)\right)& \cdots & {\phi }_{3n}\left(t,y\left(\cdot ,\gamma \right)\right)\\ \cdots & \cdots & \cdots & \cdots \\ {\phi }_{n1}\left(t,y\left(\cdot ,\gamma \right)\right)& {\phi }_{n2}\left(t,y\left(\cdot ,\gamma \right)\right)& \cdots & {\phi }_{\mathit{nn}}\left(t,y\left(\cdot ,\gamma \right)\right)\end{array}\right].\end{array}$$E23

This solution satisfies the initial condition: ${\left.y\left(t,x,\gamma \right)\right|}_{\gamma =0}=x$.

The arbitrary functions $f_{\mathit{ij}}=f_{\mathit{ij}}\left(t,x\right)$, $h_{\mathit{ij}}=h_{\mathit{ij}}\left(t,x\right)$, and ${\phi }_{\mathit{ij}}={\phi }_{\mathit{ij}}\left(t,y\left(\cdot ,\gamma \right)\right)$ are defined by the equalities $f_{\mathit{ij}}\left(t,x\right)=\frac{\partial f_i\left(t,x\right)}{\partial x_j}$, $h_{\mathit{ij}}\left(t,x\right)=\frac{\partial h_i\left(t,x\right)}{\partial x_j}$, and ${\phi }_{\mathit{ij}}\left(t,y\left(\cdot ,\gamma \right)\right)=\frac{\partial {\phi }_i\left(t,y\left(\cdot ,\gamma \right)\right)}{\partial y_j}$. Sets of functions $\left\{{\phi }_i\left(t,y\left(\cdot ,\gamma \right)\right)\right\}$ and the function $g\left(t,x\right)$ together form a class of independent functions.

Using this theorem, we can to construct SDE system of different types and ODE system. Choice of arbitrary functions allows us to construct a set of differential equations systems with the given invariant functions. Theorem (7) allows us to introduce a concept of Programmed control with probability 1 for stochastic dynamical system.

8. Programmed control with Prob. 1 for stochastic dynamical systems

Definition 1.6 [18, 19]. A PCP1 is called a control of stochastic system which allows the preservation with probability 1 of a constant value for the same function which depends on this systems position for time periods of any length $T$.

Let us consider the stochastic nonlinear jump of diffusion equations system:

$$\begin{array}{c}\mathit{dX}\left(t\right)=\left(P\left(t,X\left(t\right)\right)+Z\left(t,X\left(t\right)\right)\cdot u\left(t,X\left(t\right)\right)\right)\mathit{dt}+\left(B\left(t,X\left(t\right)\right)+K\left(t,X\left(t\right)\right)\right)\mathit{dw}\left(t\right)+\\ +{\int }_{R_{\gamma }}\left(L\left(t,X\left(t\right),\gamma \right)+\mathrm{\Lambda }\left(t,X\left(t\right)\right)\right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\end{array}$$E24

where $P\left(\cdot \right)$, $Z\left(\cdot \right)$ are given matrix functions and $B\left(\cdot \right)$, $L\left(\cdot \right)$ are the functions that may either be known or not. For such systems we construct a unit of programmed control $\left\{u\left(t,X\left(t\right)\right),K\left(t,X\left(t\right)\right),M\left(t,X\left(t\right)\right)\right\}$ which allows the system (24) to be on the given manifold $\left\{u(t,X\left(t\right))\right\}=\left\{u\left(0,x_0\right)\right\}$ with Prob. 1 (PCP1) for each $t\in \left[0,T\right]$, $T\le \infty$.

Suppose that the nonrandom function $s\left(t,X\left(t\right)\right)$ is the first integral for the same stochastic dynamical system. The PCP1 $\left\{u\left(t,X\left(t\right)\right),K\left(t,X\left(t\right)\right),M\left(t,X\left(t\right)\right)\right\}$ is the solution for the algebraic system of linear equations.

Theorem 1.4 Let a controlled dynamical system be subjected to Brownian perturbations and Poisson jumps. The unit of PCP1 $\left\{u\left(t,X\left(t\right)\right),K\left(t,X\left(t\right)\right),M\left(t,X\left(t\right)\right)\right\}$, allowing this system to remain with probability 1 on the dynamically structured integral mfd $s\left(t,X\left(t,x_o\right),\omega \right)=s\left(0,x_o\right)$, is a solution of the linear equations system (with respect to functions $\mathbf{u}\left(t,\mathbf{x}\left(t\right)\right)$), $K\left(t,X\left(t\right)\right)$, $M\left(t,X\left(t\right)\right)$ which consists of Eq. (19) and Eq. (24). The coefficients of the Eq. (19), (and the coefficients of the Eq. (24) respectively) are determined by the theorem 7. The response to the random action is defined completely.

We show how the stochastic invariants theory can be applied to solve different tasks.

9. Stochastic models with invariant function which are based on deterministic model with invariant one

In this section we consider a few examples for application of the theory above to modeling actual random processes with invariants [23]. Firstly, we consider an example of construction of a differential equation system with the given invariant. Secondly, we study a general scheme for the PCP1 determination. And finally, we show the possibility of construction of stochastic analogues for classical models described by a differential equations system with an invariant function. The suggested method of stochastization is based on both the concept of the first integral for a stochastic differentialItô equations system (SDE) and the theorem for construction of the SDE system using its first integral.

9.1 Construction of a differential equations system

It is necessary to construct a differential equations system for $\mathbf{X}\in {\mathbf{R}}^3$, $t\ge 0$ such that the equality

$$X\left(t\right)-Y_1^2\left(t\right)+Y_2\left(t\right)+e^t=0$$E25

is satisfied with Prob.1. The equality (25) means that the differential equations system has a first integral $s\left(t,X\left(t\right),Y_1\left(t\right),Y_2\left(t\right)\right)=X\left(t\right)-Y_1^2\left(t\right)+Y_2\left(t\right)+e^t$ with initial condition ${\left(\mathrm{0,1,0}\right)}^{\ast }$:

$$s\left(t,X\left(t\right),Y_1\left(t\right),Y_2\left(t\right)\right)\equiv X\left(t\right)-Y_1^2\left(t\right)+Y_2\left(t\right)+e^t=s\left(0,x\left(0\right),Y_1\left(0\right),Y_2\left(0\right)\right).$$

We have

$$B_k\left(\cdot \right)=q_{00}\left(\cdot \right)det\left[\begin{array}{ccc}{\overrightarrow{e}}_1& {\overrightarrow{e}}_2& {\overrightarrow{e}}_3\\ 1& -2y_1& 1\\ f_1\left(\cdot \right)& f_2\left(\cdot \right)& f_3\left(\cdot \right)\end{array}\right]=q_{00}\left(\cdot \right)\left[\begin{array}{c}-2y_1{f}_3\left(\cdot \right)-f_2\left(\cdot \right)\\ -f_3\left(\cdot \right)+f_1\left(\cdot \right)\\ f_2\left(\cdot \right)+2y_1{f}_1\left(\cdot \right)\end{array}\right]=\left[\begin{array}{c}{b}_{1k}\left(\cdot \right)\\ b_{2k}\left(\cdot \right)\\ b_{3k}\left(\cdot \right)\end{array}\right].$$

Therefore,

$$\begin{array}{c}\left(B_k\left(\cdot \right),{\nabla }_z\right)B_k\left(\cdot \right)=q_{00}^2\left(\cdot \right)\left[\begin{array}{ccc}\frac{\partial \left(-2Y_1{f}_3\left(\cdot \right)-f_2\left(\cdot \right)\right)}{\partial x}& \frac{\partial \left(-2Y_1{f}_3\left(\cdot \right)-f_2\left(\cdot \right)\right)}{\partial Y_1}& \frac{\partial \left(-2Y_1{f}_3\left(\cdot \right)-f_2\left(\cdot \right)\right)}{\partial Y_2}\\ \frac{\partial \left(-f_3\left(\cdot \right)+f_1\left(\cdot \right)\right)}{\partial x}& \frac{\partial \left(-f_3\left(\cdot \right)+f_1\left(\cdot \right)\right)}{\partial Y_1}& \frac{\partial \left(-f_3\left(\cdot \right)+f_1\left(\cdot \right)\right)}{\partial Y_2}\\ \frac{\partial \left(f_2\left(\cdot \right)+2Y_1{f}_1\left(\cdot \right)\right)}{\partial x}& \frac{\partial \left(f_2\left(\cdot \right)+2Y_1{f}_1\left(\cdot \right)\right)}{\partial Y_1}& \frac{\partial \left(f_2\left(\cdot \right)+2Y_1{f}_1\left(\cdot \right)\right)}{\partial Y_2}\end{array}\right]\times \\ \times \left[\begin{array}{c}-2Y_1\left(t\right)f_3\left(\cdot \right)-f_2\left(\cdot \right)\\ -f_3\left(\cdot \right)+f_1\left(\cdot \right)\\ f_2\left(\cdot \right)+2Y_1\left(t\right)f_1\left(\cdot \right)\end{array}\right]=\left[\begin{array}{c}{p}_1\left(\cdot \right)\\ p_2\left(\cdot \right)\\ p_3\left(\cdot \right)\end{array}\right]\end{array}.$$

Thus the new drift coefficients are

$$\begin{array}{c}{A}_1\left(\cdot \right)=\frac{e^t\left(h_2{f}_3-f_2{h}_3\right)+2Y_1\left(h_0{f}_3-f_0{h}_3\right)+h_0{f}_2-f_0{h}_2}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{p_1}{2},\\ A_2\left(\cdot \right)=\frac{e^t\left(h_3{f}_1-f_3{h}_1\right)+h_0{f}_3-f_0{h}_3+h_0{f}_1-f_0{h}_1}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{p_2}{2},\\ A_3\left(\cdot \right)=\frac{e^t\left(h_2{f}_1-f_2{h}_1\right)-2Y_1\left(h_0{f}_1-f_0{h}_1\right)+h_0{f}_2-f_0{h}_2}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{p_3}{2}.\end{array}$$E26

According to term 3 of Theorem 1.3, we will determine a coefficient for Poisson measure. Now we rename variables: $Z\equiv {\left(Z_1,Z_2,Z_3\right)}^{\ast }≔{\left(X,Y_1,Y_2\right)}^{\ast }$. Then, we have:

$$\begin{array}{c}u\left(t,Z\right)=Z_1-Z_2^2+Z_3+e^t,\\ u\left(t,Z\right)-u\left(t,Z+g\left(t,Z,\gamma \right)\right)\equiv u\left(t,Z\right)-u\left(t,V\right)=0,\\ V=Z+g\left(t,Z,\gamma \right),\\ g\left(t,Z,\gamma \right)=V\left(t,Z,\gamma \right)-Z,\end{array}$$

where a function $V\left(t,Z,\gamma \right)$ solves the a differential equations system:

$$\frac{\partial V\left(\cdot ,\gamma \right)}{\partial \gamma }=det\left[\begin{array}{ccc}{\overrightarrow{e}}_1& {\overrightarrow{e}}_2& {\overrightarrow{e}}_3\\ 1& -2Z_2& 1\\ {\phi }_1\left(\cdot ,\gamma \right)& {\phi }_2\left(\cdot ,\gamma \right)& {\phi }_3\left(\cdot ,\gamma \right)\end{array}\right]=\left[\begin{array}{c}-2Z_2{\phi }_3\left(\cdot ,\gamma \right)-{\phi }_2\left(\cdot ,\gamma \right)\\ {\phi }_1\left(\cdot ,\gamma \right)-{\phi }_3\left(\cdot ,\gamma \right)\\ {\phi }_2\left(\cdot ,\gamma \right)+2Z_2{\phi }_1\left(\cdot ,\gamma \right),\end{array}\right],$$

and satisfies the initial condition $V\left(\cdot ,0\right)=Z$. Then, we determine functions $g_1\left(\cdot ,\gamma \right),g_2\left(\cdot ,\gamma \right),g_3\left(\cdot ,\gamma \right)$.

Assume, that ${\phi }_1\left(\cdot ,\gamma \right)=\frac{1}{\gamma +1}$, ${\phi }_2\left(\cdot ,\gamma \right)=2\gamma$, ${\phi }_3\left(\cdot ,\gamma \right)=1$. Then, we get:

$$\begin{array}{c}{g}_1\left(t,X\left(t\right),Y\left(t,\right),\gamma \right)=-2Y_1\left(t\right)\gamma -{\gamma }^2-x\left(t\right),\\ g_2\left(t,X\left(t\right),Y\left(t,\right),\gamma \right)=ln\mid \gamma +1\mid -\gamma +1-Y_1\left(t\right),\\ g_3\left(t,X\left(t\right),Y\left(t,\right),\gamma \right)=-{\gamma }^2+2Y_{1}ln\mid \gamma +1\mid -Y_2.\end{array}$$

Finally, we have constructed three variants of differential equations system:

1. deterministic differential equations system:

$$\left\{\begin{array}{c}\frac{\mathit{dX}\left(t\right)}{\mathit{dt}}=\frac{e^t\left(h_2{f}_3-f_2{h}_3\right)+2Y_1\left(h_0{f}_3-f_0{h}_3\right)+h_0{f}_2-f_0{h}_2}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}\\ \frac{{\mathit{dY}}_1\left(t\right)}{\mathit{dt}}=\frac{e^t\left(h_3{f}_1-f_3{h}_1\right)+h_0{f}_3-f_0{h}_3+h_0{f}_1-f_0{h}_1}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}\\ \frac{{\mathit{dY}}_2\left(t\right)}{\mathit{dt}}=\frac{e^t\left(h_2{f}_1-f_2{h}_1\right)-2Y_1\left(h_0{f}_1-f_0{h}_1\right)+h_0{f}_2-f_0{h}_2}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}.\end{array}\right.$$

2. stochastic differential equations system (Itó diffusion equations):

$$\left\{\begin{array}{c}\mathit{dX}\left(t\right)=\left[\frac{e^t\left(h_2{f}_3-f_2{h}_3\right)+2Y_1\left(h_0{f}_3-f_0{h}_3\right)+h_0{f}_2-f_0{h}_2}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{1}{2}\left(-2Y_1\left(t\right)f_3\left(\cdot \right)-f_2\left(\cdot \right)\right)\right]\mathit{dt}+\\ +\left(-2y_1{f}_3\left(\cdot \right)-f_2\left(\cdot \right)\right)d{\mathrm{w}}_1\left(t\right)\\ {\mathit{dY}}_1\left(t\right)=\left[\frac{e^t\left(h_3{f}_1-f_3{h}_1\right)+h_0{f}_3-f_0{h}_3+h_0{f}_1-f_0{h}_1}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{1}{2}\left(-f_3\left(\cdot \right)+f_1\left(\cdot \right)\right)\right]\mathit{dt}+\\ +\left(-f_3\left(\cdot \right)+f_1\left(\cdot \right)\right)d{\mathrm{w}}_2\left(t\right)\\ {\mathit{dY}}_2\left(t\right)=\left[\frac{e^t\left(h_2{f}_1-f_2{h}_1\right)-2Y_1\left(h_0{f}_1-f_0{h}_1\right)+h_0{f}_2-f_0{h}_2}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{1}{2}\left(f_2\left(\cdot \right)+2Y_1\left(t\right)f_1\left(\cdot \right)\right)\right]\mathit{dt}+\\ +\left(f_2\left(\cdot \right)+2y_1{f}_1\left(\cdot \right)\right)d{\mathrm{w}}_2\left(t\right)\\ .\end{array}\right.$$

3. stochastic differential equations system (jump-diffusion Itó equations):

$$\left\{\begin{array}{c}\mathit{dX}\left(t\right)=\left[\frac{e^t\left(h_2{f}_3-f_2{h}_3\right)+2Y_1\left(h_0{f}_3-f_0{h}_3\right)+h_0{f}_2-f_0{h}_2}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{1}{2}\left(-2Y_1\left(t\right)f_3\left(\cdot \right)-f_2\left(\cdot \right)\right)\right]\mathit{dt}+\\ +\left(-2y_1{f}_3\left(\cdot \right)-f_2\left(\cdot \right)\right)d{\mathrm{w}}_1\left(t\right)+{\int }_{R_{\gamma }}\left[-2Y_1\left(t\right)\gamma -{\gamma }^2-x\left(t\right)\right]\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\mathit{dY}}_1\left(t\right)=\left[\frac{e^t\left(h_3{f}_1-f_3{h}_1\right)+h_0{f}_3-f_0{h}_3+h_0{f}_1-f_0{h}_1}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{1}{2}\left(-f_3\left(\cdot \right)+f_1\left(\cdot \right)\right)\right]\mathit{dt}+\\ +\left(-f_3\left(\cdot \right)+f_1\left(\cdot \right)\right)d{\mathrm{w}}_2\left(t\right)+{\int }_{R_{\gamma }}\left[ln|\gamma +1|-\gamma +1-Y_1\left(t\right)\right]\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\mathit{dY}}_2\left(t\right)=\left[\frac{e^t\left(h_2{f}_1-f_2{h}_1\right)-2Y_1\left(h_0{f}_1-f_0{h}_1\right)+h_0{f}_2-f_0{h}_2}{f_2{h}_3-h_2{f}_3+f_1{h}_2-h_1{g}_2-2Y_1\left(h_1{f}_3-f_1{h}_3\right)}+\frac{1}{2}\left(f_2\left(\cdot \right)+2Y_1\left(t\right)f_1\left(\cdot \right)\right)\right]\mathit{dt}+\\ +\left(f_2\left(\cdot \right)+2y_1{f}_1\left(\cdot \right)\right)d{\mathrm{w}}_2\left(t\right)+{\int }_{R_{\gamma }}\left[-{\gamma }^2+2Y_{1}ln|\gamma +1|-Y_2\right]\nu \left(\mathit{dt},\mathit{d\gamma }\right).\end{array}\right.$$E27

We choose the functions $q_{00}\left(\cdot \right)$, $f_i\left(\cdot \right)$ and $h_i\left(\cdot \right)$, $i=\mathrm{1,2,3}$, in accordance with the restriction of the task and taking into account the utility for modeling.

9.2 Transit from deterministic model with invariant to the same stochastic model

Now we describe a general scheme for application of the theory above.

The suggested method of stochastization is based on both the concept of the first integral for a stochastic differentialItô equations system (SDE) and the theorem for construction of the SDE system using its first integral.

Let us consider a classical model

$$\left\{\begin{array}{l}{\mathit{dy}}_1\left(t\right)=F_1\left(t,\mathbf{y}\left(t\right)\right)\mathit{dt},\\ {\mathit{dy}}_2\left(t\right)=F_2\left(t,\mathbf{y}\left(t\right)\right)\mathit{dt},\\ {\mathit{dy}}_3\left(t\right)=F_3\left(t,\mathbf{y}\left(t\right)\right)\mathit{dt},\end{array}\right.$$E28

with an invariant $u\left(t,\mathbf{y}\right)$.

Then we construct the GSDE system, taking into account the equality $u\left(t,\mathbf{x}\left(t\right)\right)=u\left(0,\mathbf{x}\left(0\right)\right)=C$:

$$\left\{\begin{array}{l}{\mathit{dx}}_1\left(t\right)=a_1\left(t,\mathbf{x}\left(t\right)\right)\mathit{dt}+b_1\left(t,\mathbf{x}\left(t\right)\right)d\mathrm{w}\left(t\right)+\int g_1\left(t,\mathbf{x}\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\mathit{dx}}_2\left(t\right)=a_2\left(t,\mathbf{x}\left(t\right)\right)\mathit{dt}+b_2\left(t,\mathbf{x}\left(t\right)\right)d\mathrm{w}\left(t\right)+\int g_2\left(t,\mathbf{x}\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\mathit{dx}}_3\left(t\right)=a_3\left(t,\mathbf{x}\left(t\right)\right)\mathit{dt}+b_3\left(t,\mathbf{x}\left(t\right)\right)d\mathrm{w}\left(t\right)+\int g_3\left(t,\mathbf{x}\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right).\end{array}\right.$$E29

Hence, the stochastic model has a representation

$$\left\{\begin{array}{l}{\mathit{dy}}_1\left(t\right)=a_1\left(t,\mathbf{y}\left(t\right)\right)\mathit{dt}+b_1\left(t,\mathbf{y}\left(t\right)\right)d\mathrm{w}\left(t\right)+\int g_1\left(t,\mathbf{y}\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\mathit{dy}}_2\left(t\right)=a_2\left(t,\mathbf{y}\left(t\right)\right)\mathit{dt}+b_2\left(t,\mathbf{y}\left(t\right)\right)d\mathrm{w}\left(t\right)+\int g_2\left(t,\mathbf{y}\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\mathit{dy}}_3\left(t\right)=a_3\left(t,\mathbf{y}\left(t\right)\right)\mathit{dt}+b_3\left(t,\mathbf{y}\left(t\right)\right)d\mathrm{w}\left(t\right)+\int g_3\left(t,\mathbf{y}\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ \mathbf{y}\left(0\right)={\mathbf{y}}_0.\end{array}\right.$$E30

Further, we determine complementary function which is unit of control functions for PCP1:

$$\left\{\begin{array}{l}{s}_1\left(t,\mathbf{y}\left(t\right)\right)=a_1\left(t,\mathbf{y}\left(t\right)\right)-F_1\left(t,\mathbf{y}\left(t\right)\right),\\ s_2\left(t,\mathbf{y}\left(t\right)\right)=a_2\left(t,\mathbf{y}\left(t\right)\right)-F_2\left(t,\mathbf{y}\left(t\right)\right),\\ s_3\left(t,\mathbf{y}\left(t\right)\right)=a_3\left(t,\mathbf{y}\left(t\right)\right)-F_3\left(t,\mathbf{y}\left(t\right)\right).\end{array}\right.$$E31

Finally, we have constructed stochastic analogue for classical model described by a differential equations system and having an invariant function.

9.3 The SIR (susceptible-infected-recovered) model

The SIR is a simple mathematical model of epidemic [24], which divides the (fixed) population of $N$ individuals into three” compartments” which may vary as a function of time $t$.

$S\left(t\right)$ are those susceptible but not yet infected with the disease,

$I\left(t\right)$ is the number of infectious individuals,

$R\left(t\right)$ are those individuals who have recovered from the disease and now have immunity to it,

the parameter $\lambda$ describes the effective contact rate of the disease,

the parameter $\mu$ is the mean recovery rate.

The SIR model describes the change in the population of each of these compartments in terms of two parameters:

$$\left\{\begin{array}{l}\frac{\mathit{dS}\left(t\right)}{\mathit{dt}}=-\lambda \frac{S\left(t\right)I\left(t\right)}{N},\\ \frac{\mathit{dI}\left(t\right)}{\mathit{dt}}=\lambda \frac{S\left(t\right)I\left(t\right)}{N}-\mathit{\mu I}\left(t\right),\\ \frac{\mathit{dR}\left(t\right)}{\mathit{dt}}=\mathit{\mu I}\left(t\right),\end{array}\right.$$E32

and its restrictsion is

$$S\left(t\right)+I\left(t\right)+R\left(t\right)=N.$$E33

Let the model with strong perturbation be

$$\left\{\begin{array}{l}\mathit{dS}\left(t\right)=\left(-\lambda \frac{S\left(t\right)I\left(t\right)}{N}+s_1\left(t,S\left(t\right),I\left(t\right),R\left(t\right)\right)\right)\mathit{dt}+\\ +b_1\left(t,S\left(t\right),I\left(t\right),R\left(t\right)\right)d\mathrm{w}\left(t\right)+{\int }_{R_{\gamma }}{g}_1\left(t,S\left(t\right),I\left(t\right),R\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ \mathit{dI}\left(t\right)=\left(\lambda \frac{S\left(t\right)I\left(t\right)}{N}-\mathit{\mu I}\left(t\right)+s_2\left(t,S\left(t\right),I\left(t\right),R\left(t\right)\right)\right)\mathit{dt}+\\ +b_1\left(t,S\left(t\right),I\left(t\right),R\left(t\right)\right)d\mathrm{w}\left(t\right)+{\int }_{R_{\gamma }}{g}_2\left(t,S\left(t\right),I\left(t\right),R\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ \mathit{dR}\left(t\right)=\left(\mathit{\mu S}\left(t\right)+s_3\left(t,S\left(t\right),I\left(t\right),R\left(t\right)\right)\right)\mathit{dt}+\\ +b_1\left(t,S\left(t\right),I\left(t\right),R\left(t\right)\right)d\mathrm{w}\left(t\right)+{\int }_{R_{\gamma }}{g}_3\left(t,S\left(t\right),I\left(t\right),R\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\end{array}\right.$$E34

and

$$u\left(t,S\left(t\right),I\left(t\right),R\left(t\right)\right)=S\left(t\right)+I\left(t\right)+R\left(t\right)-N\equiv 0.$$E35

Suppose that the function $u\left(t,x,y,z\right)=x+y+z-N$ is a first integral, $v\left(t,x,y,z\right)=2e^{-t}+x$ and $h\left(t,x,y,z\right)=y$ are complementary functions, and $q\left(t,x,y,z\right)=x$ is arbitrary function. The initial condition is: $x\left(0\right)=1$, $y\left(0\right)=0$, $z\left(0\right)=0$. Then constructed differential equations system has the form

$$\begin{array}{c}\left[\begin{array}{c}\mathit{dx}\left(t\right)\\ \mathit{dy}\left(t\right)\\ \mathit{dz}\left(t\right)\end{array}\right]=\left[\begin{array}{l}2e^{-t}\\ 0\\ -2e^{-t}\end{array}\right]\mathit{dt}+\left[\begin{array}{cc}0& \\ x\left(t\right)& \\ -x\left(t\right)& \end{array}\right]d\mathrm{w}\left(t\right)+{\int }_{R_{\gamma }}\left[\begin{array}{c}0\\ x\left(t\right)\gamma )\\ -x\left(t\right)\gamma \end{array}\right]\nu \left(\mathit{dt},\mathit{d\gamma }\right).\end{array}$$E36

Let us simulate a numerical solution of Eg.(36), where $N=1$ (for example). Figure 1 shows simulation for system without jumps, the Figure 2 shows the processes with jumps.

Figure 1.

Numerical solution for Eq.(36) without jumps.

Figure 2.

Numerical solution for Eq.(36) with jumps.

In such a way we could use the system of differential equations

$$\left\{\begin{array}{l}{\mathit{dy}}_1\left(t\right)=2e^{-t}\mathit{dt},\\ {\mathit{dy}}_2\left(t\right)=y_1\left(t\right)d\mathrm{w}\left(t\right)+\int \gamma y_1\left(t\right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\mathit{dy}}_3\left(t\right)=-2e^{-t}\mathit{dt}+y_1\left(t\right)d\mathrm{w}\left(t\right)-\int \gamma y_1\left(t\right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ y\left(0\right)=y_0,\end{array}\right.$$E37

as initial step for construction of stochastic SIR-model. A good choice of complementary functions $v\left(t,x,y,z\right)$ and $h\left(t,x,y,z\right)$ allows us to obtain such coefficients that ensure that the solution $\left\{x\left(t\right),y\left(t\right),z\left(t\right)\right\}$ of the differential equations system satisfy some reasonable limitations.

9.4 The predator–prey model

The Lotka - Volterra equations or the predator–prey equations used to describe the dynamics of biological systems in which two species interact, one as a predator and the other as prey.

The Lotka - Volterra model makes a number of assumptions, not necessarily realizable in nature, about the environment and evolution of the predator and prey populations:

• The prey population finds ample food at all times.

• The food supply of the predator population depends entirely on the size of the prey population.

• The rate of change of population is proportional to its size.

• During the process, the environment does not change in favor of one species, and genetic adaptation is inconsequential.

• Predators have limitless appetite.

Let us note: $N_1\left(t\right)$ is the number of prey, and $N_2\left(t\right)$ is the number of some predator, ${\epsilon }_1$, ${\epsilon }_2$, ${\eta }_1$ and ${\eta }_2$ are positive real parameters describing the interaction of the two species.

The populations change through time according to the pair of equations:

$$\left\{\begin{array}{l}{\mathit{dN}}_1\left(t\right)=N_1\left(t\right)\left({\epsilon }_1-{\eta }_1{N}_2\left(t\right)\right)\mathit{dt},\\ {\mathit{dN}}_2\left(t\right)=-N_2\left(t\right)\left({\epsilon }_2-{\eta }_2{N}_1\left(t\right)\right)\mathit{dt}.\end{array}\right.$$E38

Eq. (38) has the invariant function

$$N_1^{-{\epsilon }_2}\left(t\right)e^{{\eta }_2{N}_1\left(t\right)}={\mathit{CN}}_2^{{\epsilon }_1}\left(t\right)e^{-{\eta }_1{N}_2\left(t\right)},$$E39

where $C=\mathit{\text{const}}$.

We can introduce the stochastic model as a form

$$\left\{\begin{array}{l}{\mathit{dx}}_1\left(t\right)=\left({\epsilon }_1{x}_1\left(t\right)-{\eta }_1{x}_1\left(t\right)x_2\left(t\right)+s_1\left(t,x_1\left(t\right),x_2\left(t\right)\right)\right)\mathit{dt}+\\ +b_1\left(t,x_1\left(t\right),x_2\left(t\right)\right)d{\mathrm{w}}_1\left(t\right)+{\int }_{R_{\gamma }}{g}_1\left(t,x_1\left(t\right),x_2\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ {\mathit{dx}}_2\left(t\right)=\left(-{\epsilon }_2{x}_2\left(t\right)+{\eta }_2{x}_1\left(t\right)x_2\left(t\right)+s_2\left(t,x_1\left(t\right),x_2\left(t\right)\right)\right)\mathit{dt}+\\ +b_2\left(t,x_1\left(t\right),x_2\left(t\right)\right)d{\mathrm{w}}_2\left(t\right)+{\int }_{R_{\gamma }}{g}_2\left(t,x_1\left(t\right),x_2\left(t\right),\gamma \right)\nu \left(\mathit{dt},\mathit{d\gamma }\right),\\ x_1\left(0\right)=N_1,x_2\left(0\right)=N_2,\end{array}\right.$$E40

with condition

$$u\left(t,\mathbf{x}\left(t\right)\right)=x_1^{-{\epsilon }_2}\left(t\right)e^{{\eta }_2{x}_1\left(t\right)}-{\mathit{Cx}}_2^{{\epsilon }_1}\left(t\right)e^{-{\eta }_1{x}_2\left(t\right)}.$$E41

Let us assume that ${\epsilon }_1=2$, ${\epsilon }_2=1$, ${\eta }_1={\eta }_2=1$, and $C=1$, and initial condition is $x\left(0\right)=y\left(0\right)=1$. The function $u\left(t,x,y\right)=x^{-1}{e}^x-y^2{e}^{-2y}$ is a first integral, $h\left(t,x,y\right)=y-x+e^{-t}$ and $q\left(t,x,y\right)=x$ are complementary functions.

We cannot find an analytical solution of the differential equations system

$$\left\{\begin{array}{l}\frac{\partial z_1\left(t,x,y,\gamma \right)}{\partial \gamma }=e^{-z_2\left(t,x,y,\gamma \right)}{z}_1\left(t,x,y,\gamma \right)z_2\left(t,x,y,\gamma \right)\left(z_2\left(t,x,y,\gamma \right)-2\right),\\ \frac{\partial z_2\left(t,x,y,\gamma \right)}{\partial \gamma }=-e^{-z_1\left(t,x,y,\gamma \right)}\left(1-z_1^{-1}\left(t,x,y,\gamma \right)\right).\end{array}\right.$$

Then, the constructed SDE system includes only Wiener perturbation:

$$\begin{array}{c}\left[\begin{array}{c}\mathit{dx}\left(t\right)\\ \mathit{dy}\left(t\right)\end{array}\right]=\left[\begin{array}{c}A\left(t,x\left(t\right),y\left(t\right)\right)+B\left(t,x\left(t\right),y\left(t\right)\right)+C\left(t,x\left(t\right),y\left(t\right)\right)\\ D\left(t,x\left(t\right),y\left(t\right)\right)+E\left(t,x\left(t\right),y\left(t\right)\right)\end{array}\right]\mathit{dt}+\\ +\left[\begin{array}{c}x\left(t\right)e^{-y\left(t\right)}\left(y^2\left(t\right)-2y\left(t\right)\right)\\ \frac{e^{x\left(t\right)}}{x\left(t\right)}-e^{x\left(t\right)}\end{array}\right]d\mathrm{w}\left(t\right),\end{array}$$E42

where

$$\begin{array}{l}A\left(t,x\left(t\right),y\left(t\right)\right)=0.5e^{-y\left(t\right)}x\left(t\right){\left(y^2\left(t\right)-2y\left(t\right)e^{-y\left(t\right)}\right)}^2,\\ B\left(t,x\left(t\right),y\left(t\right)\right)=0.5e^{-y\left(t\right)}\left(x\left(t\right)e^{x\left(t\right)}-e^{x\left(t\right)}{x}^{-1}\left(t\right)\right)\left(2-4y\left(t\right)+y^2\left(t\right)\right),\\ C\left(t,x\left(t\right),y\left(t\right)\right)=-\frac{e^{-t}{e}^{-y\left(t\right)}\left(y^2\left(t\right)-2y\left(t\right)\right)}{e^{x\left(t\right)}{x}^{-1}\left(t\right)-2y\left(t\right)e^{-y\left(t\right)}-e^{x\left(t\right)}{x}^{-2}\left(t\right)+y^2\left(t\right)e^{-y\left(t\right)}},\\ D\left(t,x\left(t\right),y\left(t\right)\right)=\frac{e^{-t}{e}^{x\left(t\right)}\left(x^{-1}\left(t\right)-x^{-2}\left(t\right)\right)}{e^{x\left(t\right)}{x}^{-1}\left(t\right)-2y\left(t\right)e^{-y\left(t\right)}-e^{x\left(t\right)}{x}^{-2}\left(t\right)+y^2\left(t\right)e^{-y\left(t\right)}},\\ E\left(t,x\left(t\right),y\left(t\right)\right)=-0.5e^{x\left(t\right)}{e}^{-y\left(t\right)}\left(y^2\left(t\right)-2y\left(t\right)\right)\left(1-x^{-1}\left(t\right)+x\left(t\right)-2+2x^{-1}\left(t\right)\right).\end{array}$$E43

Finally, we have the stochastic Lotka Volterra model associated to (38) ($N\left(t\right)=\left(N_1\left(t\right),N_2\left(t\right)\right)$):

$$\begin{array}{c}\left[\begin{array}{c}{\mathit{dN}}_1\left(t\right)\\ {\mathit{dN}}_2\left(t\right)\end{array}\right]=\left[\begin{array}{c}A\left(t,N\left(t\right)\right)+B\left(t,N\left(t\right)\right)+C\left(t,N\left(t\right)\right)\\ D\left(t,N\left(t\right)\right)+E\left(t,N\left(t\right)\right)\end{array}\right]\mathit{dt}+\\ +\left[\begin{array}{c}{N}_1\left(t\right)e^{-N_2\left(t\right)}\left(N_2^2\left(t\right)-2N_2\left(t\right)\right)\\ \frac{e^{N_1\left(t\right)}}{N_1\left(t\right)}-e^{N_1\left(t\right)}\end{array}\right]d\mathrm{w}\left(t\right),\end{array}$$E44

where $A\left(t,N\left(t\right)\right)$, $B\left(t,N\left(t\right)\right)$, $C\left(t,N\left(t\right)\right)$, $D\left(t,N\left(t\right)\right)$, $E\left(t,N\left(t\right)\right)$ are determined by Eq.(43).

Figures 3 and 4 show two realizations for numerical solution of Eq. (44).

Figure 3.

Numerical simulation 1 for solution of Eq.(44).

Figure 4.

Numerical simulation 2 for solution of Eq.(44).

Another examples of a differential equation system construction and models see in [25, 26, 27, 28, 29].

10. Conclusion

The invariant method widens horizons for constructing and researching into mathematical models of real systems with the invariants that hold out under any strong random disturbances.