Construction of analytical solutions to systems of two stochastic differential equations

Zenonas Navickas; Inga Telksniene; Tadas Telksnys; Romas Marcinkevicius; Minvydas Ragulskis

doi:10.1515/math-2023-0136

Article Open Access

Construction of analytical solutions to systems of two stochastic differential equations

, , , and

Published/Copyright: November 11, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

$Open Mathematics$

From the journal Open Mathematics Volume 21 Issue 1

Abstract

A scheme for the stochastization of systems of ordinary differential equations (ODEs) based on Itô calculus is presented in this article. Using the presented techniques, a system of stochastic differential equations (SDEs) can be constructed in such a way that eliminating the stochastic component yields the original system of ODEs. One of the main benefits of this scheme is the ability to construct analytical solutions to SDEs with the use of special vector-valued functions, which significantly differs from the randomization approach, which can only be applied via numerical integration. Moreover, using the presented techniques, a system of ODEs and SDEs can be constructed from a given diffusion function, which governs the uncertainty of a particular process.

Keywords: stochastic calculus; differential equation; Itô equation

MSC 2010: 60H10; 35R60

1 Introduction

Solutions to stochastic systems governed by multi-dimensional differential equations pose a clear interest due to important applications in science and engineering. In mathematics, the theory of stochastic processes remains to be an active topic of research for both theory and applications [1]. Although the theory of one-dimensional systems governed by stochastic differential equations (SDEs) is well developed, research into analytic solutions to nonlinear multi-dimensional SDEs continues. This is due not only to the problem’s mathematical relevance but also to the plethora of applications of SDEs in various fields of research. A short review of recent works is given below.

A scheme based on SDEs is used to optimize design strategies of chemical experiments in the study by Huang et al. [2]. SDE-based models are one of the main approaches in the modeling SARS-CoV-2 pandemic [3–5] and other infectious diseases [6,7] with various levels of detail. Due to the changing underwater landscape, SDEs are required to model underwater acoustic wave propagation [8]. Since perturbations in electrical power systems can lead to reduced quality, safety, and economy of electrical supply, SDEs are required to model such perturbations and predict their effects on the entire system [9].

Despite widespread applications, the construction of solutions to SDEs is not a trivial task. A well-known approach for approximating numerical solutions to stochastic engineering problems is the stochastic finite element method (SFEM). The straightforward approach to SFEM is based on the application of deterministic algorithms finite element method. For example, Monte-Carlo simulation methods were used in the study by Schuëller and Pradlwarter [10] to admit the stochastic characteristics in the form of a series of deterministic sample inputs (a clear advantage of such an approach is that no further solution methods have to be developed).

However, the Monte-Carlo simulation often requires a prohibitively large number of samples due to its slow convergence. The main attraction of the alternative approach to SFEM is the fact that it can solve the problem only once, providing the solution for the whole continuum of different realizations of the stochastic process [11].

Clearly, the SFEM has a few drawbacks too. It solves the problem in a high-dimensional physical-stochastic product space, whereby additional nodes and degrees of freedom must be created. Thus, the complexity of the problem grows exponentially as the number of random parameters increases [12,13].

Thus, the application of reduced-order techniques, allowing the reduction of the computational effort, is clearly the priority aim in stochastic computational methods. A good overview and comparison of different model reduction techniques is provided in the study by Besselink et al. [14].

The main objective of this article is to provide a technique for the stochastization of a reduced system of two coupled nonlinear differential equations and to construct analytical solutions to these stochastic equations. The apparent simplicity of the problem is deeply misleading. First, we are not interested in developing a numerical scheme for the solution of the stochastic system. Second, we are not interested in designing a randomization scheme that could yield estimates of martingales. On the contrary, the main objective of this study is to derive the necessary and sufficient conditions for the existence of analytic Ito integrals to the coupled system, and to construct the analytic solution itself.

Building on the research presented in the one-dimensional case [15], in this article, we consider the stochastization of the following systems of ordinary differential equations (ODEs):

(1) d y k ( t ) d t = P k ( y 1 ( t ) , y 2 ( t ) ) ,

where k = 1 , 2 . Let ω ( t ) denote a Wiener process [16]. As mentioned, the objective of this article is to construct a system of SDEs with respect to functions ξ k ( t ∣ α ) ( α is a scalar parameter) of the following form:

(2) d ξ k ( t ∣ α ) = a k ( ξ 1 , ξ 2 ∣ α ) d t + σ k ( ξ k ∣ α ) d ω ( t ) ,

where k = 1 , 2 ; t ≥ 0 ; 0 < α < 1 ; and the following pointwise limits hold true:

(3) lim α → 0 σ k ( ξ k ∣ α ) = 0 ;

(4) lim α → 0 a k ( ξ 1 , ξ 2 ∣ α ) = P k ( ξ 1 , ξ 2 ) ;

(5) lim α → 0 ξ k ( t ∣ α ) = y k ( t ) ,

i.e., as parameter α tends to 0, the system of SDEs equation (2) tends to equation (1) and conversely the solution to equation (2) tends to the solution of equation (1). The trajectories ξ k are dependent on the realization of the Wiener process ω . Note that in this article, we focus on the special case of the Wiener process ω ( t ) being the same for both values of k .

2 Preliminaries

2.1 Wiener and Itô processes

Consider a stochastic Wiener process ω ( t ) that satisfies the following well-known identities [16]:

(6) lim Δ t → 0 Δ ω ( t ) = 0 , lim Δ t → 0 ( Δ ω ( t ) ) k Δ t = δ 2 k ,

where δ j k is the Kronecker delta. Note that applied to differentials, the above equalities yield that [16]:

(7) ( d ω ( t ) ) 2 = d t .

An Itô process ξ ( t ) is a stochastic process that can be expressed via the following integral:

(8) ξ ( t ) = ξ ( 0 ) + ∫ 0 t a ( s , ξ ( s ) ) d s + ∫ 0 t σ ( s , ξ ( s ) ) d ω ( s ) ,

where ξ ( 0 ) is a scalar initial condition. Note that ω ( 0 ) = 0 , which means that the value ξ ( 0 ) is not a random variable, but a scalar. Functions a ( t , ξ ( t ) ) and σ ( t , ξ ( t ) ) represent drift and diffusion, respectively.

An equivalent definition of an Itô process is obtained by the differentiation of equation (8), which results in the following SDE:

(9) d ξ ( t ) = a ( t , ξ ) d t + σ ( t , ξ ) d ω ( t ) .

One of the most well-known results in stochastic calculus is Itô’s lemma, which states that a differential of f ( t , ξ ( t ) ) , where ξ ( t ) is defined via equation (9), is given as follows [16]:

(10) d f ( t , ξ ( t ) ) = ∂ f ( t , x ) ∂ t + a ( t , x ) ∂ f ( t , x ) ∂ x + σ 2 ( t , x ) 2 ∂ 2 f ( t , x ) ∂ x 2 x = ξ d t + σ ( t , ξ ) ∂ f ( t , x ) ∂ x x = ξ d ω ( t ) .

An approach to stochastization by applying equation (10) to a single ODE has already been discussed in the study by Navickas et al. [15]. In the remainder of this article, systems of stochastic equations will be considered. While the stochastization of a system of ODEs follows a similar idea to that of a single ODE, such generalization is far from trivial, as Itô’s lemma (10) is no longer true for a pair of stochastic processes.

2.2 Two-dimensional Itô’s lemma

Let k = 1 , 2 . Suppose that two Itô processes

(11) d ξ k ( t ) = a k ( ξ 1 , ξ 2 ) d t + σ k ( ξ k ) d ω ( t )

and functions F k : R × R → R are given. Note that the parameter α present in equation (2) is not yet introduced. Then, η k ( t ) = F k ( ξ 1 ( t ) , ξ 2 ( t ) ) are Itô processes given by [17]:

(12) d η k ( t ) = ∂ F k ( x , y ) ∂ x d x + ∂ F k ( x , y ) ∂ y d y + 1 2 ∂ 2 F k ( x , y ) ∂ x 2 ( d x ) 2 + 1 2 ∂ 2 F k ( x , y ) ∂ y 2 ( d y ) 2 + ∂ 2 F k ( x , y ) ∂ x ∂ y d x d y x = ξ 1 ( t ) y = ξ 2 ( t ) ,

where d ξ i ( t ) d ξ j ( t ) ( i , j = 1 , 2 ) is computed using the rules d t d t = d t d ω ( t ) = d ω ( t ) d t = 0 and d ω ( t ) d ω ( t ) = d t . Thus, inserting equation (11) into equation (12) yields

(13) d η k ( t ) = ∂ F k ( x , y ) ∂ x a 1 ( x , y ) + ∂ F k ( x , y ) ∂ y a 2 ( x , y ) + 1 2 ∂ 2 F k ( x , y ) ∂ x 2 ( σ 1 ( x ) ) 2 + 1 2 ∂ 2 F k ( x , y ) ∂ y 2 ( σ 2 ( y ) ) 2 + ∂ 2 F k ( x , y ) ∂ x ∂ y σ 1 ( x ) σ 2 ( y ) x = ξ 1 ( t ) y = ξ 2 ( t ) d t + ∂ F k ( x , y ) ∂ x σ 1 ( x ) + ∂ F k ( x , y ) ∂ y σ 2 ( y ) x = ξ 1 ( t ) y = ξ 2 ( t ) d ω ( t ) .

Since F k are arbitrary functions, the initial conditions for η k are equal to: η k ( 0 ) = η k 0 = F k ( ξ 1 ( 0 ) , ξ 2 ( 0 ) ) .

2.3 Solutions to the systems of SDEs of type (11)

Let us first consider the following system of SDEs where drift and diffusion for both Itô processes are constant:

(14) d η k ( t ) = a ^ k d t + σ ^ k d ω ( t ) , a ^ k , σ ^ k ∈ R ; k = 1 , 2 .

In this case, solution η k ( t ) to equation (14) can be obtained directly by applying the Itô integral (8):

(15) η k ( t ) = η k 0 + a ^ k t + σ ^ k ω ( t ) , η k 0 ∈ R .

Now let us consider a generalized version of equation (14) that is defined in equation (11). Suppose that it is possible to determine a vector function F ( x , y ) = ( F 1 ( x , y ) , F 2 ( x , y ) ) and its inverse F − 1 ( x , y ) = G ( x , y ) = ( G 1 ( x , y ) , G 2 ( x , y ) ) .

The functions F and G satisfy the canonical relations for inverse mappings:

(16) F ( G ( x , y ) ) = ( x , y ) , G ( F ( x , y ) ) = ( x , y ) , ∀ x , y ∈ R .

If the vector functions F and G are considered coordinate-wise, then it can be denoted that:

(17) F ( ξ 1 ( t ) , ξ 2 ( t ) ) = ( η 1 ( t ) , η 2 ( t ) ) , G ( η 1 ( t ) , η 2 ( t ) ) = ( ξ 1 ( t ) , ξ 2 ( t ) ) .

Then, combining equations (15) and (17) yields the solution to equation (11):

(18) ξ k ( t ) = G k ( η 10 + a ^ 1 t + σ ^ 1 ω ( t ) , η 20 + a ^ 2 t + σ ^ 2 ω ( t ) ) .

2.3.1 Derivation of the functions F k ( x , y ) , k = 1 , 2

Note that the first relation in equation (17) implies equation (13), which together with equation (14) yields

(19) a ^ k = ∂ F k ( x , y ) ∂ x a 1 ( x , y ) + ∂ F k ( x , y ) ∂ y a 2 ( x , y ) + 1 2 ∂ 2 F k ( x , y ) ∂ x 2 ( σ 1 ( x ) ) 2 + 1 2 ∂ 2 F k ( x , y ) ∂ y 2 ( σ 2 ( y ) ) 2 + ∂ 2 F k ( x , y ) ∂ x ∂ y σ 1 ( x ) σ 2 ( y ) x = ξ 1 ( t ) y = ξ 2 ( t ) ;

(20) σ ^ k = ∂ F k ( x , y ) ∂ x σ 1 ( x ) + ∂ F k ( x , y ) ∂ y σ 2 ( y ) x = ξ 1 ( t ) y = ξ 2 ( t ) ,

where ξ k = ξ k ( t ) , a ^ k , σ ^ k ∈ R ; k = 1 , 2 .

2.3.1.1. Solution of equation (20)

Let us denote both terms in equation (20) by γ 1 k ( ξ 1 , ξ 2 ) and γ 2 k ( ξ 1 , ξ 2 ) , respectively, i.e.,

(21) σ ^ k = γ 1 k ( ξ 1 , ξ 2 ) + γ 2 k ( ξ 1 , ξ 2 ) ,

which yields

(22) ∂ F k ( x , y ) ∂ x x = ξ 1 ( t ) y = ξ 2 ( t ) = γ 1 k ( ξ 1 , ξ 2 ) σ 1 ( ξ 1 ) ; ∂ F k ( x , y ) ∂ y x = ξ 1 ( t ) y = ξ 2 ( t ) = γ 2 k ( ξ 1 , ξ 2 ) σ 2 ( ξ 2 ) .

Note that mixed partial derivatives of equation (22) must satisfy the relation:

(23) ∂ 2 F k ( x , y ) ∂ x ∂ y x = ξ 1 ( t ) y = ξ 2 ( t ) = ∂ 2 F k ( x , y ) ∂ y ∂ x x = ξ 1 ( t ) y = ξ 2 ( t ) ,

which yields the following PDE:

(24) ∂ γ 1 k ( x , y ) ∂ y σ 2 ( y ) x = ξ 1 ( t ) y = ξ 2 ( t ) = ∂ γ 2 k ( x , y ) ∂ x σ 1 ( x ) x = ξ 1 ( t ) y = ξ 2 ( t ) .

From equation (21), it is obtained that γ 2 k = σ ˜ k − γ 1 k and likewise γ 1 k = σ ˜ k − γ 2 k . Inserting these relations into equation (24) yields two PDEs:

(25) ∂ γ 1 k ∂ x σ 1 ( x ) + ∂ γ 1 k ∂ y σ 2 ( y ) x = ξ 1 ( t ) y = ξ 2 ( t ) = 0 ;

(26) ∂ γ 2 k ∂ x σ 1 ( x ) + ∂ γ 2 k ∂ y σ 2 ( y ) x = ξ 1 ( t ) y = ξ 2 ( t ) = 0 .

Solving the above equations via the method of characteristics yields the following result [18]:

(27) γ 1 k ( ξ 1 , ξ 2 ) = γ ^ 1 k + U k ( ξ 1 , ξ 2 ) ; γ 2 k ( ξ 1 , ξ 2 ) = γ ^ 2 k − U k ( ξ 1 , ξ 2 ) ,

where

(28) U k ( ξ 1 , ξ 2 ) = H k ∫ ξ 10 ξ 1 1 σ 1 ( u ) d u − ∫ ξ 20 ξ 2 1 σ 2 ( v ) d v ,

γ ^ 1 k , γ ^ 2 k ∈ R , γ ^ 1 k + γ ^ 2 k = σ ^ k , and H k ( x ) , k = 1 , 2 is an arbitrary differentiable function.

2.3.1.2. Solution of equation (19)

Computing second-order derivatives of equation (22) and inserting the obtained expressions into equation (19) yields:

(29) a ^ k = γ 1 k ( x , y ) σ 1 ( x ) a 1 ( x , y ) + γ 2 k ( x , y ) σ 2 ( y ) a 2 ( x , y ) + 1 2 ∂ γ 1 k ( x , y ) ∂ x σ 1 ( x ) − γ 1 k ( x , y ) d σ 1 ( x ) d x + 1 2 ∂ γ 2 k ( x , y ) ∂ y σ 2 ( y ) − γ 2 k ( x , y ) d σ 2 ( y ) d y + 1 2 ∂ γ 1 k ( x , y ) ∂ y σ 2 ( y ) + ∂ γ 2 k ( x , y ) ∂ x σ 1 ( x ) x = ξ 1 ( t ) y = ξ 2 ( t ) ,

where the term ∂ 2 F k ( x , y ) ∂ x ∂ y x = ξ 1 ( t ) y = ξ 2 ( t ) has been split into 1 2 ∂ 2 F k ( x , y ) ∂ x ∂ y + ∂ 2 F k ( x , y ) ∂ y ∂ x x = ξ 1 ( t ) y = ξ 2 ( t ) . Applying equation (27) transforms equation (29) into:

(30) a ^ k = γ ^ 1 k + U k ( x , y ) σ 1 ( x ) a 1 ( x , y ) + γ ^ 2 k − U k ( x , y ) σ 2 ( y ) a 2 ( x , y ) + 1 2 ∂ U k ( x , y ) ∂ x σ 1 ( x ) − ( γ ^ 1 k + U k ( x , y ) ) d σ 1 ( x ) d x + 1 2 − ∂ U k ( x , y ) ∂ y σ 2 ( y ) − ( γ ^ 2 k − U k ( x , y ) ) d σ 2 ( y ) d y + 1 2 ∂ U k ( x , y ) ∂ y σ 2 ( y ) − ∂ U k ( x , y ) ∂ x σ 1 ( x ) x = ξ 1 ( t ) y = ξ 2 ( t ) ,

which can be simplified as follows:

(31) a ^ k = γ ^ 1 k + U k ( x , y ) σ 1 ( x ) a 1 ( x , y ) + γ ^ 2 k − U k ( x , y ) σ 2 ( y ) a 2 ( x , y ) − 1 2 γ ^ 1 k d σ 1 ( x ) d x + γ ^ 2 k d σ 2 ( y ) d y − 1 2 U k ( x , y ) d σ 1 ( x ) d x − d σ 2 ( y ) d y x = ξ 1 ( t ) y = ξ 2 ( t ) .

Derivations presented above could be summarized as the following Lemma.

Lemma 2.1

Consider two systems of SDEs:

(32) d ξ k ( t ) = a k ( ξ 1 , ξ 2 ) d t + σ k ( ξ k ) d ω ( t ) , k = 1 , 2 ,

and

(33) d η k ( t ) = a ^ k d t + σ ^ k d ω ( t ) , a ^ k , σ ^ k ∈ R ; k = 1 , 2 .

Suppose that there exist values such as γ ^ 1 k and γ ^ 2 k and function U k ( x , y ) , k = 1 , 2 , that satisfy equation (31). Then, function F k ( x , y ) , satisfying the relation

(34) F k ( ξ 1 , ξ 2 ) = η k , k = 1 , 2 ,

reads as follows:

(35) F k ( ξ 1 , ξ 2 ) = γ ^ 1 k ∫ ξ ^ 10 ξ 1 1 σ 1 ( s ) d s + γ ^ 2 k ∫ ξ ^ 20 ξ 2 1 σ 2 ( t ) d t + Ψ k ( ξ 1 , ξ 2 ) ,

where γ ^ 1 k + γ ^ 2 k = σ ^ k , ξ ^ 10 , ξ ^ 20 , ξ 10 , ξ 20 ∈ R , and the function Ψ k ( ξ 1 , ξ 2 ) satisfies:

(36) Ψ k ( ξ 1 , ξ 2 ) = ∫ ξ 10 ξ 1 U k ( s , ξ 2 ) σ 1 ( s ) d s = − ∫ ξ 20 ξ 2 U k ( ξ 1 , t ) σ 2 ( t ) d t .

Proof

Integrating equation (22) with respect to ξ 1 and ξ 2 yields

(37) F k ( ξ 1 , ξ 2 ) = ∫ ξ 10 ξ 1 γ 1 k ( s , ξ 2 ) σ 1 ( s ) d s + C 1 ( ξ 2 ) ; F k ( ξ 1 , ξ 2 ) = ∫ ξ 20 ξ 2 γ 2 k ( ξ 1 , t ) σ 2 ( t ) d t + C 2 ( ξ 1 ) ,

where C 1 and C 2 are arbitrary univariate functions. Noting that equation (27) holds true, the above equations become:

(38) F k ( ξ 1 , ξ 2 ) = ∫ ξ 10 ξ 1 γ ^ 1 k σ 1 ( s ) d s + ∫ ξ 10 ξ 1 U k ( s , ξ 2 ) σ 1 ( s ) d s + C 1 ( ξ 2 ) ; F k ( ξ 1 , ξ 2 ) = ∫ ξ 20 ξ 2 γ ^ 2 k σ 2 ( t ) d t − ∫ ξ 20 ξ 2 U k ( ξ 1 , t ) σ 2 ( t ) d t + C 2 ( ξ 1 ) .

This yields that for arbitrary constants C ^ 1 and C ^ 2 , the following relations hold true:

(39) C 1 ( ξ 2 ) = ∫ ξ 20 ξ 2 γ ^ 2 k σ 2 ( t ) d t + C ^ 1 ;

(40) C 2 ( ξ 1 ) = ∫ ξ 10 ξ 1 γ ^ 1 k σ 1 ( s ) d s + C ^ 2 .

Inserting equations (39) and (40) into equation (38), we obtain

(41) F k ( ξ 1 , ξ 2 ) = ∫ ξ 10 ξ 1 γ ^ 1 k σ 1 ( s ) d s + ∫ ξ 20 ξ 2 γ ^ 2 k σ 2 ( t ) d t + C ^ 1 + ∫ ξ 10 ξ 1 U k ( s , ξ 2 ) σ 1 ( s ) d s ; F k ( ξ 1 , ξ 2 ) = ∫ ξ 20 ξ 2 γ ^ 2 k σ 2 ( t ) d t + ∫ ξ 10 ξ 1 γ ^ 1 k σ 1 ( s ) d s + C ^ 2 − ∫ ξ 20 ξ 2 U k ( ξ 1 , t ) σ 2 ( t ) d t .

Since constants C ^ 1 and C ^ 2 are arbitrary, they can be chosen in such a way to rewrite the first two integrals on both equations of equation (41) to have different lower bounds equal to ξ ^ 10 and ξ ^ 20 :

(42) F k ( ξ 1 , ξ 2 ) = ∫ ξ ^ 10 ξ 1 γ ^ 1 k σ 1 ( s ) d s + ∫ ξ ^ 20 ξ 2 γ ^ 2 k σ 2 ( t ) d t + ∫ ξ 10 ξ 1 U k ( s , ξ 2 ) σ 1 ( s ) d s ; F k ( ξ 1 , ξ 2 ) = ∫ ξ ^ 20 ξ 2 γ ^ 2 k σ 2 ( t ) d t + ∫ ξ ^ 10 ξ 1 γ ^ 1 k σ 1 ( s ) d s − ∫ ξ 20 ξ 2 U k ( ξ 1 , t ) σ 2 ( t ) d t .

The above equation yields that

(43) ∫ ξ 10 ξ 1 U k ( s , ξ 2 ) σ 1 ( s ) d s = − ∫ ξ 20 ξ 2 U k ( ξ 1 , t ) σ 2 ( t ) d t = Ψ k ( ξ 1 , ξ 2 ) ,

which finalizes the proof.□

3 Main results

The goal of this section is to prove the following theorem:

Theorem 3.1

Let the following system of ODEs be given as follows:

(44) d y k ( t ) d t = P k ( y 1 ( t ) , y 2 ( t ) ) , k = 1 , 2 .

Then, the system of SDEs

(45) d ξ k ( t ∣ α ) = a k ( ξ 1 , ξ 2 ∣ α ) d t + σ k ( ξ k ∣ α ) d ω ( t ) , k = 1 , 2 .

with t ≥ 0 , 0 < α < 1 , is a stochastization of equation (44) (i.e., it satisfies equations (3)–(5)) if relations (65), (67), and (69) hold true.

Moreover, the solution to the system of SDEs (45) reads as follows:

(46) ξ k ( t ∣ α ) = G k η 10 + S 1 α + S 1 ( + ) t + ( γ ^ 11 + γ ^ 21 ) ω ( t ) , η 20 + S 2 α + S 2 ( + ) t + ( γ ^ 12 + γ ^ 22 ) ω ( t ) ,

where η 10 , η 20 ∈ R , and k = 1 , 2 .

The notations required for Theorem 3.1 are introduced in the remaining parts of this section, which are dedicated entirely to the theorem’s proof.

3.1 Derivation of the functions G k ( x , y ) , k = 1 , 2

Suppose that the conditions of Lemma 2.1 hold and

(47) γ ^ 11 γ ^ 22 − γ ^ 12 γ ^ 21 ≠ 0 .

In addition, let functions g 1 ( ξ 1 ) and g 2 ( ξ 2 ) satisfy the following equations:

(48) ∫ ξ ^ 10 g 1 ( ξ 1 ) 1 σ 1 ( s ) d s = ξ 1 ; ∫ ξ ^ 20 g 2 ( ξ 2 ) 1 σ 2 ( t ) d t = ξ 2 .

Denote

(49) G k ( ξ 1 , ξ 2 ) = g k ( λ k ( ξ 1 , ξ 2 ) ) = g k ( γ ˜ 1 k ξ 1 + γ ˜ 2 k ξ 2 + Ψ ˜ k ( ξ 1 , ξ 2 ) ) ,

where Ψ ˜ k = Ψ ˜ k ( ξ 1 , ξ 2 ) and coefficients γ ˜ 1 k , γ ˜ 2 k ∈ R satisfy

(50) γ ˜ 11 γ ˜ 12 γ ˜ 21 γ ˜ 22 ⋅ γ ^ 11 γ ^ 12 γ ^ 21 γ ^ 22 = I ,

where I is an identity matrix.

Lemma 3.1

Vector function G ( ξ 1 , ξ 2 ) = ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) is an inverse of F ( ξ 1 , ξ 2 ) = ( F 1 ( ξ 1 , ξ 2 ) , F 2 ( ξ 1 , ξ 2 ) ) if the following system of equations holds true:

(51) γ ^ 1 k Ψ ˜ 1 ( ξ 1 , ξ 2 ) + γ ^ 2 k Ψ ˜ 2 ( ξ 1 , ξ 2 ) + Ψ k ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) = 0 ,

where k = 1 , 2 .

Proof

In order for the vector function G ( ξ 1 , ξ 2 ) = ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) to be an inverse of F ( ξ 1 , ξ 2 ) = ( F 1 ( ξ 1 , ξ 2 ) , F 2 ( ξ 1 , ξ 2 ) ) , the following relations must be satisfied:

(52) F k ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) = ξ k , k = 1 , 2 .

Inserting equations (35) and (49) into equation (52), we obtain

(53) γ ^ 1 k ∫ ξ ^ 10 g 1 ( λ 1 ( ξ 1 , ξ 2 ) ) 1 σ 1 ( s ) d s + γ ^ 2 k ∫ ξ ^ 20 g 2 ( λ 2 ( ξ 1 , ξ 2 ) ) 1 σ 2 ( t ) d t + Ψ k ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) = ξ k .

Using equation (48) converts equation (53) into

(54) γ ^ 1 k λ 1 ( ξ 1 , ξ 2 ) + γ ^ 2 k λ 2 ( ξ 1 , ξ 2 ) + Ψ k ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) = ξ k .

Inserting λ k ( ξ 1 , ξ 2 ) = γ ˜ 1 k ξ 1 + γ ˜ 2 k ξ 2 + Ψ ˜ k ( ξ 1 , ξ 2 ) into equation (54) yields:

(55) γ ^ 1 k ( γ ˜ 11 ξ 1 + γ ˜ 21 ξ 2 + Ψ ˜ 1 ( ξ 1 , ξ 2 ) ) + γ ^ 2 k ( γ ˜ 12 ξ 1 + γ ˜ 22 ξ 2 + Ψ ˜ 2 ( ξ 1 , ξ 2 ) ) + Ψ k ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) = ξ k .

Rearranged equation (55) reads as follows::

(56) ξ 1 ( γ ^ 1 k γ ˜ 11 + γ ^ 2 k γ ˜ 12 ) + ξ 2 ( γ ^ 1 k γ ˜ 21 + γ ^ 2 k γ ˜ 22 ) + γ ^ 1 k Ψ ˜ 1 ( ξ 1 , ξ 2 ) + γ ^ 2 k Ψ ˜ 2 ( ξ 1 , ξ 2 ) + Ψ k ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) = ξ k .

Applying equation (50) results in

(57) ξ k + γ ^ 1 k Ψ ˜ 1 ( ξ 1 , ξ 2 ) + γ ^ 2 k Ψ ˜ 2 ( ξ 1 , ξ 2 ) + Ψ k ( G 1 ( ξ 1 , ξ 2 ) , G 2 ( ξ 1 , ξ 2 ) ) = ξ k ,

which is equivalent to equation (51).□

3.2 Relation between systems of SDEs with variable and constant coefficients

In this section, a technique for the construction of two systems of SDEs:

(58) d ξ k ( t ) = a k ( ξ 1 , ξ 2 ) d t + σ k ( ξ k ) d ω ( t ) , k = 1 , 2

and

(59) d η k ( t ) = a ^ k d t + σ ^ k d ω ( t ) , a ^ k , σ ^ k ∈ R ; k = 1 , 2 ,

as well as functions F k ( x , y ) and G k ( x , y ) , satisfying

(60) F k ( ξ 1 ( t ) , ξ 2 ( t ) ) = η k ( t ) , G k ( η 1 ( t ) , η 2 ( t ) ) = ξ k ( t ) ,

is presented. Schematic diagrams of this technique are depicted in Figures 1 and 2.

$Figure 1 A schematic diagram, illustrating functions and parameters considered in the technique presented in Section 3.2. Given functions and parameters are shown in green, whereas the derived ones are depicted in red.$

Figure 1

A schematic diagram, illustrating functions and parameters considered in the technique presented in Section 3.2. Given functions and parameters are shown in green, whereas the derived ones are depicted in red.

$Figure 2 A schematic diagram of the technique presented in Section 3.2. Given functions and parameters are shown in green, whereas the derived ones are depicted in red.$

Figure 2

A schematic diagram of the technique presented in Section 3.2. Given functions and parameters are shown in green, whereas the derived ones are depicted in red.

Consider the case where the following functions and parameters are given:

Functions σ k ( ξ k ) ;
Functions U k ( ξ 1 , ξ 2 ) that satisfy the relation (36);
Parameters γ ^ 1 k , γ ^ 2 k ∈ R that satisfy the relation (47);
Parameters a ^ k ∈ R ,

where k = 1 , 2 .

Then, the results of the previous sections yield functions and parameters a k ( ξ 1 , ξ 2 ) , σ ^ k , F k ( x , y ) , and G k ( x , y ) as follows:

Functions a k ( ξ 1 , ξ 2 ) are computed by solving the system of linear equations (31), yielding:
(61) a k ( ξ 1 , ξ 2 ) = Δ k Δ , k = 1 , 2 ,
where
(62) Δ = A 1 A 2 A 3 A 4 ; Δ 1 = A 5 A 2 A 6 A 4 ; Δ 2 = A 1 A 5 A 3 A 6 ; A 1 = γ ^ 11 + U 1 ( ξ 1 , ξ 2 ) σ 1 ( ξ 1 ) ; A 2 = γ ^ 21 − U 1 ( ξ 1 , ξ 2 ) σ 2 ( ξ 2 ) ; A 3 = γ ^ 12 + U 2 ( ξ 1 , ξ 2 ) σ 1 ( ξ 1 ) ; A 4 = γ ^ 22 − U 2 ( ξ 1 , ξ 2 ) σ 2 ( ξ 2 ) ; A 5 = a ^ 1 + 1 2 γ ^ 11 d σ 1 ( ξ 1 ) d ξ 1 + γ ^ 21 d σ 2 ( ξ 2 ) d ξ 2 + 1 2 U 1 ( ξ 1 , ξ 2 ) d σ 1 ( ξ 1 ) d ξ 1 − d σ 2 ( ξ 2 ) d ξ 2 ; A 6 = a ^ 2 + 1 2 γ ^ 12 d σ 1 ( ξ 1 ) d ξ 1 + γ ^ 22 d σ 2 ( ξ 2 ) d ξ 2 + 1 2 U 2 ( ξ 1 , ξ 2 ) d σ 1 ( ξ 1 ) d ξ 1 − d σ 2 ( ξ 2 ) d ξ 2 .
Parameters σ ^ k = γ ^ 1 k + γ ^ 2 k ;
Functions F k ( x , y ) are obtained by directly applying equation (35);
Functions G k ( x , y ) are obtained by applying equation (49), where:
1. Functions g k ( ξ k ) are computed from equation (48);
2. Parameters γ ˜ 1 k , γ ˜ 2 k ∈ R are derived from equation (50);
3. Functions Ψ ˜ k ( ξ 1 , ξ 2 ) are obtained by solving the system of algebraic equations (51).

3.3 Stochastization of the systems of ODEs

Consider the following system of ODEs:

(63) d y k ( t ) d t = P k ( y 1 ( t ) , y 2 ( t ) ) ,

where k = 1 , 2 . The objective of this section is to construct a system of SDEs of the form:

(64) d ξ k ( t ∣ α ) = a k ( ξ 1 , ξ 2 ∣ α ) d t + σ k ( ξ k ∣ α ) d ω ( t ) ,

where k = 1 , 2 ; t ≥ 0 ; 0 < α < 1 , and the relations (3)–(5) hold true.

Naturally, more than one set of a k ( ξ 1 , ξ 2 ∣ α ) and σ k ( ξ k ∣ α ) exist that satisfy equations (3)–(5). In this article, the most straightforward assumption is considered: the parameter α linearly increases and decreases the stochastization intensity, resulting in an increase of both additive and multiplicative noise. Let us denote:

(65) σ k ( ξ k ∣ α ) = α h k ( ξ k ) ,

where k = 1 , 2 ; 0 < α < 1 ; and h k ( x ) is a univariate function. Higher values of α lead to a higher level of stochastization compared to lower values of the same parameter.

Note that equation (65) ensures that equation (3) is satisfied. Next, let:

(66) a ^ k = S k α + S k ( + ) ,

where S k , S k ( + ) ∈ R . Note that the singularity α = 0 need not be considered, since the SDE for ξ k ( t ) becomes ODEs (due to equations (3)–(5)) and η k ( t ) is not required.

Techniques presented in the previous section yield that

(67) a k ( ξ 1 , ξ 2 ∣ α ) = Δ k Δ = Δ ˜ k + α 2 Δ ˜ k ( + ) Δ ˜ , k = 1 , 2 ,

where

Δ ˜ = A ˜ 1 A ˜ 2 A ˜ 3 A ˜ 4 ; Δ ˜ 1 = A ˜ 5 A ˜ 2 A ˜ 6 A ˜ 4 ; Δ ˜ 2 = A ˜ 1 A ˜ 5 A ˜ 3 A ˜ 6 ; Δ ˜ 1 ( + ) = A ˜ 7 A ˜ 2 A ˜ 8 A ˜ 4 ; Δ ˜ 2 ( + ) = A ˜ 1 A ˜ 7 A ˜ 3 A ˜ 8 ;

(68) A ˜ 1 = γ ^ 11 + U 1 ( ξ 1 , ξ 2 ) h 1 ( ξ 1 ) ; A ˜ 2 = γ ^ 21 − U 1 ( ξ 1 , ξ 2 ) h 2 ( ξ 2 ) ; A ˜ 3 = γ ^ 12 + U 2 ( ξ 1 , ξ 2 ) h 1 ( ξ 1 ) ; A ˜ 4 = γ ^ 22 − U 2 ( ξ 1 , ξ 2 ) h 2 ( ξ 2 ) ; A ˜ 5 = S 1 ; A ˜ 6 = S 2 ; A ˜ 7 = S 1 ( + ) + α 2 γ ^ 11 d h 1 ( ξ 1 ) d ξ 1 + γ ^ 21 d h 2 ( ξ 2 ) d ξ 2 + α 2 U 1 ( ξ 1 , ξ 2 ) d h 1 ( ξ 1 ) d ξ 1 − d h 2 ( ξ 2 ) d ξ 2 ; A ˜ 8 = S 2 ( + ) + α 2 γ ^ 12 d h 1 ( ξ 1 ) d ξ 1 + γ ^ 22 d h 2 ( ξ 2 ) d ξ 2 + α 2 U 2 ( ξ 1 , ξ 2 ) d h 1 ( ξ 1 ) d ξ 1 − d h 2 ( ξ 2 ) d ξ 2

and U k ( ξ 1 , ξ 2 ) , γ ^ 1 k , γ ^ 2 k ∈ R satisfy the relations (31), (36), and (47).

Then, equation (4) holds true if the following condition is satisfied:

(69) Δ ˜ k Δ ˜ = P k ( ξ 1 , ξ 2 ) .

Theorem 3.1

Let the following system of ODEs be given:

(70) d y k ( t ) d t = P k ( y 1 ( t ) , y 2 ( t ) ) , k = 1 , 2 .

Then, the system of SDEs

(71) d ξ k ( t ∣ α ) = a k ( ξ 1 , ξ 2 ∣ α ) d t + σ k ( ξ k ∣ α ) d ω ( t ) , k = 1 , 2 .

with t ≥ 0 and 0 < α < 1 is a stochastization of equation (70) (i.e., it satisfies equations (3)–(5)) if relations (65), (67), and (69) hold true.

Moreover, the solution to the system of SDEs (71) reads as follows:

(72) ξ k ( t ∣ α ) = G k η 10 + S 1 α + S 1 ( + ) t + ( γ ^ 11 + γ ^ 21 ) ω ( t ) , η 20 + S 2 α + S 2 ( + ) t + ( γ ^ 12 + γ ^ 22 ) ω ( t ) ,

where η 10 , η 20 ∈ R , k = 1 , 2 .

Proof

The first result of the theorem follows directly from the derivations presented above. The solution equation (72) is obtained by inserting equation (66) as well as σ ^ k = γ ^ 1 k + γ ^ 2 k into equation (18).□

A schematic diagram of Theorem 3.1 is depicted in Figure 3.

$Figure 3 A schematic diagram of Theorem 3.1. Note, that k = 1 , 2 k=1,2 .$

Figure 3

A schematic diagram of Theorem 3.1. Note, that k = 1 , 2 .

4 Computational experiments

In this section, the stochastization of a system of ODEs is derived using the techniques presented in previous sections and compared to another approach of introducing randomness into ODE systems.

4.1 Stochastization of an ODE system

First, the following functions and parameters are selected:

The simplest nonlinear case of diffusion functions (65) – quadratic polynomials – is considered
(73) σ 1 ( ξ 1 ∣ α ) = α ( 3 − ξ 1 ) ( ξ 1 − 2 ) ; σ 2 ( ξ 2 ∣ α ) = α ( 2 − ξ 2 ) ( ξ 2 − 1 ) ,
where 0 < α < 1 . Note that equation (3) is satisfied.
Then, parameters S k = 1 and S k ( + ) = 0 , k = 1 , 2 , are selected, which when inserted into equation (66) yields
(74) a ^ k = 1 α , k = 1 , 2 .
The next step is selecting functions U 1 ( ξ 1 , ξ 2 ) and U 2 ( ξ 1 , ξ 2 ) do satisfy the relation (36) for some constants ξ 10 and ξ 20 . Note that this selection is not unique, and the following is an example:
(75) U 1 ( ξ 1 , ξ 2 ) = ( ξ 1 − 2 ) ( 2 − ξ 2 ) ( 3 − ξ 1 ) ( ξ 2 − 1 ) ; U 2 ( ξ 1 , ξ 2 ) = 0 .
The above functions satisfy the relation (36) when ξ 1 0 = ξ 2 0 = 2 .
As a final step, four constants γ ^ j k must be selected. As in the previous case, this selection is not unique:
(76) γ ^ 11 = 13 ; γ ^ 12 = 4 ; γ ^ 21 = 2 ; γ ^ 22 = 1 .

Next, functions and parameters shown in red in Figure 1 are derived as described in Sections 3.2 and 3.3:

Drift functions (67) read as follows:
(77) a 1 ( ξ 1 , ξ 2 ∣ α ) = ( ξ 1 − 2 ) ( − 3 + ξ 1 ) ( 20 ξ 2 − 30 ) ξ 1 − 50 ξ 2 + 70 ( 20 α 2 ξ 1 2 ξ 2 − 30 α 2 ξ 1 2 − 100 α 2 ξ 1 ξ 2 + 145 α 2 ξ 1 + 125 α 2 ξ 2 − 175 α 2 + 2 ξ 1 − 2 ξ 2 − 2 ) ; a 2 ( ξ 1 , ξ 2 ∣ α ) = 20 ( − 2 + ξ 2 ) ( ξ 2 − 1 ) ( 20 ξ 1 − 50 ) ξ 2 − 30 ξ 1 + 70 × α 2 ( ξ 1 − 5 ⁄ 2 ) ξ 2 2 + − 3 α 2 ξ 1 + 29 α 2 4 − ξ 1 + 29 10 ξ 2 + 9 ⁄ 4 α 2 ξ 1 − 21 α 2 4 + 11 ξ 1 10 − 31 10 .
Parameters σ ^ k , k = 1 , 2 are computed as σ ^ 1 = γ ^ 11 + γ ^ 21 = 15 and σ ^ 2 = γ ^ 12 + γ ^ 22 = 5 .
Selecting ξ ^ 10 = ξ ^ 20 = − ∞ (note that these are distinct from ξ 10 and ξ 20 ) and applying equation (35) yields functions F k ( x , y ) :
(78) F 1 ( x , y ) = 1 α ln ( x − 2 ) 13 ( y − 1 ) 2 ( 3 − x ) 13 ( 2 − y ) 2 + ( − 2 + y ) ( x − 2 ) ( y − 1 ) α ( − 3 + x ) ; F 2 ( x , y ) = 1 α ln ( x − 2 ) 4 ( y − 1 ) ( 3 − x ) 4 ( 2 − y ) .
Applying equation (49) yields functions G k ( x , y ) :
(79) G 1 ( x , y ) = 3 e 1 5 x α e − 2 5 y α e − 1 5 W ( e α ( x − 3 y ) ) + 2 1 + e 1 5 x α e − 2 5 y α e − 1 5 W ( e α ( x − 3 y ) ) G 2 ( x , y ) = 2 e − 4 5 x α e 13 y α 5 e 4 5 W ( e α ( x − 3 y ) ) + 1 1 + e − 4 5 x α e 13 y α 5 e 4 5 W ( e α ( x − 3 y ) ) ,
where W ( ⋅ ) is Lambert W function (the solution of equation y exp y = x with respect to y for a given x ) [19].

The functions P 1 and P 2 that define the ODE system can be obtained by using condition (4):

(80) P 1 ( ξ 1 , ξ 2 ) = lim α → 0 a 1 ( ξ 1 , ξ 2 ∣ α ) = ( ξ 1 − 2 ) ( − 3 + ξ 1 ) ( ξ 1 − ξ 2 − 1 ) ( 10 ξ 2 − 15 ) ξ 1 − 25 ξ 2 + 35 ; P 2 ( ξ 1 , ξ 2 ) = lim α → 0 a 2 ( ξ 1 , ξ 2 ∣ α ) = − ( 10 ξ 1 ξ 2 − 11 ξ 1 − 29 ξ 2 + 31 ) ( ξ 2 − 1 ) ( − 2 + ξ 2 ) 10 ξ 1 ξ 2 − 15 ξ 1 − 25 ξ 2 + 35 .

Then, according to Theorem 3.1, the system of SDEs

(81) d ξ k ( t ∣ α ) = a k ( ξ 1 , ξ 2 ∣ α ) d t + σ k ( ξ k ∣ α ) d ω ( t ) , k = 1 , 2 ,

is a stochastization of the system of ODEs:

(82) d y k ( t ) d t = P k ( y 1 ( t ) , y 2 ( t ) ) , k = 1 , 2 ,

where a k ( ξ 1 , ξ 2 ∣ α ) , σ k ( ξ k ∣ α ) , and P k ( y 1 ( t ) , y 2 ( t ) ) , k = 1 , 2 , are defined by equations (77), (73), and (80), respectively, and t ≥ 0 and 0 < α < 1 .

Moreover, the analytical solution to the system of SDEs (81) can be obtained by applying (72):

(83) ξ 1 ( t ∣ α ) = 3 e η 10 α 5 e − t 5 e ω ( t ) α e − 2 η 20 α 5 e − W ( e ( η 10 − 3 η 20 ) α − 2 t ) 5 + 2 1 + e η 10 α 5 e − t 5 e ω ( t ) α e − 2 η 20 α 5 e − W ( e ( η 10 − 3 η 20 ) α − 2 t ) 5 ; ξ 2 ( t ∣ α ) = 2 e − 4 η 10 α 5 e 9 t 5 e ω ( t ) α e 13 η 20 α 5 e 4 W ( e ( η 10 − 3 η 20 ) α − 2 t ) 5 + 1 1 + e − 4 η 10 α 5 e 9 t 5 e ω ( t ) α e 13 η 20 α 5 e 4 W ( e ( η 10 − 3 η 20 ) α − 2 t ) 5 ,

where η 10 , η 20 ∈ R are values satisfying the following relation:

(84) y k ( 0 ) = ξ k ( 0 ∣ 0 ) = G k ( η 10 , η 20 ) , k = 1 , 2 .

Selecting α = 1 5 , η 10 = η 20 = − 100 yields the results depicted in Figure 4. Note that not only the deterministic, but also the stochastic solutions have fixed limits as t → ± ∞ . As can be seen from the solution expressions (83), the influence of the Wiener process is decreasing as t → + ∞ , since the rate of change of the deterministic part of the exponent outpaces the Wiener process.

$Figure 4 Solutions to ODEs (82) ( y 1 {y}_{1} and y 2 {y}_{2} depicted by dark red and dark blue lines respectively) and SDEs (81) ( ξ 1 {\xi }_{1} and ξ 2 {\xi }_{2} depicted by light red and light blue lines respectively) for different realizations of the Wiener process ω ( t ) \omega \left(t) . The parameter values are α = 1 5 \alpha =\frac{1}{5} , η 10 = η 20 = ‒ 100 {\eta }_{10}={\eta }_{20}=‒100 , y 1 ( 0 ) = 2.036 {y}_{1}\left(0)=2.036 , and y 2 ( 0 ) = 1.001 {y}_{2}\left(0)=1.001 .$

Figure 4

Solutions to ODEs (82) ( y 1 and y 2 depicted by dark red and dark blue lines respectively) and SDEs (81) ( ξ 1 and ξ 2 depicted by light red and light blue lines respectively) for different realizations of the Wiener process ω ( t ) . The parameter values are α = 1 5 , η 10 = η 20 = ‒ 100 , y 1 ( 0 ) = 2.036 , and y 2 ( 0 ) = 1.001 .

4.2 Comparison to the randomization approach

The presented technique can be compared to the most straightforward approach to introducing stochasticity into a system of ODE: randomization. In this approach, the ODE system (82) is integrated numerically, but a random number from a Gaussian distribution is added to the system function at each step.

Consider the system (82), a scaling variable ε > 0 , and two samples θ 1 ( k ) , θ 2 ( k ) , … , θ n ( k ) , k = 1 , 2 , of Gaussian random variables with a mean of 0 and a variance of 1. The randomization process can be described as follows: a numerical integration algorithm is used to integrate (82), but at each integration step, the ODE system is modified to have the following expression:

(85) d y k ( t ) d t = P k ( y 1 ( t ) , y 2 ( t ) ) + ε θ n ( k ) , k = 1 , 2 ,

where n is the number of integration step. This process yields randomized solutions ξ ^ 1 and ξ ^ 2 to system (82), as shown in Figure 5.

$Figure 5 Randomized solutions ξ ^ 1 {\widehat{\xi }}_{1} and ξ ^ 2 {\widehat{\xi }}_{2} of (82) ( y 1 {y}_{1} and y 2 {y}_{2} depicted by dark red and dark blue lines, respectively) depicted by light red and light blue lines, respectively. The initial conditions are set to y 1 ( 0 ) = 2.036 {y}_{1}\left(0)=2.036 and y 2 ( 0 ) = 1.001 {y}_{2}\left(0)=1.001 , as shown in Figure 4. Randomization was repeated 100 times.$

Figure 5

Randomized solutions ξ ^ 1 and ξ ^ 2 of (82) ( y 1 and y 2 depicted by dark red and dark blue lines, respectively) depicted by light red and light blue lines, respectively. The initial conditions are set to y 1 ( 0 ) = 2.036 and y 2 ( 0 ) = 1.001 , as shown in Figure 4. Randomization was repeated 100 times.

As can be seen from Figure 5, some randomized solutions diverge. In fact, there are three typical cases for randomized solutions, as shown in Figure 6. In Figure 6(a), the solution ξ ^ 2 diverges, while ξ ^ 1 increases past the maximum of the non-stochastic solution y 1 . In Figure 6(b), randomized solutions are bounded within the range of values of non-stochastic solutions y 1 and y 2 . In Figure 6(c), the solution ξ ^ 2 is destabilized initially but is then pushed back from divergence by a subsequent random effect.

$Figure 6 Typical cases of randomized solutions ξ ^ 1 {\widehat{\xi }}_{1} and ξ ^ 2 {\widehat{\xi }}_{2} of (82) ( y 1 {y}_{1} and y 2 {y}_{2} depicted by dark red and dark blue lines, respectively) depicted by light red and light blue lines, respectively. The initial conditions are set to y 1 ( 0 ) = 2.036 {y}_{1}\left(0)=2.036 and y 2 ( 0 ) = 1.001 {y}_{2}\left(0)=1.001 , as shown in Figure 4.$

Figure 6

Typical cases of randomized solutions ξ ^ 1 and ξ ^ 2 of (82) ( y 1 and y 2 depicted by dark red and dark blue lines, respectively) depicted by light red and light blue lines, respectively. The initial conditions are set to y 1 ( 0 ) = 2.036 and y 2 ( 0 ) = 1.001 , as shown in Figure 4.

There are notable differences between the approach presented in this article and randomization. First, the solutions to the SDE obtained via the presented approach cannot diverge, and the degree of randomness in them is related to the derivative value of the non-stochastic solution. This is visible in Figure 4, where the solutions to the SDE are near the solutions to the ODE at small and large values of t . The same is not true for the randomization approach. In fact, many solutions diverge quite rapidly (Figure 5).

5 Concluding remarks

Extending the work on the one-dimensional case in the study by Navickas et al. [15], a scheme for the stochastization of a system of ODEs is presented in this article. Note that the results of this stochastization scheme are completely distinct from randomization: instead of random effects introduced during numerical integration, the presented scheme allows us to observe analytical solutions for both the non-stochastic ODEs and obtained stochastic SDEs. Furthermore, the introduced parameter α allows control over the level of stochastization of a particular system or solution: as the parameter tends to zero, stochastic solutions coincide with deterministic solutions, while the system of SDEs becomes a system of ODEs.

As shown in the Section 3, the presented scheme can also be used to construct both the system of ODEs and SDEs given the diffusion functions that govern the stochastic part of the model. This has a potential impact on applications, as supplying the diffusion functions from experimental results could be used to construct a model of differential equations describing the process. Further investigation of these possibilities is a definite objective of future research.

Another objective of future research pertains to the generalization of the presented stochastization techniques into n dimensions ( n > 2 ). While the approach presented here is robust for two dimensions, higher-order systems would require a significant adjustment. In fact, investigating whether dimension n has an impact on the required derivation techniques is one of the first goals of further studies.

Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors declare that there are no conflicts of interest.

References

[1] J. Blath, P. Imkeller, and S. Rœlly, Surveys in Stochastic Processes, Vol. 4, European Mathematical Society, Switzerland 2011. 10.4171/072Search in Google Scholar

[2] C. Huang, F. Cattani, and F. Galvanin, An optimal experimental design strategy for improving parameter estimation in stochastic models, Comput. Chem. Eng. 170 (2023), 108133. 10.1016/j.compchemeng.2023.108133Search in Google Scholar

[3] T. Khan, R. Ullah, B. Al Alwan, Y. El-Khatib, and G. Zaman, Correlated stochastic epidemic model for the dynamics of SARS-CoV-2 with vaccination, Sci. Rep. 12 (2022), no. 1, 16105. 10.1038/s41598-022-20059-0Search in Google Scholar PubMed PubMed Central

[4] J. R Mihaljevic, S. Borkovec, S. Ratnavale, T. D. Hocking, K. E. Banister, J. E. Eppinger, et al. SPARSEMODr: Rapidly simulate spatially explicit and stochastic models of COVID-19 and other infectious diseases, Biol. Methods Protoc. 7 (2022), no. 1, bpac022. 10.1093/biomethods/bpac022Search in Google Scholar PubMed PubMed Central

[5] A. Omame, M. Abbas, and A. Din, Global asymptotic stability, extinction and ergodic stationary distribution in a stochastic model for dual variants of SARS-CoV-2, Math. Comput. Simulation 204 (2023), 302–336. 10.1016/j.matcom.2022.08.012Search in Google Scholar PubMed PubMed Central

[6] Z. Shi and D. Jiang, Environmental variability in a stochastic HIV infection model, Commun. Nonlinear Sci. Numer. Simul. 120 (2023), 107201. 10.1016/j.cnsns.2023.107201Search in Google Scholar

[7] D. H. Serrano, J. Villarroel, J. Hernández-Serrano, and A. Tocino, Stochastic simplicial contagion model, Chaos Solitons Fractals 167 (2023), 113008. 10.1016/j.chaos.2022.113008Search in Google Scholar

[8] Y. Haiyang, W. Haiyan, Z. Zhichen, X. Yong, and J. Kurths, A stochastic nonlinear differential propagation model for underwater acoustic propagation: Theory and solution, Chaos Solitons Fractals 150 (2021), 111105. 10.1016/j.chaos.2021.111105Search in Google Scholar

[9] H. Verdejo, A. Awerkin, W. Kliemann, and C. Becker, Modelling uncertainties in electrical power systems with stochastic differential equations, Int. J. Electr. Power Energy Syst. 113 (2019), 322–332. 10.1016/j.ijepes.2019.05.054Search in Google Scholar

[10] G. I. Schuëller and H. J. Pradlwarter, Advances in stochastic structural dynamics under the perspective of reliability estimation, EURODYN 99 (1999), 267–272. Search in Google Scholar

[11] R. G Ghanem and P. D. Spanos, Stochastic Finite Elements: A Spectral Approach, Courier Corporation, N.Y., 2003. Search in Google Scholar

[12] M. Eldred and J. Burkardt, Comparison of non-intrusive polynomial chaos and stochastic collocation methods for uncertainty quantification, 47th AIAA Aerospace Sciences Meeting Including The New Horizons Forum And Aerospace Exposition, 2009. 10.2514/6.2009-976Search in Google Scholar

[13] I. Babussska, F. Nobile, and R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal. 45 (2007), no. 3, 1005–1034. 10.1137/050645142Search in Google Scholar

[14] B. Besselink, U. Tabak, A. Lutowska, N. Van de Wouw, H. Nijmeijer, D. J Rixen, et al., A comparison of model reduction techniques from structural dynamics, numerical mathematics and systems and control, J. Sound Vib. 332 (2013), no. 19, 4403–4422. 10.1016/j.jsv.2013.03.025Search in Google Scholar

[15] Z. Navickas, I. Timofejeva, T. Telksnys, R. Marcinkevicius, and M. Ragulskis, Construction of special soliton solutions to the stochastic Riccati equation, Open Math. 20 (2022), no. 1, 829–844. 10.1515/math-2022-0051Search in Google Scholar

[16] R. Durrett, Stochastic Calculus: A Practical Introduction, CRC Press, Florida, 2018. 10.1201/9780203738283Search in Google Scholar

[17] B. Øksendal, Stochastic Differential Equations, Springer, Berlin, 2003, pp. 65–84. 10.1007/978-3-642-14394-6_5Search in Google Scholar

[18] L. C. Evans, Partial Differential Equations: Graduate Studies in Mathematics, Vol. 19, American Mathematical Society, Rhode Island, 2009. Search in Google Scholar

[19] I. Mezo, The Lambert W Function: Its Generalizations and Applications, CRC Press, Florida, 2022. 10.1201/9781003168102Search in Google Scholar

Received: 2023-03-09

Revised: 2023-08-28

Accepted: 2023-09-26

Published Online: 2023-11-11

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/math-2023-0136

Keywords for this article

stochastic calculus; differential equation; Itô equation

Creative Commons

BY 4.0