Approximation of an Integral Markov Process Arising in the Approximation of Stochastic Differential Equation

Mohammad Rahman

Department of Mathematics & Statistics, University of North Florida, Jacksonville, FL-32224, USA.

DOI: 10.4236/apm.2022.121003                

Abstract

We provide the derivation of a new formula for the approximation of an integral Markov process arising in the approximation of stochastic differential equations. This formula extends an existing formula derived in [1]. We have shown numerically that the leading order approximation of the differential equation with noise by solving an associated averaged problem and estimating the difference between them and the result is illustrated through some examples.

Keywords

Variance, Markov Process, Parametric Noise, Differential Equation, Approximations

1. Introduction

Nonlinear ordinary and partial differential equations arise in various fields of sciences, particularly fluid mechanics, solid state physics, plasma physics, nonlinear optics, and mathematical biology. Powerful numerical methods and its implementation can be obtained in [2] - [7]. However, nonlinear differential equations (DE) with parametric noise play a significant role in a range of application areas, including engineering, physics, mechanics, epidemiology, and neuroscience. It is important to mention that noisy systems can be modeled in several ways: for example, Langevin’s equation describes a linear physical system to which white noise is added, and the linear theory for it has been extended to nonlinear stochastic differential equations with additive white noise [8]. Another approach is to derive models, such as Markov chains, for the systems state variables as being random, and then use the method of probability for analysis, see [9]. Third approach was derived by averaging nonlinear oscillatory systems. In this approach, parameters in the system are allowed to be random processes, and the method based on averaging and ergodic theory provides useful predictions from the model [10]. Indeed, solutions of differential or dynamical systems are functions, say of time. If they are in addition random, we must describe both randomness and time dependence simultaneously. Thus we refer to the system as either a random process, or a stochastic process. A complete understanding of DE theory with perturbed noise requires familiarity with advanced probability and stochastic processes (see [8] [11] ).

In this paper, our approach to modeling randomness in differential or dynamical system is through allowing parameters b in a system to be a random process. For example, in the case of differential equation we write $\stackrel{\dot{}}{x}\left(t\right)=f\left(t\mathrm{,}\omega \mathrm{,}x\left(t\right)\mathrm{,}b\left(y\left(t\right)\right)\right)$, indicating that the parameters can change with some random process $y\left(t\right)$. The resulting solutions $x\left(t\mathrm{,}b\right)$ will also be a random process. This approach will be based on assuming noise process faster time scale than the system time scale. This work has been dedicated to the question in particular when y is a discrete-space Markov process. Most recent contributions have aimed in general at relaxing the assumptions on y [1] [12]. This result is known as Functional Central Limit Theorem [13] [14] [15].

In Section 2 we reviewed and discussed limit of variance for discrete-space Markov Chain [1]. In Section 3, a new equivalent expression for a limit of a variance is given for circulant n-state markov chain. In Section 3, stochastic approximation and numerical simulation were discussed. It is observed that computer simulations of this type of stochastic ordinary differential equation with standard methods have some issues needed to be explored so that the reader will be benefitted while solving these type problems numerically.

2. Discrete-Space Markov Chains

We consider the case where $Y=\left\{y_1\mathrm{,}\cdots \mathrm{,}{y}_n\right\}$, and let

$\underset{\_}{y}=\left[\begin{array}{c}{y}_1\\ ⋮\\ y_n\end{array}\right],\phi \left(\underset{\_}{y}\right)=\left[\begin{array}{c}\phi \left(y_1\right)\\ ⋮\\ \phi \left(y_n\right)\end{array}\right],\text{and}\underset{\_}{e}=\left[\begin{array}{c}1\\ ⋮\\ 1\end{array}\right].$

The transition probability $P\left(\Delta t\mathrm{,}y\mathrm{,}d{y}^{\prime }\right)$ can be represented by an $n\times n$ matrix P such that $P_{i\mathrm{,}j}=P\left(\Delta t\mathrm{,}{y}_i\mathrm{,}{y}_j\right)$, denotes the probability that $y\left(\Delta t\right)=y_j$, if $y\left(0\right)=y_i$ (P depends on $\Delta t$ but this dependence is omitted from the notation for simplicity). The matrix P satisfies the following properties.

· P is nonnegative (its entries are probabilities).

· P is stochastic i.e., $P\underset{\_}{e}=\underset{\_}{e}$ (one of the $y_j$ must be the outcome of a transition from $y_i$ ).

Then ${\lambda }_1=1$, is an eigenvalue of P and $1\le \rho \left(P\right)\le {\Vert P\Vert }_{\infty }=1$, shows that all other eigenvalues have modulus at most 1. We shall assume that

· P is irreducible, i.e., any state $y_j$ can be eventually be reached in a finite number of steps with a nonzero probability (this is the case if $P>0$ ). This implies that ${\lambda }_1=1=\rho \left(P\right)$, has multiplicity one (e.g., see [16] ).

· P is aperiodic (or acyclic, e.g., not a permutation). This implies that $\left|{\lambda }_j\right|<1$ for $j=2,\cdots ,n$ (e.g., see [9]; the chain is then called regular).

The above conditions guarantee the existence of a unique vector $\underset{\_}{v}>0$ such that

${\underset{\_}{v}}^{\text{T}}\underset{\_}{e}=1,$ (2.1)

and ${\lim }_{N\to \infty }{P}^N=\underset{\_}{e}{\underset{\_}{v}}^{\text{T}}$. The vector $\underset{\_}{v}$ is the unique positive left eigenvector associated to ${\lambda }_1=1$, satisfying (2.1). It is natural to consider the (spectral) decomposition

$P=\underset{\_}{e}{\underset{\_}{v}}^{\text{T}}+S,$ (2.2)

where

${\underset{\_}{v}}^{\text{T}}S={\underset{\_}{0}}^{\text{T}},S\underset{\_}{e}=\underset{\_}{0},$ (2.3)

and $\rho \left(\tilde{P}\right)<1$. From Chapman-Kolmogorov equation and the homogeneity property (see, [1] ) we have

$P\left(2\Delta t,y_i,y_j\right)={\displaystyle \underset{k=1}{\overset{n}{\sum }}}P\left(\Delta t,y_i,y_k\right)P\left(\Delta t,y_k,y_j\right)={\displaystyle \underset{k=1}{\overset{n}{\sum }}}{P}_{i,k}{P}_{k,j}=P_{i,j}^2,$

and by induction

$P\left(N\Delta t,y_i,y_j\right)=P_{i,j}^N={\left(\underset{\_}{e}{\underset{\_}{v}}^{\text{T}}+S^N\right)}_{i,j}\approx {\left(\underset{\_}{e}{\underset{\_}{v}}^{\text{T}}\right)}_{i,j}=v_j,$

as $N\to \infty$, independently of i (i.e., $y_i$ ). Thus

$\underset{\_}{v}={\rho }_{\Delta t}\left(\underset{\_}{y}\right)\mathrm{,}$ (2.4)

defines the limit distribution. An explicit expression of the coefficients of $\underset{\_}{v}$ in terms of the coefficients of P can be found in [ [17], p. 21].

The relation of the expected value yields

$\begin{array}{c}{E}_{y\left(0\right)=y_j}\left[\phi \left(y\left(N\Delta t\right)\right)\right]={\displaystyle \underset{j=1}{\overset{n}{\sum }}}\phi \left(y_j\right)P\left(N\Delta t,y_j,y_i\right)\\ ={\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\left(P^N\right)}_{i,j}\phi \left(y_j\right)\\ ={\left(P^N\phi \left(\underset{\_}{y}\right)\right)}_j,\end{array}$ (2.5)

i.e.,

$E_{y\left(0\right)=\underset{\_}{y}}\left[\phi \left(y\left(N\Delta t\right)\right)\right]=\left[\begin{array}{c}{E}_{y\left(0\right)=y_1}\left[\phi \left(y\left(N\Delta t\right)\right)\right]\\ ⋮\\ E_{y\left(0\right)=y_n}\left[\phi \left(y\left(N\Delta t\right)\right)\right]\end{array}\right]=P^N\phi \left(\underset{\_}{y}\right).$

Then the zero average condition on $\phi$ becomes

$0={\displaystyle {\int }_Y}\phi \left(y\right)\rho \left(\text{d}y\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}}\phi \left(y_j\right)\rho \left(y_j\right)=\rho {\left(\underset{\_}{y}\right)}^{\text{T}}\phi \left(\underset{\_}{y}\right)={\underset{\_}{v}}^{\text{T}}\phi \left(\underset{\_}{y}\right).$ (2.6)

The relation (2.6) implies $P\phi \left(y\right)=\phi \left(y\right)$. Therefore

${\mathcal{R}}^{\Delta t,1}\phi \left(\underset{\_}{y}\right)=\Delta t{\displaystyle \underset{N=1}{\overset{\infty }{\sum }}}{E}_{y\left(0\right)=\underset{\_}{y}}\left[\phi \left(y\left(N\Delta t\right)\right)\right]$

$=\Delta t{\displaystyle \underset{N=1}{\overset{\infty }{\sum }}}{P}^N\phi \left(\underset{\_}{y}\right)$ (2.7)

$=\Delta t{\displaystyle \underset{N=1}{\overset{\infty }{\sum }}}{S}^N\phi \left(\underset{\_}{y}\right)$ (2.8)

$=\Delta t{\left(I-S\right)}^{-1}S\phi \left(\underset{\_}{y}\right),$ (2.9)

where $\frac{1}{\sqrt{\epsilon }}{\mathcal{R}}^{\Delta t\mathrm{,1}}\phi \left(\underset{\_}{y}\right)$ can be interpreted as the expected value of the random variable $\frac{1}{\sqrt{\epsilon }}{\mathcal{R}}^{\Delta t\mathrm{,1}}\phi \left(\underset{\_}{y}\right)$ obtained after one transition probability applied to the random variable y (for more detail, see [1] ). Note that (2.7) implies

${\mathcal{R}}^{\Delta t\mathrm{,1}}\phi \left(\underset{\_}{y}\right)-P{\mathcal{R}}^{\Delta t\mathrm{,1}}\phi \left(\underset{\_}{y}\right)=\Delta tP\phi \left(\underset{\_}{y}\right)\mathrm{.}$

Because of (2.3) we also have

${\underset{\_}{v}}^{\text{T}}{\mathcal{R}}^{\Delta t\mathrm{,1}}\phi \left(\underset{\_}{y}\right)=\mathrm{0,}$ (2.10)

With $V=\text{diag}\left(\underset{\_}{v}\right)$ we obtain

$\begin{array}{c}{\sigma }_{\Delta t}^2=2{\displaystyle \underset{j=1}{\overset{n}{\sum }}}\phi \left(y_j\right)\rho \left(y_j\right){\left({\mathcal{R}}^{\Delta t\mathrm{,1}}\phi \left(\underset{\_}{y}\right)\right)}_j+\Delta t{\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\phi }^2\left(y_j\right)\rho \left(y_j\right)\\ =2\Delta t\phi {\left(\underset{\_}{y}\right)}^{\text{T}}V{\left(I-S\right)}^{-1}S\phi \left(\underset{\_}{y}\right)+\Delta t\phi {\left(\underset{\_}{y}\right)}^{\text{T}}V\phi \left(\underset{\_}{y}\right)\\ =\Delta t\phi {\left(\underset{\_}{y}\right)}^{\text{T}}V{\left(I-S\right)}^{-1}\left(I+S\right)\phi \left(\underset{\_}{y}\right)\mathrm{,}\end{array}$ (2.11)

as $t\to \infty$. The expression (2.11) represents the limit of a variance and thus expected to be nonnegative.

Two-State Markov Chain

Let $n=2$ and

$P=P\left(\Delta t\right)=I+\Delta tQ,Q=\left[\begin{array}{cc}-a& a\\ b& -b\end{array}\right],$ (2.12)

with $0<a,b<1$, $0<\Delta t<\min \left(\frac{1}{a},\frac{1}{b}\right)\le \frac{2}{a+b}$. Then (2.2) holds with

$\begin{array}{l}\underset{\_}{v}=\frac{1}{a+b}\left[\begin{array}{c}b\\ a\end{array}\right]\\ \text{and}S=\frac{1-\Delta t\left(a+b\right)}{a+b}\left[\begin{array}{cc}a& -a\\ -b& b\end{array}\right].\end{array}$

As a result $V=\frac{1}{a+b}\left[\begin{array}{cc}b& \\ & a\end{array}\right]$ and

$\Delta tV{\left(I-S\right)}^{-1}\left(I+S\right)=\frac{1}{{\left(a+b\right)}^3}\left[\begin{array}{cc}b\left(2a+\Delta t\left(b^2-a^2\right)\right)& -2ab\left(1-\Delta t\left(a+b\right)\right)\\ -2ab\left(1-\Delta t\left(a+b\right)\right)& a\left(2b+\Delta t\left(a^2-b^2\right)\right)\end{array}\right]$ (2.13)

is (symmetric) positive definite for any choice $0<a,b<1$ and $0<\Delta t<\frac{2}{a+b}$. A simplified expression for (2.11) is obtained using (2.6), i.e.,

$\phi \left(\underset{\_}{y}\right)=\left[\begin{array}{c}1\\ -\frac{b}{a}\end{array}\right]\phi \left(y_1\right).$

We obtain

${\sigma }_{\epsilon }^2=\frac{b\left(2-\Delta t\left(a+b\right)\right)}{a\left(a+b\right)}{\left(\phi \left(y_1\right)\right)}^2$ (2.14)

$=\frac{ab\left(2-\Delta t\left(a+b\right)\right)}{{\left(a+b\right)}^3}{\left(\phi \left(y_1\right)-\phi \left(y_2\right)\right)}^2.$

3. Circulant n-State Markov Chain

If P is doubly stochastic (i.e., P and P^T are stochastic) and irreducible aperiodic then

$\underset{\_}{v}=\frac{1}{n}\underset{\_}{e}\mathrm{,}$

i.e., the limit distribution is uniform and

$D=\frac{1}{n}I.$ (3.16)

Stochastic Toeplitz, hence circulant, matrices constitute an example of doubly stochastic matrices. Let

$P=\left[\begin{array}{cccc}{a}_1& a_2& \ddots & a_n\\ a_n& \ddots & \ddots & \ddots \\ \ddots & \ddots & \ddots & a_2\\ a_2& \ddots & a_n& a_1\end{array}\right],$

with $a_i\ge 0$ and ${\displaystyle {\sum }_{j=1}^n}{a}_j=1$. A sufficient condition for P to be irreducible and aperiodic is for two consecutive $a_j$ to be nonzero (i.e., positive, see e.g. [ [18], p. 5]). The symbol of P is the polynomial

$p\left(z\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}}{a}_j{z}^{j-1}.$

Since

$P\left[\begin{array}{c}1\\ z\\ ⋮\\ z^{n-1}\end{array}\right]=p\left(z\right)\left[\begin{array}{c}1\\ z\\ ⋮\\ z^{n-1}\end{array}\right]\mathrm{,}$

for $z^n=1$, P admits the spectral decomposition

$P=V\Lambda V^{-1}=V\Lambda V^H,$

with

$V=\frac{1}{\sqrt{n}}\left[\begin{array}{cccc}1& 1& \cdots & 1\\ 1& z_1& \cdots & z_{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ 1& z_1^{n-1}& \cdots & z_{n-1}^{n-1}\end{array}\right]=\left[\begin{array}{cc}\frac{1}{\sqrt{n}}\underset{\_}{e}& \tilde{V}\end{array}\right],$ (3.17)

$\Lambda =\left[\begin{array}{cccc}1& & & \\ & p\left(z_1\right)& & \\ & & \ddots & \\ & & & p\left(z_{n-1}\right)\end{array}\right]=\left[\begin{array}{cc}1& \\ & \tilde{\Lambda }\end{array}\right],$ (3.18)

where $z_j={\text{e}}^{\frac{2i\pi j}{n}}=z_1^j$. Multiplication of a vector by the matrix $V^H$ performs a (normalized) discrete Fast Fourier Transform (FFT), while multiplication by V results in the inverse (normalized) discrete FFT.

We write $P=\frac{1}{n}\underset{\_}{e}{\underset{\_}{e}}^{\text{T}}+\tilde{V}\tilde{\Lambda }{\tilde{V}}^H=\frac{1}{n}\underset{\_}{e}{\underset{\_}{e}}^{\text{T}}+\tilde{P}$, with $\tilde{P}=\tilde{V}\tilde{\Lambda }{\tilde{V}}^H$. Since $\left|p\left(z_k\right)\right|<1$ for $k=1,\cdots ,n-1$ the matrix $I-\tilde{\Lambda }$ is nonsingular. The matrix $\tilde{V}{\tilde{V}}^H$ represents the (orthogonal) projection onto $\text{Span}{\left\{\underset{\_}{v}\right\}}^{\perp }=\text{Span}{\left\{\underset{\_}{e}\right\}}^{\perp }$. Hence ${\underset{\_}{v}}^{\text{T}}\phi \left(\underset{\_}{y}\right)=0$ implies $\phi \left(\underset{\_}{y}\right)=\tilde{V}{\tilde{V}}^H\phi \left(\underset{\_}{y}\right)$. Then

$\begin{array}{c}\left(I+\tilde{P}\right)\phi \left(\underset{\_}{y}\right)=\left(I+\tilde{V}\tilde{\Lambda }{\tilde{V}}^H\right)\tilde{V}{\tilde{V}}^H\phi \left(\underset{\_}{y}\right)\\ =\tilde{V}\left(I+\tilde{\Lambda }\right){\tilde{V}}^H\phi \left(\underset{\_}{y}\right)\\ =\tilde{V}\left(I-\tilde{\Lambda }\right){\left(I-\tilde{\Lambda }\right)}^{-1}\left(I+\tilde{\Lambda }\right){\tilde{V}}^H\phi \left(\underset{\_}{y}\right)\\ =\left(I-\tilde{P}\right)\tilde{V}{\left(I-\tilde{\Lambda }\right)}^{-1}\left(I+\tilde{\Lambda }\right){\tilde{V}}^H\phi \left(\underset{\_}{y}\right)\mathrm{.}\end{array}$ (3.19)

Now

$\begin{array}{c}{\left({\tilde{V}}^H\phi \left(\underset{\_}{y}\right)\right)}_k=\frac{1}{\sqrt{n}}\left[\begin{array}{cccc}1& \overline{z_k}& \cdots & {\overline{z_k}}^{n-1}\end{array}\right]\left[\begin{array}{c}\phi \left(y_1\right)\\ \phi \left(y_2\right)\\ ⋮\\ \phi \left(y_n\right)\end{array}\right]\\ =\frac{1}{\sqrt{n}}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}\phi \left(y_j\right){\overline{z_k}}^{j-1}\\ =\overline{\Phi \left(z_k\right)},\end{array}$

for $k=1,\cdots ,n-1$ ( $\phi \left(y\right)$ is real), with

$\Phi \left(z\right)=\frac{1}{\sqrt{n}}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}\phi \left(y_j\right)z^{j-1}.$ (3.20)

From (3.19) and (3.16) we obtain

${\sigma }_{\Delta t}^2=\frac{\Delta t}{n}{\left({\tilde{V}}^H\phi \left(\underset{\_}{y}\right)\right)}^H{\left(I-\tilde{\Lambda }\right)}^{-1}\left(I+\tilde{\Lambda }\right)\left({\tilde{V}}^H\phi \left(\underset{\_}{y}\right)\right)$

$=\frac{\Delta t}{n}{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}}\frac{1+p\left(z_k\right)}{1-p\left(z_k\right)}{\left|\Phi \left(z_k\right)\right|}^2$ (3.21)

$=\frac{\Delta t}{n}{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}}\Re \left(\frac{1+p\left(z_k\right)}{1-p\left(z_k\right)}\right){\left|\Phi \left(z_k\right)\right|}^2$ (3.22)

$=\frac{\Delta t}{n}{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}}\frac{1-{\left|p\left(z_k\right)\right|}^2}{{\left|1-p\left(z_k\right)\right|}^2}{\left|\Phi \left(z_k\right)\right|}^2,$ (3.23)

since ${\sigma }_{\Delta t}^2$ is real and $\overline{p\left(z_k\right)}=p\left(\overline{z_k}\right)$, for $k=1,\cdots ,n-1$. The condition $\left|p\left(z_k\right)\right|<1$ for $k=1,\cdots ,n-1$ clearly shows that (3.23) is nonnegative.

3.1. Example: Uniform Distribution

If $a_j=\frac{1}{n}$, $j=1,\cdots ,n$, (i.e., $P=\frac{1}{n}\underset{\_}{e}{\underset{\_}{e}}^{\text{T}}$ ) we obtain

$p\left(z_k\right)=\frac{1}{n}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}{z}_k^{j-1}=\frac{1}{n}\frac{1-z_k^n}{1-z_k}=0$

for $k=1,\cdots ,n-1$. Parseval’s identity yields

$\begin{array}{c}{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}}{\left|\Phi \left(z_k\right)\right|}^2={\displaystyle \underset{k=0}{\overset{n-1}{\sum }}}{\left|\Phi \left(z_k\right)\right|}^2={\displaystyle \underset{k=0}{\overset{n-1}{\sum }}}{\left|{\left(V^H\phi \left(\underset{\_}{y}\right)\right)}_k\right|}^2\\ =\phi {\left(\underset{\_}{y}\right)}^{\text{T}}V{V}^H\phi \left(\underset{\_}{y}\right)=\phi {\left(\underset{\_}{y}\right)}^{\text{T}}\phi \left(\underset{\_}{y}\right)={\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\left|\phi \left(y_j\right)\right|}^2.\end{array}$

Then (3.23) reduces to

${\sigma }_{\Delta t}^2=\frac{\Delta t}{n}{\displaystyle \underset{k=1}{\overset{n-1}{\sum }}}{\left|\Phi \left(z_k\right)\right|}^2=\frac{\Delta t}{n}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\left|\phi \left(y_j\right)\right|}^2\approx \Delta t{\displaystyle {\int }_Y}{\phi }^2\left(y\right){\rho }_{\Delta t}\left(\text{d}y\right),\text{as}n\to \infty .$

3.2. Example: Big World Transition Probability

Assume now that the probability $y_i\to y_j$, $j\ne i$, is independent of j but distinct from the probability $y_i\to y_i$, i.e.,

$a_1=1-a\Delta t,a_2=\cdots =a_n=\frac{a\Delta t}{n-1},$

with $0<a<1/\Delta t$. Then

$p\left(z\right)=1-a\Delta t+\frac{a\Delta t}{n-1}\left(z+\cdots +z^{n-1}\right)=1-a\Delta t+\frac{a\Delta t}{n-1}\frac{z-z^n}{1-z}$

for $z\ne 1$. In particular

$p\left(z_k\right)=1-a\Delta t+\frac{a\Delta t}{n-1}\frac{z_k-1}{1-z_k}=1-\frac{na\Delta t}{n-1},$

for $k=1,\cdots ,n-1$. Therefore, by Parseval again,

$\begin{array}{c}{\sigma }_{\Delta t}^2=\frac{2-\frac{na\Delta t}{n-1}}{\frac{na}{n-1}}\frac{1}{n}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\left|\phi \left(y_j\right)\right|}^2\\ \approx \frac{2-\Delta ta}{a}{\displaystyle {\int }_Y}{\phi }^2\left(y\right){\rho }_{\Delta t}\left(\text{d}y\right)\text{as}n\to \infty \\ \approx \frac{2}{a}{\displaystyle {\int }_Y}{\phi }^2\left(y\right){\rho }_{\Delta t}\left(\text{d}y\right)\text{as}\Delta t\to 0.\end{array}$ (3.24)

3.3. Example: Small World Transition Probability

A case of particular interest in the study of the transmission of a signal around a cyclic biological chain corresponds to

$a_1=1-a\Delta t,a_2=a\Delta t,a_3=\cdots =a_n=0,$

with $0<a<1/\Delta t$. The quantity $a_2$ represents the transition probability of a state $y_j$ to the next state $y_{j+1}$. Then $p\left(z\right)=1-a\Delta t+az\Delta t$, and

$\frac{1-{\left|p\left({\text{e}}^{iy}\right)\right|}^2}{{\left|1-p\left({\text{e}}^{iy}\right)\right|}^2}=\frac{1-{\left(1-a\Delta t+a\Delta t\right)}^2-{\left(a\Delta t\right)}^2}{{\left(a\Delta t-a\Delta t\right)}^2+{\left(a\Delta t\right)}^2}=\frac{1-a\Delta t}{a\Delta t},$

for $0<y<2\pi$. Therefore

${\sigma }_{\Delta t}^2=\frac{1-a\Delta t}{a}\frac{1}{n}{\displaystyle \underset{j=1}{\overset{n}{\sum }}}{\left|\phi \left(y_j\right)\right|}^2$ (3.25)

$\approx \frac{1-a\Delta t}{a}{\displaystyle {\int }_Y}{\phi }^2\left(y\right){\rho }_{\Delta t}\left(\text{d}y\right)\text{as}n\to \infty$

$\approx \frac{1}{a}{\displaystyle {\int }_Y}{\phi }^2\left(y\right){\rho }_{\Delta t}\left(\text{d}y\right),\text{as}\Delta t\to 0.$ (3.26)

For $n=2$, (3.25) reduces to $\frac{1-\Delta ta}{a}\frac{\phi {\left(y_1\right)}^2+\phi {\left(y_2\right)}^2}{2}$. Note that (3.26) is half

of (3.24), as could be expected from a less dispersive signal.

In the following section our approach to modeling randomness in dynamical systems is through allowing parameters in a system to be a random process. This approach will determine when solution of this type of stochastic problem do or do not persist when the system is perturbed.

4. Mathematical Derivation for Numerical Approximation of Stochastic Differential Equations (SDE)

For the simplicity, consider

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=f\left(t,\omega ,x\left(t\right),y\right),\\ y=y\left(t/\epsilon \right),\\ x\left(0\right)=x_0,\end{array}$ (27)

The average system is defined by the differential equation

$\{\begin{array}{l}\stackrel{\dot{}}{\bar{x}}\left(t\right)=f\left(t\mathrm{,}\omega \mathrm{,}\bar{x}\left(t\right)\right)\mathrm{,}\\ \bar{x}\left(0\right)=x_0\mathrm{,}\end{array}$ (28)

We want to compute the deviation of the perturbed system to be average one. We consider $\tilde{x}=x-\bar{x}$, then

$\begin{array}{c}\stackrel{\dot{}}{\tilde{x}}=\stackrel{\dot{}}{x}-\stackrel{\dot{}}{\bar{x}}\\ =f\left(t\mathrm{,}\omega \mathrm{,}\bar{x}+\tilde{x}\mathrm{,}y\right)-\bar{f}\left(t\mathrm{,}\omega \mathrm{,}\bar{x}\right)\\ =f\left(t\mathrm{,}\omega \mathrm{,}\bar{x}\mathrm{,}y\right)+f\left(t\mathrm{,}\omega \mathrm{,}\bar{x}\mathrm{,}y\right)\tilde{x}+\cdots -\bar{f}\left(t\mathrm{,}\omega \mathrm{,}\bar{x}\right)\\ \simeq {\bar{f}}_x\left(t\mathrm{,}\omega \mathrm{,}\bar{x}\right)\tilde{x}+\left(f\left(t\mathrm{,}\omega \mathrm{,}\bar{x}\mathrm{,}y\right)-\bar{f}\left(t\mathrm{,}\omega \mathrm{,}\bar{x}\right)\right)\mathrm{.}\end{array}$

Note that $\tilde{x}\left(0\right)=0$, $\Rightarrow$

$\tilde{x}={\displaystyle {\int }_0^t{\bar{f}}_x\left(s\mathrm{,}\omega \mathrm{,}\bar{x}\left(s\right)\right)\tilde{x}\text{d}s}+{\displaystyle {\int }_0^t\left[f\left(s\mathrm{,}\omega \mathrm{,}\bar{x}\left(s\right)\mathrm{,}y\left(s/\epsilon \right)\right)-\bar{f}\left(s\mathrm{,}\omega \mathrm{,}\bar{x}\left(s\right)\right)\right]\text{d}s}\mathrm{.}$

It follows from limit theorem of stochastic processes (see, [2] ),

$\frac{1}{\sqrt{\epsilon }}{\displaystyle {\int }_0^t}\left[f\left(s\mathrm{,}\omega \mathrm{,}\bar{x}\left(s\right)\mathrm{,}y\left(s/\epsilon \right)\right)-\bar{f}\left(s\mathrm{,}\omega \mathrm{,}\bar{x}\left(s\right)\right)\right]\text{d}s\simeq N\left(\mathrm{0,}{\sigma }^2\left(t\right)\right)\mathrm{.}$

Detailed derivation of ${\sigma }^2\left(t\right)$ is shown in Section 2. The stochastic processes

$\frac{x\left(t\right)-\bar{x}\left(t\right)}{\sqrt{\epsilon }}\approx \tilde{x}\left(t\right)\mathrm{,}$

converge in expected sense to the solution $\tilde{x}\left(t\right)$ that is the solution to the integral equation

$\tilde{x}\left(t\right)={\displaystyle {\int }_0^t}{\bar{f}}_x\left(s\mathrm{,}\omega \mathrm{,}\bar{x}\left(s\right)\right)\tilde{x}\left(s\right)\text{d}s+\sigma \left(t\right)\mathrm{.}$

The distribution of $x\left(t\right)$ is close to the distribution of the stochastic processes $\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)$ in the sense that

$E\left(x\left(t\right)\right)=E\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\mathrm{.}$

Thus $x\left(t\right)\approx \bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)+o\left(\sqrt{\epsilon }\right)$, as $\epsilon \to 0$. Thus the nature of this convergence and the sense in which the error in the formula of the expansion x is small are in the sense of distributions over some time interval.

4.1. Forward Euler Scheme for SDE

The Euler method to approximate the analytic solution of the IVP

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=f\left(t\mathrm{,}x\mathrm{,}y\left(t/\epsilon \right)\right)\mathrm{,}\\ x\left(0\right)=x_0\mathrm{.}\end{array}$ (4.29)

We can derive an entire family of discrete numerical methods (including the Euler method) by truncating the Taylor series and utilizing Taylor’s Theorem. First, we rewrite the Taylor series expansion in differential form,

$\begin{array}{c}x\left(t+h\right)=x\left(t\right)+\frac{h}{\mathrm{1\; <mo>\; !\; </mo>}}{x}^{\prime }\left(t\right)+\frac{h^2}{\mathrm{2\; <mo>\; !\; </mo>}}{y}^{″}\left(t\right)+\cdots +\frac{h^n}{n\mathrm{<mo>\; !\; </mo>}}{y}^{″}\left(t\right)+\cdots \\ =x\left(t\right)+hf\left(t\mathrm{,}x\mathrm{,}y\left(t/\epsilon \right)\right)+\frac{h^2}{2}{f}^{\prime }\left(t\mathrm{,}x\mathrm{,}y\left(t/\epsilon \right)\right)f\left(t\mathrm{,}x\mathrm{,}y\left(t/\epsilon \right)\right)+\cdots \\ +\frac{h^3}{6}\left(f^{″}\left(f\mathrm{,}f\right)+f^{\prime }\left(f^{\prime }\left(f\right)\right)\right)+\cdots .\end{array}$

The Euler method is found by truncating Taylor series at the first derivative, giving

$x\left(t+h\right)=x\left(t\right)+hf\left(t\mathrm{,}x\mathrm{,}y\left(t/\epsilon \right)\right)+o\left(h^2\right)\mathrm{,}$ (4.30)

where $o\left(h^2\right)$ is the error term, or Taylors remainder term, which is of order $h^2$. So if we define $x_n=x\left(t+h\right)$, and $x_{n+1}=x\left(t+h\right)$, we obtain

$x_{n+1}=x_n+hf\left(n\mathrm{,}{x}_n\mathrm{,}y\left(n/\epsilon \right)\right)\mathrm{,}$

which has local error term (error at each step) of $o\left(h^2\right)$, giving a global error $o\left(h\right)$.

4.2. Discrete Approximation for SDE

Numerical schemes for solving SDE can be classified into either explicit or implicit methods. Explicit methods compute approximations that are dependent on previous approximations only, whereas the implicit methods compute approximations that are dependent on previous and current approximations. The Euler method presented earlier is also known as the explicit Euler method. The explicit Euler method is defined as

$x_{i+1}=x_i+hf\left(t_i,x_i,y_{i/\epsilon }\right),$ (4.31)

where $x_i=x\left(t_i\right)$ and $x_{i+1}=x\left(t_{i+h}\right)$. Here h is the step size. This method has a global error $\approx o\left(h\right)$. The implicit method is defined as

$x_{i+1}=x_i+hf\left(t_i,x_{i+1},y_{\left(i+1\right)/\epsilon }\right),$ (4.32)

where $x_i=x\left(t_i\right)$, and $x_{i+1}=x\left(t_i+h\right)$. This method has a global error $\approx o\left(h\right)$.

4.3. Example 1

We write the algorithm for the following equation, which we will be used as a test equation since it has an analytic solution.

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=y\left(t/\epsilon \right)\mathrm{,}\\ x\left(0\right)=0,\end{array}$ (4.33)

over $\left[\mathrm{0,1}\right]$ and then

$x\left(t+h\right)=x\left(t\right)+hy\left(t/\epsilon \right)+o\left(h\right)\mathrm{.}$ Now,

1) choose a step $h=\left(1-0\right)/N$. Set $x_n=0+nh$, $n=0:N$.

2) Generate approximation $x_n$ to from the following recursion:

$\begin{array}{c}{x}_n=0+h{\displaystyle {\sum }_{k=0}^{n-1}y\left(kh/\epsilon \right)}\\ =h{\displaystyle {\sum }_{k=0}^{n-1}y\left(kh/\epsilon \right)}\\ =hn\bar{y}+h{\displaystyle {\sum }_{l=0}^{h\left(n-1\right)/\epsilon }\left(y\left(l\right)-\bar{y}\right)}\\ =hn\bar{y}+H\sqrt{\epsilon }\frac{1}{1/\sqrt{\epsilon }}{\displaystyle {\sum }_{l=0}^{h\left(n-1\right)/\epsilon }\left(y\left(l\right)-\bar{y}\right)},\end{array}$

for $n=0:N-1$.

As $\epsilon \to 0$,

$x_n\to hn\bar{y}+H\sqrt{\epsilon }W\left(hn\right)+Error\mathrm{,}$

where $h=H\epsilon$ in our numerical simulation and W is realized by the normally distributed random number whose expectation and variance are 0 and 1, respectively. A computer simulation of the distribution of solution of $x\left(\mathrm{1,}\epsilon \right)$ using Euler without max step size, ode15s without max step, and ode15s using max step size are given in Figure 1.

Figure 1. Horizontal bar plot: x vs. bins. This figure shows the distribution of $x\left(\mathrm{1,}\epsilon \right)$ represented by the vertical axis of 500 sample paths versus bins in the horizontal axis at the final time t = 1 of the system (4.33). In this case, $x\left(0\right)=1$, =0.67, $0\le t\le 1$, and $H=0.1$, $\epsilon =0.01$. The left figure shows the distribution of $x\left(\mathrm{1,}\epsilon \right)$ using Euler without max step size. The middle figure shows the distribution of $x\left(\mathrm{1,}\epsilon \right)$ using ode15s solver without max step size. The right figure shows the distribution of $x\left(\mathrm{1,}\epsilon \right)$ using ode15s using max step size.

4.4. Example 2

We consider a stochastic process $x_{\epsilon }\left(t\right)$ is defined by the differential equation

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=-tx\left(y\left(\frac{t}{\epsilon }\right)-0.5\right)\equiv f\left(t\mathrm{,}\epsilon \mathrm{,}x\mathrm{,}y\right)\mathrm{,}\\ x\left(0\right)=1,\end{array}$ (4.34)

with

$Y=\left\{\mathrm{0,1}\right\}\mathrm{,}$

$P=\left(\begin{array}{cc}0.8& 0.2\\ 0.1& 0.9\end{array}\right),$

$Q=\left(\begin{array}{cc}-0.2& 0.2\\ 0.1& -0.1\end{array}\right),$

and limiting distribution

$\rho =\left\{1/3\mathrm{,}2/3\right\}$

We want to solve (4.34) numerically over time interval $\left[\mathrm{0,1}\right]$ using

1) Forward Euler with smaller time scale.

2) evaluate $\bar{x}+\sqrt{\epsilon }\tilde{x}$, where

$\stackrel{\dot{}}{\bar{x}}=-tx\left(\frac{1}{3}{y}_1+\frac{2}{3}{y}_2-0.5\right)\equiv \bar{f}\left(t\mathrm{,}\epsilon \mathrm{,}x\right)\mathrm{,}$ (4.35)

and

$\text{d}\tilde{x}={\bar{f}}_x\left(t\mathrm{,}\epsilon \mathrm{,}x\right)\tilde{x}\text{d}t+\sigma \left(t\right)\text{d}W\mathrm{.}$ (4.36)

3) compare 2. with 1., where $\bar{x}\left(t\right)$ and $\tilde{x}\left(t\right)$ are the solutions of (4.35) (using large time scale) and (4.36) respectively. And $\sigma \left(t\right)$ is computed using (2.15) along with $\rho \mathrm{,}f\mathrm{,}\bar{f}$, and $\text{d}W=\sqrt{\text{d}t}N\left(\mathrm{0,1}\right)$. Comparison of the averaged system to the perturbed system as well as error convergence are illustrated in Figure 2 and Figure 3.

4.5. Example 3

We consider a stochastic process $x_{\epsilon }\left(t\right)$ is defined by the differential equation

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=1/2+\cos \left(y\left(\frac{t}{\epsilon }\right)\right)\cos \left(x\right)\equiv f\left(t,\epsilon ,x,y\right),\\ x\left(0\right)=1,\end{array}$ (4.37)

with

$Y=\left\{\mathrm{0,1}\right\}\mathrm{,}$

Figure 2. Left figure: $x$ and $\bar{x}$ vs. t. Comparison of the averaged system to the perturbed system. Red curve shows the solution of the averaged system; other curves are solutions of the perturbed system with different $\epsilon$. Right figure: horizontal bar plot, where the vertical axis represents the distribution of solutions for 500 sample paths of the system (4.34) at the final time t = 1. In this case, $\epsilon =0.01$, $0\le t\le 1$, $y$ is generated with two state Markov process with $\bar{y}$ = 0.67.

Figure 3. Error plots: $E\left|x\left(t\right)-\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\right|$ represents vertical axis vs $\epsilon$, which represents horizontal axis. Dashed line is the line with the reference slope 1/2. Graphs are drawn on log-log scale.

$P=\left(\begin{array}{cc}0.8& 0.2\\ 0.1& 0.9\end{array}\right),$

$Q=\left(\begin{array}{cc}-0.2& 0.2\\ 0.1& -0.1\end{array}\right),$

and limiting distribution

$\rho =\left\{1/3\mathrm{,}2/3\right\}$

We want to solve (4.37) numerically over time interval $\left[\mathrm{0,1}\right]$ using

1) Forward Euler with smaller time scale.

2) evaluate $\bar{x}+\sqrt{\epsilon }\tilde{x}$ (using smaller time scale, h), where

$\stackrel{\dot{}}{\bar{x}}=1/2+\left(\frac{1}{3}\cos 0+\frac{2}{3}\cos 1\right)\cos \left(\bar{x}\right)\equiv \bar{f}\left(t,\epsilon ,x\right),$ (4.38)

and

$\text{d}\tilde{x}={\bar{f}}_x\left(t,\epsilon ,x\right)\tilde{x}\text{d}t+\sigma \left(t\right)\text{d}W.$ (4.39)

3) compare 2. with 1., where $\bar{x}\left(x\right)$ and $\tilde{x}\left(t\right)$ is the solution of (4.38) (using large time scale), and (4.39) respectively. And $\sigma \left(t\right)$ is computed using (2.15) along with $\rho \mathrm{,}f\mathrm{,}\bar{f}$ and $\text{d}W=\sqrt{\text{d}t}N\left(\mathrm{0,1}\right)$. Comparison of the averaged system to the perturbed system as well as error convergence are illustrated in Figure 4 and Figure 5.

4.6. Example 4

We consider a stochastic process $x_{\epsilon }\left(t\right)$ is defined by the differential equation

Figure 4. Left figure: $x$ and $\bar{x}$ vs. t. Comparison of the averaged system to the perturbed system. Red curve shows the solution of the averaged system; other curves are solutions of the perturbed system with different $\epsilon$. Right figure: horizontal bar plot, where the vertical axis represents the distribution of solutions for 500 sample paths of the system (4.37) at the final time t = 1. Here $\epsilon =0.01$, $0\le t\le 1$, $y$ is generated with two sate Markov process with $\bar{y}$ = 0.67.

Figure 5. Error plots: $E\left|x\left(t\right)-\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\right|$ represents vertical axis vs $\epsilon$, which represents horizontal axis. Dashed line is the line with slope 1/2. Graphs are drawn on log-log scale.

$\{\begin{array}{l}\stackrel{\dot{}}{x}\left(t\right)=1+y\left(\frac{t}{\epsilon }\right)-\left(3+y\left(\frac{t}{\epsilon }\right)\right)\cos \left(x\right)\equiv f\left(t,\epsilon ,x,y\right)\\ x\left(0\right)=1,\end{array}$ (4.40)

with

$Y=\left\{0,1\right\},$

$P=\left(\begin{array}{cc}0.8& 0.2\\ 0.1& 0.9\end{array}\right),$

$Q=\left(\begin{array}{cc}-0.2& 0.2\\ 0.1& -0.1\end{array}\right),$

and limiting distribution

$\rho =\left\{1/3,2/3\right\}.$

We want to solve (4.40) numerically over time interval $\left[\mathrm{0,1}\right]$ using

1) Forward Euler with smaller time scale.

2) evaluate $\bar{x}+\sqrt{\epsilon }\tilde{x}$ (using smaller time scale, h), where

$\stackrel{\dot{}}{\bar{x}}=1+y-\left(3+y\right)\cos \left(x\right)\equiv \bar{f}\left(t,\epsilon ,x\right).$ (4.41)

and

$\text{d}\tilde{x}={\bar{f}}_x\left(t,\epsilon ,x\right)\cos \left(\bar{x}\right)\tilde{x}\text{d}t+\sigma \left(t\right)\text{d}W.$ (4.42)

3) compare 2. with 1. where $\bar{x}\left(t\right)$ and $\tilde{x}\left(t\right)$ is the solution of (4.41) (using large time scale) and (4.42) respectively. And $\sigma \left(t\right)$ is computed using (2.15) along with $\rho \mathrm{,}f\mathrm{,}\bar{f}$, and $\text{d}W=\sqrt{\text{d}t}N\left(\mathrm{0,1}\right)$. Comparison of the averaged system to the perturbed system as well as error convergence are illustrated in Figure 6 and Figure 7.

4.7. Strong Convergence

In the examples above, directly simulated solution $x_{\epsilon }\left(t\right)$ with smaller time scale matches more closely to the solution $\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)$, (where $\bar{x}\left(t\right)$ solved using larger time scale and $\tilde{x}\left(t\right)$ is the solution of $\text{d}\tilde{x}={\bar{f}}_x\left(t,y,\bar{x}\right)\tilde{x}\text{d}t+\sigma \left(t\right)\text{d}w$, using Euler with large time scale H) as $\epsilon$ is decreased, the convergence seems to take place. Using $E\left|x_{\epsilon }\left(t\right)-\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\right|$, where E denotes the expected value, leads the concept of strong convergence. A method is said to have strong order of convergence equal to m if there exists a constants K such that

$E\left|x_{\epsilon }\left(t\right)-\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\right|\le K{\left(\epsilon \right)}^m\mathrm{,}$ (4.43)

and $\epsilon$ is sufficiently small. It can be shown that perturbation method has strong order of convergence $m=1$. In our numerical tests, we will focus on the error at the end point $t=tfinal$, so let

$Erro{r}^{strong}=E\left|x_{\epsilon }\left(t\right)-\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\right|.$ (4.44)

If the bound in (4.43) holds with $m=1$ at any fixed point in $\left[\mathrm{0,}tfinal\right]$, then it certainly holds at the end point, so we have

Figure 6. Left figure: $x$ and $\bar{x}$ vs. t. Comparison of the averaged system to the perturbed system. Red curve shows the solution of the averaged system; other curves are solutions of the perturbed system with different $\epsilon$. Right figure: horizontal bar plot, where the vertical axis represents the distribution of solutions for 500 sample paths of the system (4.40) at the final time t = 1. In this case, $\epsilon =0.01$, $0\le t\le 1$, $y$ is generated with two state Markov process with $\bar{y}$ = 0.67.

Figure 7. Error plots: The vertical axis represents $E\left|x\left(t\right)-\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\right|$ vs $\epsilon$, which represents the horizontal axis. Dashed line is the line with the reference slope 1/2. Graphs are drawn on log-log scale.

$Erro{r}^{strong}\le K{\epsilon }^1\mathrm{,}$ (4.45)

for sufficiently small $\epsilon$. It is shown that perturbation method has strong order of convergence $m=1$. While experimenting the error $Erro{r}^{strong}$, we implicitly assumed that number of other sources of error are negligible, including error arising from approximating an expected value by sample mean, inherent error in the random generator, and floating point roundoff errors. For a typical computation the sampling error is likely to be the most significant of these three. In preparing the programs for these simulations we found that some experimentation is required to make the number of samples sufficiently large and the time step is sufficiently small for the predicted order of convergence to be observable. The sampling error decays like $1/\sqrt{n}$, where n is the number of sample paths used. A study in ( [19] ) indicates that as step size decreases, the lack of independence in the samples from a random generator typically degrades the computation before rounding errors becomes significant.

Although the definition of strong convergence involves an expected value, it has implications for individual simulations. The Markov inequality says that if a random variable $\xi$ has a finite expected value, then for any $a>0$ the probability that $\left|\xi \right|\ge a$ is bounded above by $E\left|\xi \right|/a$, that is,

$P\left(\left|\xi \right|>a\right)\le \frac{E\left|\xi \right|}{a}.$

Hence taking $a={\epsilon }^{1/2}$, we see that perturbation method’s strong convergence of order $m=1$ is

$P\left(\left|x_{\epsilon }\left(t\right)-\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\right|\ge {\epsilon }^{1/2}\right)\le K{\epsilon }^{1/2},$

or, equivalently,

$P\left(\left|x_{\epsilon }\left(t\right)-\left(\bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)\right)\right|<{\epsilon }^{1/2}\right)\ge 1-K{\epsilon }^{1/2}.$

This shows that the error at a fixed point in $\left[\mathrm{0,}tfinal\right]$ is small with probability close to 1.

4.8. Conclusion

Stochastic approximations for the process $\stackrel{\dot{}}{x}\left(t\right)$ with parametric noise can be used to analyze aspects of noise. The result of the analysis is an approximation of the form $x\left(t\right)\approx \bar{x}\left(t\right)+\sqrt{\epsilon }\tilde{x}\left(t\right)+o\left(\sqrt{\epsilon }\right)$ for $x\left(t\right)$ of the system. This can be used to evaluate the impact of parametric noise in the neural network. First of all, the system can be averaged. The system for $\bar{x}$ may or may not be analytically solvable but we develop numerical method for its solution. Second, the next order term solves a linear system forced by a Gaussian process, whose statistics depends on the nature of noise in the model.

Our systems are characterized by some system components which combine very fast and very slow behavior. These systems require adaptable step-size, as only in certain phases they require very small step size. It is important to use integration method that allows an efficient step size control. A system is called stiff when integrated with an explicit algorithm and a local error tolerance ${10}^{-n}$, the step size of the algorithm is forced down to below a value indicated by the local error estimate due to constraints imposed on it by the limited size of the numerical stable region. Ode15s is a variable-order solver based on numerical differentiations formulas, optionally uses the backward differentiations formulas (also known a Gear’s method) like ode113, ode15s is multi-step solver. If one suspects that the problem is stiff or if ode45 fails or is very inefficient, try ode15s. But this slow vs. fast time scale problem, if we do not use max step size for the ode15s, the method does not give the correct result which is reflected in Figure 1.

Future work will address systems involving noisy neural network and the impact of noise on a node’s information processing capability which is determined by its signal-to-noise ratio which can be estimated by spectral methods.

Acknowledgements

The author would like to thank Prof. Samir K. Bhowmik for his opinion and discussion. The author also thanks the unanimous referees for their careful reading the manuscript and useful suggestions that improved the paper significantly.

Conflicts of Interest

The author declares no conflict of interest regarding the publication of the paper.

References

[1]	Rahman, M. and Welfert, B. (2013) Functional Central Limit Theorem for Markov Processes and Chains. Journal of Probability and Statistical Science (JPSS), 11, 111-127.
[2]	Battal Gazi Karakoc, S. and Zeybek, H. (2016) Solitary Wave Solutions of the GRLW Equation Using Septic B Spline Collocation Method. Applied Mathematics and Computation, 289, 159-172. https://doi.org/10.1016/j.amc.2016.05.021
[3]	Zeybek, H. and Battal Gazi Karakoc, S. (2017) Application of the Collocation Method with B-Splines to the GEW Equation. Electronic Transactions on Numerical Analysis, 46, 77-88.
[4]	Battal Gazi Karakoc, S., Geyikli, T. and Bashana, A. (2013) A Numerical Solution of the Modified Regularized Long Wave MRLW Equation Using Quartic B Splines. TWMS Journal of Applied and Engineering Mathematics, 3, 231-244. https://doi.org/10.1186/1687-2770-2013-27
[5]	Turgut, A.K. and Battal Gazi Karakoc, S. (2018) A Numerical Technique Based on Collocation Method for Solving Modified Kawahara Equation. Journal of Ocean Engineering and Science, 3, 67-75. https://doi.org/10.1016/j.joes.2017.12.004
[6]	Turgut, A.K., Battal Gazi Karakoc, S. and Triki, H. (2016) Numerical Simulation for Treatment of Dispersive Shallow Water Waves with Rosenau KdV Equation. The European Physical Journal Plus, 131, 1-15. https://doi.org/10.1140/epjp/i2016-16356-3
[7]	Bhowmik, S. and Rahman, M. (2019) Stability and Accuracy Analysis of Theta Scheme for a Convolutional Integro-Differential Equation. Differential Equations and Dynamical Systems, 28, 633-646. https://doi.org/10.1007/s12591-019-00476-w
[8]	Kloeden, P.E. and Platen, E. (1999) Numerical Solution of Stochastic Differential Equations, Applications of Mathematics. Vol. 23, Corrected Third Printing, Springer-Verlag, Berlin.
[9]	Karlin, S. and Taylor, H.W. (1975) A First Course in Stochastic Processes. Academic Press, New York. https://doi.org/10.1016/B978-0-08-057041-9.50005-2
[10]	Gikhmann (1973) Differential Equations with Random Functions. AMS Transl. 12.
[11]	Borodin, A.N. and Salminen, P. (2002) Handbook Brownian Motion-Facts and Formulae. 2nd Edition, Probability and Its Applications, Birkhäuser. https://doi.org/10.1007/978-3-0348-8163-0
[12]	Rahman, M. (2018) Asymptotic Estimate of Variance with Applications to Stochastic Differential Equations Arises in Mathematical Neuroscience. Communications in Statistics—Theory and Methods, 27, 289-306. https://doi.org/10.1080/03610926.2017.1303729
[13]	Skorokhod, A.V. (2000) On Randomly Perturbed Linear Oscillating Mechanical Systems. Ukrainian Mathematical Journal, 52, 1483-1495. https://doi.org/10.1023/A:1010392421925
[14]	Skorokhod, A.V., Hoppensteadt, F.C. and Salehi, H. (2002) Random Perturbations Methods with Applications in Science and Engineering, Applied Mathematical Sciences. Vol. 150, Springer-Verlag, New York. https://doi.org/10.1007/b98905
[15]	Bhattacharya, R.N. (1982) Functional Central Limit Theorem for Markov Processes. Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 60, 185-201. https://doi.org/10.1007/BF00531822
[16]	Berman, A. and Plemmons, R.J. (1994) Nonnegative Matrices in the Mathematical Sciences. Classics in Applied Mathematics, SIAM, Philadelphia. https://doi.org/10.1137/1.9781611971262
[17]	Romanovsky, V. (1970) Discrete Markov Chain. Wolters-Noordhoff Publishing, Groningen.
[18]	Gordin, M.I. and Lifsic, B.A. (1978) The Central Limit Theorem for Stationary Processes. Soviet Mathematics—Doklady, 19, 392-394.
[19]	Komory, Y., Sato, Y. and Mitsui, T. (1994) Some Issues in Discrete Approximate Solution for Stochastic Differential Equations. Computers & Mathematics with Applications, 28, 269-278. https://doi.org/10.1016/0898-1221(94)00197-9