Convergence and Stability of the Split-Step 𝜽-Milstein Method for Stochastic Delay Hopfield Neural Networks

Abstract

A new splitting method designed for the numerical solutions of stochastic delay Hopfield neural networks is introduced and analysed. Under Lipschitz and linear growth conditions, this split-step θ-Milstein method is proved to have a strong convergence of order 1 in mean-square sense, which is higher than that of existing split-step θ-method. Further, mean-square stability of the proposed method is investigated. Numerical experiments and comparisons with existing methods illustrate the computational efficiency of our method.

1. Introduction

Hopfield neural networks, which originated with Hopfield in the 1980s [1], have been successfully applied in many areas such as combinatorial optimization [2, 3], signal processing [4], and pattern recognition [5, 6]. In the last decade, neural networks in the presence of signal transmission delay and stochastic perturbations, also named as stochastic delay Hopfield neural networks (SDHNNs), have gained considerable research interest (see, e.g., [7–9] and the references therein). It is noticed that, so far, most works on SDHNNs focus mainly on the stability analysis of the analytical solutions, including mean-square exponential stability [7], global asymptotic stability [9], and so forth. Not only simulation is an important tool to explore interesting dynamics of kinds of Hopfield neural networks (HNNs) (see, e.g., [10] and the references therein), but also parameter estimation in dynamical systems based on HNNs (see, e.g., [11]) needs to solve HNNs numerically. Moreover, because most of SDHNNs do not have explicit solutions, the numerical analysis of SDHNNs recently stirred some initial research attention. For example, Li et al. [12] investigated the exponential stability of the Euler method and the semi-implicit Euler method for SDHNNs. Rathinasamy [13] introduced a split-step 𝜃-method (SST) for SDHNNs and analysed the mean-square stability of this method, and the SST is only given for the commensurable delay case. To the best of our current knowledge, the authors mainly discussed the stability of numerical solutions for stochastic Hopfield neural networks with discrete time delays but skipped the details of convergence analysis.

The split-step Euler method for stochastic differential equations (SDEs) was proposed by Higham et al. [14], further, the splitting Euler-type algorithms have been derived for stochastic delay differential equations (SDDEs) [15, 16]. In this paper, we will present a splitting method with higher order convergence for SDHNNs. To be specific, we will go into detail about the convergence analysis and comparing the stability with split-step 𝜃-method given in [13].

The rest of this paper is organized as follows. In Section 2, we recall the stochastic delay neural networks model and present a split-step 𝜃-Milstein method. In Section 3, we derive the convergence results of the split-step 𝜃-Milstein method for the model. In Section 4, the numerical stability analysis is performed. In Section 5, some numerical examples are given to confirm the theory. In the last Section, we draw some conclusions.

2. Model and the Split-Step 𝜽-Milstein Method

2.1. Model

Consider the stochastic delay Hopfield neural networks of the form d𝐱(𝑡)=[−𝐶𝐱(𝑡)+𝐴𝐟(𝐱(𝑡))+𝐵𝐠(𝐳(𝑡))]dt+Φ(𝐱(𝑡))d𝑊(𝑡), (1) where 𝐱(𝑡)=[𝑥₁(𝑡),…,𝑥_𝑛(𝑡)]^𝑇∈ℝ^𝑛 is the state vector associated with the 𝑛 neurons, 𝐳(𝑡)=[𝑥₁(𝑡−𝜏₁),…,𝑥_𝑛(𝑡−𝜏_𝑛)]^𝑇, the diagonal matrix 𝐶=diag(𝑐₁,𝑐₂,…,𝑐_𝑛) has positive entries, and 𝑐_𝑖 represents the rate at which the 𝑖th unit will reset its potential to the resting state in isolation when discounted from the network and the external stochastic perturbation. The matrices 𝐴=(𝑎_𝑖𝑗)_𝑛×𝑛 and 𝐵=(𝑏_𝑖𝑗)_𝑛×𝑛 are the connection weight matrix and the discretely delayed connection weight matrix, respectively. Furthermore, the vector functions 𝐟(𝐱(𝑡))=[𝑓₁(𝑥₁(𝑡)),…,𝑓_𝑛(𝑥_𝑛(𝑡))]^𝑇 and 𝐠(𝐳(𝑡))=[𝑔₁(𝑥₁(𝑡−𝜏₁)),…,𝑔_𝑛(𝑥_𝑛(𝑡−𝜏_𝑛))]^𝑇 denote the neuron activation functions with the conditions 𝑓_𝑖(0)=0,𝑔_𝑖(0)=0 for all positive 𝜏_𝑗.

On the initial segment [−𝜏,0] the state vector satisfies 𝐱(𝑡)=𝜓(𝑡), where 𝜓(𝑡)=[𝜓₁(𝑡),…,𝜓_𝑛(𝑡)]^𝑇 is a given function in 𝐶([−𝜏,0],ℝ^𝑛) and 𝜏 stands for max_{1≤𝑖≤𝑛}{𝜏_𝑖}.

Moreover, Φ(𝐱(𝑡))=diag(𝜙₁(𝑥₁(𝑡)),…,𝜙_𝑛(𝑥_𝑛(𝑡))) is a diagonal matrix with 𝜙_𝑖(0)=0 and 𝑊(𝑡)=[𝑊₁(𝑡),…,𝑊_𝑛(𝑡)]^𝑇∈ℝ^𝑛 is an 𝑛-dimensional Wiener process defined on the complete probability space (Ω,ℱ,{ℱ_𝑡}_𝑡≥0,ℙ) with a filtration {ℱ_𝑡}_𝑡≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while ℱ₀ contains all ℙ-null sets).

Let 𝑓_𝑖 and 𝑔_𝑖 be functions in 𝐶²(𝐷;ℝ)⋂ℒ²([0,𝑇];ℝ) and 𝜙_𝑖 be in 𝐶¹(𝐷;ℝ)∩ℒ²([0,𝑇];ℝ). Here 𝐶^𝑙(𝐷;ℝ) denotes the family of continuously 𝑙-times differentiable real-valued function defined on 𝐷, while ℒ^𝑙([0,𝑇];ℝ) denotes the family of all real-valued measurable {ℱ_𝑡}-adapted stochastic processes {𝑓(𝑡)}_{𝑡∈[0,𝑇]} such that ∫^𝑇₀|𝑓(𝑡)|^𝑙𝑑𝑡<+∞.

2.2. Numerical Scheme

We define the mesh with a uniform step-size ℎ(0<ℎ<1) on the interval [0,𝑇]; that is, 𝑡_𝑘=𝑘⋅ℎ(𝑘=0,1,…,𝐾) and 𝑇=𝐾ℎ.

Let Δ𝑊^𝑘_𝑖=𝑊_𝑖(𝑡_𝑘+1)−𝑊_𝑖(𝑡_𝑘) denote the increment of the Wiener process. The split-step 𝜃-Milstein (SSTM) scheme for the solution of SDEs (1) is given by 𝑌^𝑘=𝑋^𝑘+−(1−𝜃)𝐶𝑋^𝑘−𝜃𝐶𝑌^𝑘+𝐴𝐟𝑋^𝑘+𝐵𝐠𝑍^𝑘𝑋ℎ,_𝑘+1=𝑌^𝑘+Φ𝑌^𝑘Δ𝑊^𝑘+12Φ𝑌^𝑘Δ𝑊^𝑘∘Δ𝑊^𝑘−𝟙ℎ, (2) where the merging parameter 𝜃 satisfies 0≤𝜃≤1, 𝑋^𝑘=[𝑋^k₁,…,𝑋^𝑘_𝑛]^𝑇 is an approximation to 𝐱(𝑡_𝑘), and for 1≤𝑞_𝑗∈ℤ⁺𝑍^𝑘_𝑗=𝜓_𝑗𝑡_𝑘−𝜏_𝑗,𝑡_𝑘−𝜏_𝑗𝑋≤0,_{𝑘−𝑞𝑗𝑗},0<𝑡_𝑘−𝜏_𝑗∈𝑡_{𝑘−𝑞𝑗},𝑡_{𝑘−𝑞𝑗+1}. (3) Moreover, we adopt the symbols Φ(𝑌^𝑘) = diag(𝜙^′₁(𝑌^𝑘₁)𝜙₁(𝑌^𝑘₁),…,𝜙^′_𝑛(𝑌^𝑘_𝑛)𝜙_𝑛(𝑌^𝑘_𝑛)) and Δ𝑊^𝑘 = [Δ𝑊^𝑘₁,…,Δ𝑊^𝑘_𝑛]^𝑇, the Hadamard product Δ𝑊^𝑘∘Δ𝑊^𝑘 means [(Δ𝑊^𝑘₁)²,…,(Δ𝑊^𝑘_𝑛)²]^𝑇, and 𝟙=[1,…,1]^𝑇∈ℝ^𝑛. When 𝑡_𝑘≤0, we define 𝑋^𝑘=𝜓(𝑡_𝑘).

Then scheme (2) can be written in equivalent form as𝑌^𝑘=𝑋^𝑘−(𝐼+𝜃𝐶ℎ)⁻¹𝐶ℎ𝑋^𝑘+(𝐼+𝜃𝐶ℎ)⁻¹ℎ𝐴𝐟𝑋^𝑘+𝐵𝐠𝑍^𝑘,𝑋^𝑘+1=𝑌^𝑘+Φ𝑌^𝑘Δ𝑊^𝑘+12Φ𝑌^𝑘Δ𝑊^𝑘∘Δ𝑊^𝑘−𝟙ℎ. (4a)(4b)Substituting (4a) into (4b), we have a stochastic explicit single-step method with an increment function Λ(𝜉,𝜂,ℎ,Δ𝑊^𝑘)∈ℝ^𝑛; that is, 𝑋^𝑘+1=𝑋^𝑘+Λ𝑋^𝑘,𝑍^𝑘,ℎ,Δ𝑊^𝑘. (5)

3. Order and Convergence Results for SSTM

In this section we consider the global error of SSTM (2) as applied to SDHNNs (1) with initial condition. In what follows, ‖⋅‖ denotes Euclidean norm in ℝ^𝑛.

For convergence purpose we make the following standard assumptions.

Assumption 1. Assume that 𝐟, 𝐠, Φ, and Φ satisfy the Lipschitz condition 𝑓_𝑖(𝑎)−𝑓_𝑖(𝑏)≤𝛼_𝑖||,𝑔𝑎−𝑏_𝑖(𝑎)−𝑔_𝑖(𝑏)≤𝛽_𝑖||,(𝑎−𝑏𝑎,𝑏∈ℝand[𝜙𝑡∈0,𝑇])_i(𝑎)−𝜙_𝑖(𝑏)≤𝛾_𝑖||,𝜙𝑎−𝑏_𝑖𝜙^′_𝑖(𝑎)−𝜙_𝑖𝜙^′_𝑖(𝑏)≤𝜅_𝑖|𝑎−𝑏|, (6) for every 𝑖 and the linear growth condition ‖𝐟(𝐱)‖²∨‖𝐠(𝐱)‖²∨‖Φ(𝐱)⋅𝟙‖²∨Φ(𝐱)⋅𝟙²≤𝐿^∼1+‖𝐱‖²,∀𝐱∈ℝ^𝑛, (7) where 𝐿^∼ is a positive constant and ∨ is the maximal operator. We also define 𝐿_𝑖 as 𝐿_𝑖=max{𝛼_𝑖,𝛽_𝑖,𝛾_𝑖,𝜅_𝑖}  (𝑖=1,2,…,𝑛).

We also need the following assumption on the initial condition.

Assumption 2. Assume that the initial function 𝜓(𝑡) is Lipschitz continuous from [−𝜏,0] to ℝ^𝑛, that is, there is a positive constant 𝐿_𝜓 satisfying 𝜓(𝑡)−𝜓(𝑠)≤𝐿_𝜓(𝑡−𝑠)if−𝜏≤𝑠<𝑡<0. (8)

Now we give the definition of local and global errors.

Definition 1. Let 𝐱(𝑡) denote the exact solution of (1). The local approximate solution 𝐱^∼(𝑡_𝑘+1) starting from 𝐱(𝑡_𝑘) by SSTM (2) given by 𝐱^∼𝑡_𝑘+1∶=𝐱𝑡_𝑘+Λ𝐱𝑡_𝑘,𝑍^∼𝑡_𝑘,ℎ,Δ𝑊^𝑘, (9) where 𝑍^∼(𝑡_𝑘) denotes the evaluation of (3) using the exact solution, yields the difference 𝛿^𝑘+1∶=𝐱𝑡_𝑘+1−𝐱^∼𝑡_𝑘+1. (10) Then the local error of SSTM is defined by ‖𝛿^𝑘+1‖, whereas its global error means ‖𝜖^𝑘‖ where 𝜖^𝑘∶=𝐱(𝑡_𝑘)−𝑋^𝑘.

Definition 2. If the global error satisfies 𝐸𝜖^𝑘2≤Γℎ^2𝑝∀ℎ∈0,ℎ₀ (11) with positive constants ℎ₀ and Γ and a finite 𝑝, then we say that the order of mean-square convergence accuracy of the method is 𝑝. Here 𝐸 is the expectation with respect to ℙ.

We then give the following lemmas that are useful in deriving the convergence results.

Lemma 3 (see also [17]). Let the linear growth condition (7) hold, and the initial function 𝜓(𝑡) is assumed to be ℱ₀-measurable and right continuous. And one puts 𝐸_𝜓∶=𝐸(sup_{−𝜏≤𝑡≤0}‖𝜓(𝑡)‖²)<∞. For any given positive 𝑇, there exist positive numbers Γ_𝜓 and Γ₂ such that the solution of (1) satisfies 𝐸sup_{−𝜏≤𝑠≤𝑇}‖𝐱(𝑠)‖²≤Γ_𝜓, (12) where the constant Γ_𝜓 is independent of step-size ℎ but dependent on 𝑇. Moreover, for any 0≤𝑠<𝑡≤𝑇,  𝑡−𝑠<1, the estimation 𝐸(‖𝐱(𝑡)−𝐱(𝑠)‖)²≤Γ₂(𝑡−𝑠) (13) holds.

The Jensen inequality derives 𝐸sup_{−𝜏≤𝑠≤𝑇}‖𝐱(𝑠)‖≤Γ_𝜓 (14) from (12).

Lemma 4. For 𝑠∈[𝑡_𝑘,𝑡_𝑘+ℎ], one has 𝐸𝐳(𝑠)−𝑍^∼𝑡_𝑘²≤Γ_𝜏ℎ. (15) Here the constant Γ_𝜏 is independent of step-size ℎ.

Proof. If 𝑡_𝑘−𝜏_𝑗≤0 and 𝑠−𝜏_𝑗≤0, under Assumption 2 we have 𝐸𝑥_𝑗𝑠−𝜏_𝑗−𝑍^∼_𝑗𝑡_𝑘²=𝐸𝜓_𝑗𝑠−𝜏_𝑗−𝜓_𝑗𝑡_𝑘−𝜏_𝑗²≤𝐿_𝜓𝑠−𝑡_𝑘²≤𝐿_𝜓ℎ². (16) If 𝑡_𝑘−𝜏_𝑗≤0 and 𝑠−𝜏_𝑗>0, with (13) we obtain 𝐸𝑥_𝑗𝑠−𝜏_𝑗−𝑍^∼_𝑗𝑡_𝑘²=𝐸𝑥_𝑗𝑠−𝜏_𝑗−𝜓_𝑗𝑡_𝑘−𝜏_𝑗²≤2𝐸𝑥_𝑗𝑠−𝜏_𝑗−𝑥_𝑗(0)²+2𝐸𝜓_𝑗(0)−𝜓_𝑗𝑡_𝑘−𝜏_𝑗²≤2Γ₂+𝐿_𝜓ℎℎ. (17) If 𝑡_𝑘−𝜏_𝑗>0 and 𝑠−𝜏_𝑗>0, we assume 𝑡_𝑘−𝜏_𝑗∈[𝑡_{𝑘−𝑞𝑗},𝑡_{𝑘−𝑞𝑗+1}) without loss of generality. Hence, 𝐸𝑥_𝑗𝑠−𝜏_𝑗−𝑍^∼_𝑗𝑡_𝑘²=𝐸𝑥_𝑗𝑠−𝜏_𝑗−𝑥_𝑗𝑡_{𝑘−𝑞𝑗}²≤2𝐸𝑥_𝑗𝑠−𝜏_𝑗−𝑥_𝑗𝑡_𝑘−𝜏_𝑗²+2𝐸𝑥_𝑗𝑡_{𝑘−𝑞𝑗}−𝑥_𝑗𝑡_𝑘−𝜏_𝑗²≤4Γ₂ℎ (18) by using inequality (13).

Lemma 5. Let 𝐱(𝑡) denote the exact solution of (1). One assumes conditions (6) and (7). Then for the local intermediate value 𝐲(𝑡_𝑘)∶=𝐱(𝑡_𝑘)−(𝐼+𝜃𝐶ℎ)⁻¹𝐶ℎ𝐱(𝑡_𝑘)+(𝐼+𝜃𝐶ℎ)⁻¹ℎ[𝐴𝐟(𝐱(𝑡_𝑘))+𝐵𝐠(𝑍^∼(𝑡_𝑘))], one has the estimation 𝐸𝐱(𝑡_𝑘)−𝐲(𝑡_𝑘)²≤Γ₃ℎ². (19)

Proof. The difference between the 𝑖th components of 𝐱(𝑡_𝑘) and 𝐲(𝑡_𝑘) leads to 𝑥_𝑖𝑡_𝑘−𝑦_𝑖𝑡_𝑘²=ℎ²1+𝑐_𝑖𝜃ℎ²𝑐_𝑖𝑥_𝑖𝑡_𝑘−^𝑛_𝑗=1𝑎_𝑖𝑗𝑓_𝑗𝑥_𝑗𝑡_𝑘−^𝑛_𝑗=1𝑏_𝑖𝑗𝑔𝑍_𝑗^∼_𝑗𝑡_𝑘², (20) whose expectation, together with 𝑐_𝑖>0, 𝑓_𝑖(0)=𝑔_𝑖(0)=0, the Lipschitz condition (6), and the estimation (12), gives 𝐸𝐱(𝑡_𝑘)−𝐲(𝑡_𝑘)²≤^𝑛_𝑖=13ℎ²1+𝑐_𝑖𝜃ℎ²×𝑐²_𝑖𝐸𝑥_𝑖𝑡_𝑘²+𝑛^𝑛_𝑗=1𝑎²_𝑖𝑗𝛼²_𝑗𝐸𝑥_𝑗𝑡_𝑘²+𝑛^𝑛_𝑗=1𝑏²_𝑖𝑗𝛽²_𝑗𝐸𝑍^∼_𝑗𝑡_𝑘²≤3Γ_𝜓^𝑛_𝑖=1𝑐²_𝑖+𝑛^𝑛_𝑗=1𝑎²_𝑖𝑗𝛼²_𝑗+𝑛^𝑛_𝑗=1𝑏²_𝑖𝑗𝛽²_𝑗ℎ². (21)

Now we discuss local error estimates.

Theorem 6. When one assumes Assumptions 1 and 2 and the conditions of Lemma 3, there exist positive constants Γ₀ and Γ₁, such that max_{0≤𝑘≤𝐾−1}𝐸𝛿^𝑘+1≤Γ₀ℎ²,max_{0≤𝑘≤𝐾−1}𝐸𝛿^𝑘+12≤Γ₁ℎ³ (22)(23) as ℎ↓0.

Proof. The Itô integral form of the 𝑖th component of (1) on [𝑡_𝑘,𝑡] implies 𝑥_𝑖(𝑡)−𝑥_𝑖𝑡_𝑘=^𝑡_𝑡𝑘−𝑐_𝑖𝑥_𝑖(𝑠)+^𝑛_𝑗=1𝑎_𝑖𝑗𝑓_𝑗𝑥_𝑗(𝑠+_𝑛𝑗=1)𝑏_𝑖𝑗𝑔_𝑗𝑥_𝑗𝑠−𝜏_𝑗d𝑠+^𝑡_𝑡𝑘𝜙_𝑖𝑥_𝑖(𝑠)d𝑊_𝑖(𝑠). (24) By utilizing the previous identity, the 𝑖th component of the difference 𝛿^𝑘+1 introduced in Definition 1 can be calculated as 𝛿_𝑖^𝑘+1=𝑥_𝑖𝑡_𝑘+1−𝑥_𝑖𝑡_𝑘+𝑐_𝑖ℎ1+𝑐_𝑖𝑥𝜃ℎ^𝑖𝑡_𝑘−ℎ1+𝑐_𝑖𝜃ℎ^𝑛_𝑗=1𝑎_𝑖𝑗𝑓_𝑗𝑥_𝑗𝑡_𝑘+^𝑛_𝑗=1𝑏_𝑖𝑗𝑔𝑍_𝑗^∼_𝑗𝑡_𝑘−𝜙_𝑖𝑦_𝑖𝑡_𝑘Δ𝑊^𝑘_𝑖−𝜙^′_𝑖𝑦_𝑖𝑡_𝑘×𝜙_𝑖𝑦_𝑖𝑡_𝑘Δ𝑊^𝑘_𝑖²−ℎ2=^𝑡𝑘+1_𝑡𝑘𝑐_𝑖𝑥_𝑖𝑡_𝑘−𝑥_𝑖(𝑠)d𝑐𝑠−^2𝑖ℎ²𝜃1+𝑐_𝑖𝑥𝜃ℎ^𝑖𝑡_𝑘+^𝑡𝑘+1_𝑡𝑘^𝑛_𝑗=1𝑎_𝑖𝑗𝑓_𝑗𝑥_𝑗(𝑠)−𝑓_𝑗𝑥_𝑗𝑡_𝑘+^𝑛_𝑗=1𝑏_𝑖𝑗𝑔_𝑗𝑥_𝑗𝑠−𝜏_𝑗−𝑔𝑍_𝑗^∼_𝑗𝑡_𝑘d𝑠+𝑐_𝑖𝜃ℎ²1+𝑐_𝑖𝜃ℎ^𝑛_𝑗=1𝑎_𝑖𝑗𝑓_𝑗𝑥_𝑗𝑡_𝑘+^𝑛_𝑗=1𝑏_𝑖𝑗𝑔𝑍_𝑗^∼_𝑗𝑡_𝑘+^𝑡𝑘+1_𝑡𝑘𝜙_𝑖𝑥_𝑖(𝑠)d𝑊_𝑖(𝑠)−^𝑡𝑘+1_𝑡𝑘𝜙_𝑖𝑦_𝑖𝑡_𝑘d𝑊_𝑖(𝑠)−^𝑡𝑘+1_𝑡𝑘^𝑠_𝑡𝑘𝜙^′_𝑖𝑦_𝑖𝑡_𝑘𝜙_𝑖𝑦_𝑖𝑡_𝑘d𝑊_𝑖(𝑟)d𝑊_𝑖(𝑠), (25) where Δ𝑊^𝑘_𝑖=∫^𝑡𝑘+1_𝑡𝑘d𝑊_𝑖(𝑠) and ((Δ𝑊^𝑘_i)²−ℎ)/2=∫^𝑡𝑘+1_𝑡𝑘∫^𝑠_𝑡𝑘d𝑊_𝑖(𝑟)d𝑊_𝑖(𝑠).
Taking expectations of both sides of (25), 𝐸𝛿_𝑖^𝑘+1=𝐸^𝑡𝑘+1_𝑡𝑘𝑐_𝑖^𝑠_𝑡𝑘−𝜇_𝑖𝑥_𝑖(𝑟)d𝑟d𝑠−𝑐²_𝑖ℎ²𝜃1+𝑐_𝑖𝐸𝑥𝜃ℎ^𝑖𝑡_𝑘+𝐸^𝑡𝑘+1_𝑡𝑘^𝑛_𝑗=1𝑎_𝑖𝑗^𝑠_𝑡𝑘𝑓^′_𝑗𝑥_𝑗(𝑟)𝜇_𝑗𝑥_𝑗(𝑟)+12𝑓_𝑗^′′𝑥_𝑗(𝑟)𝜙²_𝑗𝑥_𝑗(𝑟)d𝑟d𝑠+𝐸^𝑡𝑘+1_𝑡𝑘^𝑛_𝑗=1𝑏_𝑖𝑗^{𝑠−𝜏𝑗}_{𝑡𝑗𝑘−𝑞}𝑔^′_𝑗𝑥_𝑗(𝑟)𝜇_𝑗𝑥_𝑗(𝑟+1)2𝑔_𝑗^′′𝑥_𝑗(𝑟)𝜙²_𝑗𝑥_𝑗(𝑟)d𝑟d𝑠+𝑐_𝑖𝜃ℎ²1+𝑐_𝑖𝜃ℎ^𝑛_𝑗=1𝑎_𝑖𝑗𝐸𝑓_𝑗𝑥_𝑗𝑡_𝑘+^𝑛_𝑗=1𝑏_𝑖𝑗𝐸𝑔𝑍_𝑗^∼_𝑗𝑡_𝑘 (26) by (24) and Itô formula, where 𝜇_𝑖(𝑥_𝑖(𝑟))  := −𝑐_𝑖𝑥_𝑖(𝑟)+∑^𝑛_𝑗=1𝑎_𝑖𝑗𝑓_𝑗(𝑥_𝑗(𝑟))+∑^𝑛_𝑗=1𝑏_𝑖𝑗𝑔_𝑗(𝑥_𝑗(𝑟−𝜏_𝑗)). Under conditions of this theorem, we have |𝐸(𝛿_𝑖^𝑘+1)|≤Γ^∼₀ℎ² by (14), Jensen inequality |𝐸(𝑋)|≤𝐸(|𝑋|), triangle inequality, and properties of definite integral. Then we have max_{0≤𝑘≤𝐾−1}‖𝐸(𝛿^𝑘+1)‖≤Γ₀ℎ² from the relation between |𝐸(𝛿_𝑖^𝑘+1)| and ‖𝐸(𝛿^𝑘+1)‖.
Now we prove (23). By Itô formula, ^𝑡𝑘+1_𝑡𝑘𝜙_𝑖𝑥_𝑖(𝑠)d𝑊_𝑖(𝑠)−^𝑡𝑘+1_𝑡𝑘𝜙_𝑖𝑦_𝑖𝑡_𝑘d𝑊_𝑖(𝑠)−^𝑡𝑘+1_𝑡𝑘^𝑠_𝑡𝑘𝜙^′_𝑖𝑦_𝑖𝑡_𝑘𝜙_𝑖𝑦_𝑖𝑡_𝑘d𝑊_𝑖(𝑟)d𝑊_𝑖(𝑠)=^𝑡𝑘+1_𝑡𝑘𝜙_𝑖𝑥_𝑖𝑡_𝑘−𝜙_𝑖𝑦_𝑖𝑡_𝑘d𝑊_𝑖(𝑠)+^𝑡𝑘+1_𝑡𝑘^𝑠_𝑡𝑘𝜙^′_𝑖𝑥_𝑖(𝑟)𝜇_𝑖𝑥_𝑖(𝑟)+12𝜙_𝑖^′′𝑥_𝑖(𝑟)𝜙_𝑖𝑥_𝑖(𝑟)²d𝑟d𝑊_𝑖(𝑠)+^𝑡𝑘+1_𝑡𝑘^𝑠_𝑡𝑘𝜙^′_𝑖𝑥_𝑖(𝑟)𝜙_𝑖𝑥_𝑖(𝑟)−𝜙^′_𝑖𝑦_𝑖𝑡_𝑘×𝜙_𝑖𝑦_𝑖𝑡_𝑘d𝑊_𝑖(𝑟)d𝑊_𝑖(𝑠). (27) From (25) and (27), we have 𝐸𝛿_𝑖^𝑘+12≤9𝑐_𝑖ℎ^𝑡𝑘+1_𝑡𝑘𝐸𝑥_𝑖𝑡_𝑘−𝑥_𝑖(𝑠)²d𝑠_{Cauchy−Schwartzinequality}+𝑛^𝑛_𝑗=1𝑎²_𝑖𝑗𝛼²_𝑗ℎ^𝑡𝑘+1_𝑡𝑘𝐸𝑥_𝑗(𝑠)−𝑥_𝑗𝑡_𝑘²d𝑠_{Cauchy−SchwartzandḦ‌olderinequality}+𝑛^𝑛_𝑗=1𝑏²_𝑖𝑗𝛽²_𝑗ℎ^𝑡𝑘+1_𝑡𝑘𝐸𝑥_𝑗𝑠−𝜏_𝑗−𝑍^∼_𝑗𝑡_𝑘²d𝑠_{Cauchy−SchwartzandḦ‌olderinequality}+𝑐⁴_𝑖ℎ⁴𝜃²1+𝑐_𝑖𝜃ℎ²𝐸𝑥_𝑖𝑡_𝑘²+𝑛𝑐²_𝑖ℎ⁴𝜃²1+𝑐_𝑖𝜃ℎ^2𝑛_𝑗=1𝑎²_𝑖𝑗𝐸𝑓_𝑗𝑥_𝑗𝑡_𝑘²+^𝑛_𝑗=1𝑏²_𝑖𝑗𝐸𝑔𝑍_𝑗^∼_𝑗𝑡_𝑘²+^𝑡𝑘+1_𝑡𝑘𝛾²_𝑗𝐸𝑥_𝑖𝑡_𝑘−𝑦_𝑖𝑡_𝑘²d𝑠_ItôIsometry+^𝑡𝑘+1_𝑡𝑘(𝑠−𝑡_𝑘)^𝑠_𝑡𝑘𝐸𝜙^′_𝑖𝑥_𝑖(𝑟)𝜇_𝑖𝑥_𝑖(𝑟)+(1/2)𝜙_𝑖^′′𝑥_𝑖(𝑟)𝜙_𝑖𝑥_𝑖(𝑟)²²d𝑟d𝑠_{Cauchy−SchwartzinequalityandItôIsometry}+^𝑡𝑘+1_𝑡𝑘^𝑠_𝑡𝑘𝐸𝜙^′_𝑖𝑥_𝑖(𝑟)𝜙_i𝑥_𝑖(𝑟)−𝜙^′_𝑖𝑦_𝑖𝑡_𝑘𝜙_𝑖𝑦_𝑖𝑡_𝑘²d𝑟d𝑠_ItôIsometry≤9ℎ𝑐_𝑖+𝑛^𝑛_𝑗=1𝑎²_𝑖𝑗𝛼²_𝑗^𝑡𝑘+1_𝑡𝑘Γ₂ℎ₍₁₃₎d𝑠+ℎ𝑛^𝑛_𝑗=1𝑏²_𝑖𝑗𝛽²_𝑗^𝑡𝑘+1_𝑡𝑘Γ_𝜏ℎ₍₁₅₎d𝑠+ℎ⁴𝑐²_𝑖𝜃²1+𝑐_𝑖𝜃ℎ²Γ_𝜓𝑐²_𝑖+𝑛^𝑛_𝑗=1𝑎²_𝑖𝑗𝛼²_j+𝑛^𝑛_𝑗=1𝑏²_𝑖𝑗𝛽²_{𝑗(6)and(12)}+^𝑡𝑘+1_𝑡𝑘𝛾²_𝑗Γ₃ℎ²₍₁₉₎d𝑠+^𝑡𝑘+1_𝑡𝑘Γ₄𝑠−𝑡_𝑘²_{momentbounded}d𝑠+^𝑡𝑘+1_𝑡𝑘^𝑠_𝑡𝑘Γ₅ℎd𝑟d𝑠_{TriangleInequality,(6),(13),(19)}≤Γ^∼₁ℎ³. (28) Finally, it is easy to prove 𝐸(‖𝛿^𝑘+1‖²)≤Γ₁ℎ³.

Thanks to Theorem 1 in [18], we can conclude that 𝐸𝜖^𝑘2≤Γℎ²∀ℎ∈(0,1), (29) that is, the mean-square order of global error of the SSTM is 1.

4. Stability of SSTM

We are concerned with the stability of SSTM solution. Since (1) has an equilibrium solution 𝐱(𝑡)≡𝟎, we will discuss whether the SSTM solution 𝑋_𝑘 with a positive step-size can attain a similar stability when 𝑘 goes to infinity. First we give a sufficient condition for the exponential stability in mean-square sense of the equilibrium solution. The references [13, 19] give the condition as ^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗+^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗≤𝛼_𝑖^𝑛_𝑗=1𝑎_𝑗𝑖+𝛽_𝑖^𝑛_𝑗=1𝑏_𝑗𝑖,−2𝑐_𝑖+^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗+^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗+𝛼_𝑖^𝑛_𝑗=1𝑎_𝑗𝑖+𝛽_𝑖^𝑛_𝑗=1𝑏_𝑗𝑖+𝛾²_𝑖<0 (30) for every 𝑖(𝑖=1,…,𝑛).

Definition 7. A numerical method is said to be mean-square stable (MS-stable) if there exists an ℎ₀>0 such that any application of the method to problem (1) generates numerical approximations 𝑋^𝑘_𝑖, which satisfies lim_𝑘→∞𝐸𝑋^𝑘_𝑖²=0,𝑖=1,2,…,𝑛, (31) for all ℎ∈(0,ℎ₀).

Theorem 8. Assume (6), (30), and (1−𝜃)𝑐_𝑖ℎ≤1 are satisfied; then the SSTM (2) are mean-square stable if ℎ∈(0,ℎ₀), where ℎ₀=min_{1≤𝑖≤𝑛}1,ℎ_𝑖𝑓𝑜𝑟𝜃=1,min_{1≤𝑖≤𝑛}11,(1−𝜃)𝑐_𝑖,ℎ_𝑖𝑓𝑜𝑟𝜃∈[0,1). (32) Here ℎ_𝑖 is the smallest positive root of the cubic equation with respect to 𝑧 given by 𝒜_𝑖𝑧³+ℬ_𝑖𝑧²+𝒞_𝑖𝑧+𝒟_𝑖=0, (33) where the coefficients mean 𝒜_𝑖=𝜅²_𝑖2−(1−𝜃)𝑐_𝑖+^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗+^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗²,ℬ_𝑖=𝛾²_𝑖−(1−𝜃)𝑐_𝑖+^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗+^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗²+𝜅²_𝑖^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗+^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗−𝜅²_𝑖(1−𝜃)𝑐_𝑖,𝒞_𝑖=−(1−𝜃)𝑐_𝑖+^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗+^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗²+2𝛾²_𝑖−(1−𝜃)𝑐_𝑖+^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗+^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗−𝑐²_𝑖𝜃²+𝜅²_𝑖2,𝒟_𝑖=2−𝑐_𝑖+^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗+^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗+𝛾²_𝑖. (34) for 𝑖=1,2,…,𝑛.

Proof. Squaring on both sides of (4b) and (4a), we have 𝑋_𝑖^𝑘+12=𝑌^𝑘_𝑖²+𝜙_𝑖𝑌^𝑘_𝑖Δ𝑊^𝑘_𝑖²+𝜙^′_𝑖𝑌^𝑘_𝑖𝜙_𝑖𝑌^𝑘_𝑖²Δ𝑊^𝑘_𝑖²−ℎ2²+𝜙_𝑖𝑌^𝑘_𝑖Δ𝑊^𝑘_𝑖+𝑌^𝑘_𝑖×𝜙^′_𝑖𝑌^𝑘_𝑖𝜙_𝑖𝑌^𝑘_𝑖Δ𝑊^𝑘_𝑖²−ℎ+2𝑌^𝑘_𝑖𝜙_𝑖𝑌^𝑘_𝑖Δ𝑊^𝑘_𝑖,1+𝜃𝑐_𝑖ℎ²𝑌^𝑘_𝑖²=1−(1−𝜃)𝑐_𝑖ℎ𝑋^𝑘_𝑖²+ℎ^2𝑛_𝑗=1𝑎_𝑖𝑗𝑓_𝑗𝑋^𝑘_𝑗²+ℎ^2𝑛_𝑗=1𝑏_𝑖𝑗𝑔_𝑗𝑍^𝑘_𝑗²+^𝑛_𝑗=12ℎ1−(1−𝜃)𝑐_𝑖ℎ𝑋^𝑘_𝑖×𝑎_𝑖𝑗𝑓_𝑗𝑋^𝑘_𝑗+𝑏_𝑖𝑗𝑔_𝑗𝑍^𝑘_𝑗+2ℎ^2𝑛_𝑗=1𝑎_𝑖𝑗𝑓_𝑗𝑋^𝑘_𝑗^𝑛_𝑗=1𝑏_𝑖𝑗𝑔_𝑗𝑍^𝑘_𝑗. (35)(36) Taking expectations of both sides of (35), we can get 𝐸𝑋_𝑖^𝑘+12=𝐸𝑌^𝑘_𝑖²+ℎ𝐸𝜙_𝑖𝑌^𝑘_𝑖²+ℎ²2𝐸𝜙^′_𝑖𝑌^𝑘_𝑖𝜙_𝑖𝑌^𝑘_𝑖²≤1+𝛾²_𝑖ℎℎ+²2𝜅²_𝑖𝐸𝑌^𝑘_𝑖². (37) Together with (36), we have 𝐸𝑋_𝑖^𝑘+12≤1+𝛾²_𝑖ℎℎ+²2𝜅²_𝑖𝐸𝑌^𝑘_𝑖²≤1+𝛾²_𝑖ℎ+ℎ²/2𝜅²_𝑖1+𝜃𝑐_𝑖ℎ²×1−(1−𝜃)𝑐_𝑖ℎ²𝐸𝑋^𝑘_𝑖²+ℎ^2𝑛_𝑟=1𝑎_𝑖𝑟𝛼_𝑟^𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗𝐸𝑋^𝑘_𝑗²+ℎ^2𝑛_𝑟=1𝑏_𝑖𝑟𝛽_𝑟^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗𝐸𝑍^𝑘_𝑗²+^𝑛_𝑗=1ℎ1−(1−𝜃)𝑐_𝑖ℎ×𝑎_𝑖𝑗𝛼_𝑗𝐸𝑋^𝑘_𝑖²+𝐸𝑋^𝑘_𝑗²+𝑏_𝑖𝑗𝛽_𝑗𝐸𝑋^𝑘_𝑖²+𝐸𝑍^𝑘_𝑗²+2ℎ^2𝑛_𝑗=1𝑎_𝑖𝑗𝛼_𝑗×^𝑛_𝑗=1𝑏_𝑖𝑗𝛽_𝑗𝐸𝑋^𝑘_𝑗²+𝐸𝑍^𝑘_𝑗². (38) Thus, we attain 𝐸𝑋_𝑖^𝑘+12≤𝑃(ℎ)+^𝑛_𝑗=1𝑄_𝑗(ℎ)+^𝑛_𝑗=1𝑅_𝑗(ℎ)×max_{1≤𝑗≤𝑛}𝐸𝑋^𝑘_𝑖²,𝐸𝑋^𝑘_𝑗²,𝐸𝑍^𝑘_𝑗², (39) where 𝑃(ℎ)=1+𝛾²_𝑖ℎ+ℎ²𝜅²_𝑖/21+𝜃𝑐_𝑖ℎ²×1−(1−𝜃)𝑐_𝑖ℎ²+ℎ^𝑛_𝑗=11−(1−𝜃)𝑐_𝑖ℎ×𝑎_𝑖𝑗𝛼_𝑗+𝑏_𝑖𝑗𝛽_𝑗,𝑄_𝑗(ℎ)=1+𝛾²_𝑖ℎ+ℎ²𝜅²_𝑖/21+𝜃𝑐_𝑖ℎ²×ℎ²𝑎_𝑖𝑗𝛼_𝑗^𝑛_𝑗=1𝑎_𝑖𝑟𝛼_𝑟+𝑏_𝑖𝑟𝛽_𝑟+ℎ1−(1−𝜃)𝑐_𝑖ℎ𝑎_𝑖𝑗𝛼_𝑗,𝑅_𝑗(ℎ)=1+𝛾²_𝑖ℎ+ℎ²𝜅²_𝑖/21+𝜃𝑐_𝑖ℎ²×ℎ²𝑎_𝑖𝑗𝛼_𝑗+𝑏_𝑖𝑗𝛽_𝑗^𝑛_𝑟=1𝑏_𝑖𝑟𝛽_𝑟+ℎ1−(1−𝜃)𝑐_𝑖ℎ𝑏_𝑖𝑗𝛽_𝑗. (40)(41)(42) Note that the assumption of Theorem implies the nonnegativity of 1−(1−𝜃)𝑐_𝑖ℎ.
Obviously, when 𝑘→+∞, 𝐸(|𝑋^𝑘_𝑖|²)→0 if the inequality 𝑃(ℎ)+∑^𝑛_𝑗=1𝑄_𝑗(ℎ)+∑^𝑛_𝑗=1𝑅_𝑗(ℎ)<1 holds, which is equivalent to the inequality 𝒜_𝑖ℎ³+ℬ_𝑖ℎ²+𝒞_𝑖ℎ+𝒟_𝑖<0.
Furthermore, it is easy to prove 𝒜_𝑖>0 and 𝒟_𝑖<0 by virtue of (30). By Vieta's formulas, the product of three roots of (33) satisfies 𝑧₁𝑧₂𝑧₃=−(𝒟_𝑖/𝒜_𝑖)>0. This means that (33) has at least one positive root. Therefore, let ℎ_𝑖 denote the smallest positive root of the equation. Moreover, we note that at the origin the right-hand side polynomial of (33) is negative. This completes the proof.

5. Numerical Results

Now, we apply the introduced SSTM method to two test cases of SDHNNs in order to compare their performance with the split-step 𝜃-method in [13], which has strong convergence order 0.5.

The mean-square error 𝜀 of numerical approximations at time 𝑇 versus the step-size is depicted in log-log diagrams, where 𝜀∶=(1/2000)∑²⁰⁰⁰_𝑟=1‖𝑋^𝐾_𝜔𝑟−𝐱(𝑇)‖². Here 𝐱(𝑇) stands for the value of explicit solution of (1) at time 𝑇=3 and 𝑋^𝐾_𝜔𝑟 is its numerical approximation along the 𝑟th sample path {𝜔_𝑟∶𝑟=1,2,…,2000}. We compute the numerical solution using the split-step 𝜃-Milstein method (2) with step-size ℎ=2⁻¹², and we will call this the “exact solution.”

Example 9. Consider the following two-dimensional stochastic delay Hopfield neural networks of the form 𝑑𝑥₁(𝑡)𝑥₂(𝑡)𝑥=−𝐶¹(𝑡)𝑥₂(𝑡)d𝑓𝑥𝑡+𝐴¹(𝑡)𝑓𝑥₂(𝑡)d𝑡𝑔𝑥+𝐵¹𝑡−𝜏₁𝑔𝑥₂𝑡−𝜏₂d𝑡+Φd𝑊(𝑡) (43) on 𝑡⩾0 with the initial condition 𝑥₁(𝑡)=𝑡+1,𝑡∈[−1,0] and 𝑥₂(𝑡)=𝑡+1,𝑡∈[−2,0].

Case 1. Let 𝑓(𝑥)=sin(𝑥),𝑔(𝑥)=arctan(𝑥), ,𝑥𝐶=200020,𝐴=4−563𝐵=−6431,Φ=¹(𝑡)00−√5𝑥₂(𝑡). (44)

Case 2. Let 𝑓(𝑥)=𝑥,𝑔(𝑥)=tanh(𝑥), 𝐶=,100010,𝐴=0000𝐵=5−115,Φ=1.5𝑥₁(𝑡)00−1.5𝑥₂(𝑡). (45)

In Figure 1, SSTM is applied with 7 different step-sizes: ℎ_𝑚=2^𝑚−12 for 𝑚=1,2,…,7. Two pairs of time delays (𝜏₁,𝜏₂) are set to (1,2) and (1.13,2.31). The first pair has common factor ℎ_𝑚; however, the second pair is incommensurable by ℎ_𝑚. The computation errors 𝜀 versus step-sizes ℎ are plotted on a log-log scale and the reference lines of slope 1 are added. It illustrates that SSTM has raised the strong order of the split-step 𝜃-method at least to 1 for SDHNNs [13].

(a) Case 1 with 𝜏1=1,𝜏2=2

(b) Case 2 with 𝜏1=1,𝜏2=2

(c) Case 1 with 𝜏1=1.13,𝜏2=2.31

(d) Case 2 with 𝜏1=1.13,𝜏2=2.31

Figure 1 

Errors versus step-sizes for the SDHNNs (43).

Next, Table 1 shows a comparison of stability intervals between the SST and the SSTM for (43). Two sets of the interval in the Table are calculated through Theorem 8 in this paper and Theorem 5.1 in [13]. It is easy to see that the stability intervals of the two methods are similar.

Table 1 

Calculated intervals of step-size for mean-square stability of the numerical schemes.

We know that Theorem 5.1 in [13] and Theorem 8 in this paper only give sufficient conditions of mean-square stability. Therefore the stability intervals given by these theorems are only subsets of real ones. To confirm the situation, we calculated the sample moments of the approximate solution and plotted them along the time 𝑡. Here the sample moment 𝜂^𝑖 means (1/2000)∑²⁰⁰⁰_𝑟=1‖𝑋^𝑖_𝜔𝑟‖² for the numerical solution 𝑋^𝑖_𝜔𝑟 approximating 𝐱(𝑡_𝑖) along the 𝑟th sample path. Figures 2, 3, 4, 5, and 6 depict the results by SST and SSTM in the log-scaled vertical axis. All the figures can give a rough estimate of the stability interval in each case.

(a) SST

(b) SSTM

Figure 2 

Mean-square stability comparison of SST and SSTM for the SDHNNs (43), Case 1, 𝜃=0.

(a) SST

(b) SSTM

Figure 3 

Mean-square stability comparison of SST and SSTM for the SDHNNs (43), Case 2, 𝜃=0.

(a) SST

(b) SSTM

Figure 4 

Mean-square stability comparison of SST and SSTM for the SDHNNs (43), Case 1, 𝜃=0.2.

(a) SST

(b) SSTM

Figure 5 

Mean-square stability comparison of SST and SSTM for the SDHNNs (43), Case 2, 𝜃=0.2.

(a) Case 1

(b) Case 2

Figure 6 

Mean-square stability comparison of SST and SSTM for the SDHNNs (43).

6. Concluding Remarks

We introduce the split-step 𝜃-Milstein method (SSTM), which exhibits higher strong convergence rate than the split-step 𝜃-method (SST, see [13]) for a stochastic delay Hopfield neural networks, and the scheme proposed in this paper can deal with incommensurable time delays which were not considered in [13]. We give the proof of convergence results, which has generally been omitted in the previous works on the same subject. By comparing the stability intervals of step size for the SST and SSTM for a test example, we find they exhibit similar mean-square stability.

In this paper, we have found a delay-independent sufficient condition for mean-square stability of split-step 𝜃-Milstein method applied to nonlinear stochastic delay Hopfield neural networks. Further, Figure 6 suggests that the value of ℎ₀, the right end-point of the stability interval, given by Theorem 5.1 in [13] and Theorem 8 in this paper is much smaller than the true value when 𝜃 is close to unity. In this case, we need other techniques for stability analysis in this kind of stochastic delay differential system. To the best of our knowledge, the works in [20, 21] put forward good attempts. On the other hand, with respect to stochastic delay differential equations, some other types of stability have been successfully discussed for the Euler-type scheme, for example, mean-square exponential stability [12], delay-dependent stability [22], delay-dependent exponential stability [23], and almost sure exponential stability [24]. To Milstein-type scheme, in view of more sophisticated derivations, these issues would be challenging for future research.

Acknowledgment

The authors would like to thank Dr. E. Buckwar and Dr. A. Rathinasamy for their valuable suggestions. The authors also thank the editor and the anonymous referees, whose careful reading and constructive feedback have improved this work. The first author is partially supported by E-Institutes of Shanghai Municipal Education Commission (no. E03004), National Natural Science Foundation of China (no. 10901106), Natural Science Foundation of Shanghai (no. 09ZR1423200), and Innovation Program of Shanghai Municipal Education Commission (no. 09YZ150). The third author is partially supported by the Grants-in-Aid for Scientific Research (no. 24540154) supplied by Japan Society for the Promotion of Science.