A New Approach to Approximate Solutions of Single Time-Delayed Stochastic Integral Equations via Orthogonal Functions

Abstract

This paper proposes a new numerical method for solving single time-delayed stochastic differential equations via orthogonal functions. The basic principles of the technique are presented. The new method is applied to approximate two kinds of stochastic differential equations with additive and multiplicative noise. Excellence computational burden is achieved along with a $O(h^2)$ convergence rate, which is better than former methods. Two examples are examined to illustrate the validity and efficiency of the new technique.

1. Introduction

Many real-world problems are mathematically formulated in terms of differential equations (DEs), integro-differential equations (IDEs) or partial differential equations (PDEs). Symmetry plays an important role in the study of many types of these equations [1,2,3]. Stochastic differential equations (SDEs) are often adopted to model the time evolution of systems in the areas of biology, physics, economics, and engineering, among others. SDEs have an important role in explaining some symmetry phenomena, namely symmetry breaking in molecular vibration [4]. Due to physical constraints, it is often necessary to introduce time delays into the equations in order to take into account the fact that systems’ response to inputs or perturbations is not instantaneous, leading to the so-called stochastic delayed differential equations (SDDEs). Indeed, energy transfer or signal transmission might be the origin of delay terms. Recurrently, the calculation of the reaction of time-delay systems turns out to be incredibly demanding. However, regardless of the challenges, addressing such issues ends up being unavoidable.

Solving SDDEs is an open issue, and different numerical methods have been proposed. Among them, orthogonal functions can reduce the SDDEs to a linear system of algebraic equations, whose solution describes approximately the behavior of the original system.

In the follow-up, we consider:

$$\left\{\begin{array}{c}j\left(\zeta \right)=j_0\left(\zeta \right)+\underset{0}{\overset{\zeta }{\int }}{k}_1\left(s,\zeta \right)·j(s-\tau )ds+\underset{0}{\overset{\zeta }{\int }}{k}_2\left(s,\zeta \right)·j(s-\tau )dW\left(s\right),\hspace{1em}\hspace{1em}\zeta \in \left[0,\zeta \right),\hfill \\ j\left(\zeta \right)=\mu \left(\zeta \right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\zeta \in \left[-\tau ,0\right),\hfill \end{array}\right.$$

(1)

where $j\left(\zeta \right)$, $j_0\left(\zeta \right)$, $k_1(s,\zeta )$, $k_2(s,\zeta )$, and the 1-D Wiener process $W\left(\zeta \right)$, for $s,\zeta \in [0,Z)$, are stochastic processes defined on the same probability space $(\Omega ,ϝ,P)$, with a filtration $F_{\zeta }$. These satisfy the usual conditions, namely: (i) right-continuous filtration ${\left(F_{\zeta }\right)}_{\zeta ⩾0}$; (ii) each $F_{\zeta },\zeta ⩾0$ contains all P-null sets in F; and (iii) $j\left(\zeta \right)$ is an unknown function, and ${\int }_0^{\zeta }{k}_2(s,\zeta )·j(s-\tau )dW\left(s\right)$ is a stochastic integral that is interpreted in the Itô sense [5].

We are interested in the response to the initial function $j\left(\zeta \right)=\mu \left(\zeta \right)$ or the parameters contained in the SDDEs. In fact, extra data are required to deal with a system of delayed differential equations (DDEs). Since the derivative of Equation (1) depends upon the solution at former time $\zeta -\tau$, it is essential to give an initial history function to determine the value of the solution before time $\zeta =0$. In many usual models, the history is a constant vector. However, non-constant history functions are experienced regularly. The initial time is a jump derivative discontinuity in most problems. Moreover, any solution or derivative discontinuity in the history function at preceding points to the initial time should be dealt with suitably since such discontinuities can spread to after-times.

Numerical analysis of SDDEs based on the Euler–Maruyama scheme was suggested by Buckwar [6]. Triangular functions (TFs) were proposed for the analysis of dynamic systems by Deb et al. [7]. Numerical and computational methods for solving SDEs based on TF and block pulse function (BPF) methods were presented by Khodabin et al. [8]. Operational matrices of BPFs for Taylor series and Bernstein polynomials were computed by Marzban et al. [9] and by Behroozifar et al. [10].

In this paper, we will be keen on obtaining SDDEs approximations via the TFs method. The new approach leads to $O(h^2)$ convergence rate, which is better than alternative methods [11]. The rest of the paper is organized as follows. In Section 2, the basics of TFs and operational delayed matrices are introduced. In Section 3, the new method is presented. In Section 4, two numerical examples are studied. In Section 5, the main conclusion remarks are given.

2. Basic Concepts of Triangular Functions

The TFs originate from a simple dissection of BPFs, where ${\varphi }_i\left(\zeta \right)=T_i^1\left(\zeta \right)+T_i^2\left(\zeta \right)$ is the i-th BPF such that:

$$\begin{gathered} T_i^1\left(\zeta \right)=\left\{\begin{array}{c}1-\frac{\zeta -ih}{h}\hspace{1em}\hspace{1em}ih⩽\zeta <(i+1)h,\hfill \\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{otherwise}.\hfill \end{array}\right. \\ {T}_i^2\left(\zeta \right)=\left\{\begin{array}{c}\frac{\zeta -ih}{h},\hspace{1em}ih⩽\zeta <(i+1)h,\hfill \\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{otherwise},\hfill \end{array}\right. \end{gathered}$$

where $h=\frac{Z}{m}$ with $Z=1$ is considered herein in the interval $\zeta \in [0,Z)$. Therefore, some basic properties, such as orthogonality, finiteness, disjointedness and orthonormality are shared by both.

Two m-set TF vectors are defined in $[0,Z)$ as $T^1\left(\zeta \right)={\left[T_0^1\left(\zeta \right),T_1^1\left(\zeta \right),\dots ,T_{m-1}^1\left(\zeta \right)\right]}^T$, $T^2\left(\zeta \right)={\left[T_0^2\left(\zeta \right),T_1^2\left(\zeta \right),\dots ,T_{m-1}^2\left(\zeta \right)\right]}^T,$ and the 1-D TF vector is:

$$T\left(\zeta \right)={\left[T_i^1\left(\zeta \right),T_i^2\left(\zeta \right)\right]}^T,$$

(2)

where ${\left[\dots \right]}^T$ denotes transposition for $i=0,1,\dots ,m-1$. From the above representation, it follows that:

$$\begin{gathered} T\left(\zeta \right)·T^T\left(\zeta \right)=diag\left(T\left(\zeta \right)\right)=\widehat{T}\left(\zeta \right), \\ T\left(\zeta \right)·T^T\left(\zeta \right)·J\simeq \tilde{J}·T\left(\zeta \right), \\ {T}^T\left(\zeta \right)·B·T\left(\zeta \right)\simeq \stackrel{⌢}{B}·T\left(\zeta \right), \end{gathered}$$

such that $\widehat{T}\left(\zeta \right)$ is a $2m\times 2m$ diagonal matrix, and $\stackrel{⌢}{B}$ is a 2m-vector with elements equal to the diagonal entries of matrix B. Moreover,

$$\begin{array}{c}\hfill \underset{0}{\overset{\zeta }{\int }}{T}_i^p\left(s\right)·T_j^q\left(s\right)ds={\delta }_{ij}\times \left\{\begin{array}{c}\frac{h}{3}\hspace{1em}\hspace{1em}\hspace{1em}p=q\in \left\{1,2\right\},\hfill \\ \frac{h}{6}\hspace{1em}\hspace{1em}\hspace{1em}p\ne q.\hfill \end{array}\right.\end{array}$$

There are some details that make the TF-based methods more efficient than others. Indeed, determining the coefficient vectors of j₀(ζ) in the TF method needs only samples instead of the usual integration and scaling. Another option is a TF piece-wise linear solution, which in most cases leads to a smaller error than the staircase solution of the BPF piece-wise constant. A square integrable function j₀(ζ) over [0, Z) can be broadened into an m-set TF series as:

$$j\left(\zeta \right)\simeq \widehat{j}\left(\zeta \right)\simeq \sum _{i=0}^{m-1}{J}_i^1·T_i^1\left(\zeta \right)+\sum _{i=0}^{m-1}{J}_i^2·T_i^2\left(\zeta \right)=\left[\begin{array}{c}{J}^1\hfill \\ J^2\hfill \end{array}\right]\left[T^1\left(\zeta \right),T^2\left(\zeta \right)\right]={\left[J^T·T\left(\zeta \right)\right]}_{2m\times 2m}\zeta \in \left[0,Z\right),$$

(3)

where $J_i^1=j(i·h)$ and $J_i^2=j(\left(i+1\right)·h)$ for $i=0,1,\dots ,m-1$, and $J^1,J^2$ are coefficient vectors. The operational matrix for integration of TFs, for $i=0,1,\dots ,m-1$, is derived as upper triangular matrices:

$$\underset{0}{\overset{\zeta }{\int }}T\left(s\right)ds=\left[\begin{array}{c}\underset{0}{\overset{\zeta }{\int }}{T}_1\left(s\right)ds\hfill \\ \underset{0}{\overset{\zeta }{\int }}{T}_2\left(s\right)ds\hfill \end{array}\right]=\left[\begin{array}{c}{P}_1·T^1\left(\zeta \right)+P_2·T^2\left(\zeta \right)\hfill \\ P_1·T^1\left(\zeta \right)+P_2·T^2\left(\zeta \right)\hfill \end{array}\right]=P·T\left(\zeta \right),$$

where

$$P=\left(\begin{array}{cc}{P}_1& P_2\\ P_1& P_2\end{array}\right)$$

and

$$P_1=\frac{h}{2}{\left(\begin{array}{ccccc}0& 1& 1& \cdots & 1\\ 0& 0& 1& \cdots & 1\\ 0& 0& 0& \cdots & 1\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 0\end{array}\right)}_{m\times m}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{P}_2=\frac{h}{2}{\left(\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 0& 1& 1& \cdots & 1\\ 0& 0& 1& \cdots & 1\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1\end{array}\right)}_{m\times m},$$

with

$$\underset{0}{\overset{\zeta }{\int }}T\left(s\right)·dW\left(s\right)\simeq P_S·T\left(\zeta \right),$$

where

$$P_S=\left(\begin{array}{cc}{P}_S^1& P_S^1\\ \\ P_S^2& P_S^2\end{array}\right),$$

$$\begin{array}{c}{P}_S^1={\left[\begin{array}{ccccc}\alpha \left(0\right)& \beta \left(0\right)& \beta \left(0\right)& \cdots & \beta \left(0\right)\\ 0& \alpha \left(1\right)& \beta \left(1\right)& \cdots & \beta \left(1\right)\\ 0& 0& \alpha \left(2\right)& \cdots & \beta \left(2\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & \beta \left(m-2\right)\\ 0& 0& 0& \cdots & \alpha \left(m-1\right)\end{array}\right]}_{m\times m},\hfill \end{array}$$

$$\begin{array}{c}{P}_S^2={\left[\begin{array}{ccccc}\gamma \left(0\right)& \rho \left(0\right)& \rho \left(0\right)& \cdots & \rho \left(0\right)\\ 0& \gamma \left(1\right)& \rho \left(1\right)& \cdots & \rho \left(1\right)\\ 0& 0& \gamma \left(2\right)& \cdots & \rho \left(2\right)\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & \rho \left(m-2\right)\\ 0& 0& 0& \cdots & \gamma \left(m-1\right)\end{array}\right]}_{m\times m},\hfill \\ \\ \alpha \left(i\right):=\left(i+1\right)\left[W\left(\left(i+0.5\right)h\right)-W\left(ih\right)\right]-\frac{1}{h}\underset{i·h}{\overset{\left(i+0.5\right)h}{\int }}sdW\left(s\right),\hfill \\ \beta \left(i\right):=\left(i+1\right)\left[W\left(\left(i+1\right)·h\right)-W(i·h)\right]-\frac{1}{h}\underset{i·h}{\overset{\left(i+1\right)h}{\int }}sdW\left(s\right),\hfill \\ \gamma \left(i\right):=-i\left[W\left(\left(i+0.5\right)h\right)-W\left(ih\right)\right]-\frac{1}{h}\underset{ih}{\overset{\left(i+0.5\right)h}{\int }}sdW\left(s\right),\hfill \\ \rho \left(i\right):=-i\left[W\left(\left(i+1\right)·h\right)-W(i·h)\right]-\frac{1}{h}\underset{i·h}{\overset{\left(i+1\right)·h}{\int }}sdW\left(s\right).\hfill \end{array}$$

A function of two variables $k\left(s,\zeta \right)$ can be expanded with respect to TFs as $k\left(s,\zeta \right)=T^T\left(\zeta \right)KT\left(s\right),$ where $T\left(s\right)$ and $T\left(\zeta \right)$ are $2m_1-\mathrm{D}$ and $2m_2-\mathrm{D}$ triangular vectors, and K is a $2m_1\times 2m_2$ coefficient matrix of TFs. For simplicity, we put $m_1=m_2=m$. Therefore, K can be written as:

$$K={\left(\begin{array}{cc}{\left(K_{11}\right)}_{m\times m}& {\left(K_{12}\right)}_{m\times m}\\ {\left(K_{21}\right)}_{m\times m}& {\left(K_{22}\right)}_{m\times m}\end{array}\right)}_{2m\times 2m},$$

where $K_{11},K_{12},K_{21},K_{22}$ can be computed by sampling $k\left(s,\zeta \right)$ at points $s_i$ and ${\zeta }_i$ such that $s_i=t_i=i·h$ for $i=0,1,\dots ,m-1$. Therefore, the following approximations can be obtained:

$$\begin{array}{c}{\left(K_{11}\right)}_{i,j}=k\left(s_i,{\zeta }_j\right)\left(i,j=0,1,\dots ,m-1\right),\hfill \\ {\left(K_{12}\right)}_{i,j}=k\left(s_i,{\zeta }_j\right)\left(i=0,1,\dots ,m-1\right),\left(j=1,\dots ,m\right),\hfill \\ {\left(K_{21}\right)}_{i,j}=k\left(s_i,{\zeta }_j\right)\left(i=1,\dots ,m\right)\left(j=0,1,\dots ,m-1\right),\\ {\left(K_{22}\right)}_{i,j}=k\left(s_i,{\zeta }_j\right)\left(i,j=1,\dots ,m\right),\hfill \\ \hfill \end{array}$$

and

$$k_l(s_i,{\zeta }_j)=\frac{1}{h^2}\underset{(i-1)·h}{\overset{ih}{\int }}\underset{(j-1)h}{\overset{jh}{\int }}{k}_l(s,\zeta )d\zeta ds,l=1,2.$$

As noticed in [12],

$$j_0(\zeta -\tau )\cong j_0^T·T(\zeta -\tau )=j_0^T\times \left[\begin{array}{c}{T}^1(\zeta -\tau )\hfill \\ T^2(\zeta -\tau )\hfill \end{array}\right],$$

(4)

where $J_0^T=\left[J_0^1,J_0^2,\dots ,J_0^m\right]$ is the coefficient vector defined as:

$$J_0^i=\frac{1}{h}\underset{(i-1)·h}{\overset{ih}{\int }}{j}_0\left(\zeta -\tau \right)d\zeta .$$

The delay function $T(\zeta -\tau )$ in (3) is the transfer of $T\left(\zeta \right)$ defined in (1) over the time axis by $(-\tau )$. The general definition is given by:

$$T(\zeta -\tau )={\left[\begin{array}{c}{T}^1(\zeta -\tau )\hfill \\ T^2(\zeta -\tau )\hfill \end{array}\right]}_{2m\times 1}={\left[\begin{array}{c}{Q}_{m\times m}·T^1\left(\zeta \right)\hfill \\ Q_{m\times m}·T^2\left(\zeta \right)\hfill \end{array}\right]}_{2m\times 1}={\stackrel{⌢}{Q}}_{2m\times 2m}·{\left[\begin{array}{c}{T}^1\left(\zeta \right)\hfill \\ T^2\left(\zeta \right)\hfill \end{array}\right]}_{2m\times 1},\hspace{1em}\zeta >\tau$$

where ${\stackrel{⌢}{Q}}_{2m\times 2m}$ is defined as the Kronecker delta:

$$\begin{array}{c}\hfill {\stackrel{⌢}{Q}}_{2m\times 2m}\end{array}=I_l\otimes {\left[Q\right]}_{m\times m}={\left(\begin{array}{cc}{\left[Q\right]}_{m\times m}& \overline{0}\\ \overline{0}& {\left[Q\right]}_{m\times m}\end{array}\right)}_{2m\times 2m}\hspace{1em}\hspace{1em}i=1,2,\dots ,l,$$

where $I_l$ is the l-D identity matrix, and ⊗ denotes the Kronecker product [13], while Q is the delay operational matrix of the TF domain. We can derive the delay matrix as follows:

$$\begin{array}{c}{T}_i^1(\zeta -\tau )=T_i^1(\zeta -kh),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}k=1,2,\dots ,m-1,\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\left\{\begin{array}{c}1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(i-1\right)\frac{Z}{m}⩽\zeta -kh<i\frac{Z}{m},\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{elsewhere},\hspace{1em}\hspace{1em}\hfill \end{array}\right.\hspace{1em}\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\left\{\begin{array}{c}1\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\left(i+k-1\right)\frac{Z}{m}⩽\zeta <\left(i+k\right)\frac{Z}{m},\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\mathit{elsewhere},\hspace{1em}\hspace{1em}\hfill \end{array}\right.\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=T_{i+k}^1\left(\zeta \right).\hfill \end{array}$$

Thus, one has:

$$Q={\left[\begin{array}{cc}\hspace{1em}{\left[\begin{array}{c}0\\ ⋮\\ 0\end{array}\right]}_{k\times 1}& \hspace{1em}{\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]}_{k\times k}\\ {\left[\begin{array}{c}0\end{array}\right]}_{1\times 1}& {\left[\begin{array}{ccc}0& \cdots & 0\end{array}\right.}_{1\times k}\end{array}\right]}_{m\times m}.$$

3. Problem Statement

Many systems are impacted by positive and negative feedback. Such mechanisms push the system to a new state of equilibrium or back to its primary state [14].

Consider a simple delayed feedback by $j^{\prime }\left(\zeta \right)=\alpha ·j(\zeta -\tau ),$ where $\alpha$ is a real value and $\tau \ge 0$. The cases $\alpha >0$ and $\alpha <0$ correspond to negative and positive feedback, respectively.

When $\tau =0$, we recover a simple ODE. We prescribe $j\left(\zeta \right)$ for $-\tau \le \zeta <0$ as initial data for the upper DE. If such an ODE involves a Gaussian white noise process $\xi \left(\zeta \right)=\frac{dW\left(\zeta \right)}{d\zeta }$, then there will be two states, with additive noise $\alpha +b\xi \left(\zeta \right)$ and multiplicative noise $j(\zeta -\tau )+b·\xi \left(\zeta \right)$, in which b is the noise intensity. For further details, please see [15].

3.1. Solving SDDEs Based on Triangular Function Operational Matrices with Additive Noise

Below, we consider a linear stochastic Volterra integral equation with single time delay and intend to obtain the $j\left(t\right)$ coefficient based upon a TF in the interval $\zeta \in \left[0,Z\right]$:

$$\left\{\begin{array}{c}j\left(\zeta \right)=j_0\left(\zeta \right)+\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta )·j(s-\tau )ds+\underset{0}{\overset{t}{\int }}{k}_2(s,\zeta )dW\left(s\right),\tau \in \left[0,Z),\right.\hfill \\ j\left(0\right)=j_0,\hfill \\ j\left(\zeta \right)=\mu \left(\zeta \right),-\tau \le \zeta <0,\hfill \end{array}\right.$$

(5)

where $j_0$ is a constant specified vector, $k_1(s,\zeta )$ and $k_2(s,\zeta )$ are known matrices, and $\mu \left(\zeta \right)$ is an arbitrary, initial-history-known function. We approximate the $j\left(\zeta \right),j_0\left(\zeta \right),k_1(s,\zeta )$ and $k_2(s,\zeta )$ functions by TFs as mentioned below from (6)–(9):

$$j\left(\zeta \right)\approx T^T\left(\zeta \right)·J=J^T·T\left(\zeta \right),$$

(6)

$$j_0\left(\zeta \right)\approx T^T\left(\zeta \right)·J_0={J_0}^T·T\left(\zeta \right),$$

(7)

$$\left\{\begin{array}{c}{k}_1(s,\zeta )\simeq T^T\left(\zeta \right)·K_1·T\left(s\right),\hfill \\ k_2(s,\zeta )\simeq T^T\left(\zeta \right)·K_2·T\left(s\right),\hfill \end{array}\right.$$

(8)

$$j(\zeta -\tau )=\left\{\begin{array}{c}\mu (\zeta -\tau ),-\tau \le \zeta -\tau <0,\hfill \\ J^T·T(\zeta -\tau )=T^T(\zeta -\tau )·J,0\le \zeta -\tau <T-\tau ,\hfill \end{array}\right.$$

(9)

where the vectors $J,J_0$ and matrices $K_1,K_2$ are TF coefficients of $j\left(\zeta \right),j_0\left(\zeta \right),k_1(s,\zeta )$ and $k_2(s,\zeta )$, respectively. Replacing the above approximations into (5), we arrive at:

$$\begin{array}{c}{T}^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+\underset{0}{\overset{\tau }{\int }}{T}^T\left(\zeta \right)·K_1·T\left(s\right)·\mu (s-\tau )ds+\underset{\tau }{\overset{\zeta }{\int }}{T}^T\left(\zeta \right)·K_1·T\left(s\right)·T^T(s-\tau )·Jds\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\underset{0}{\overset{\zeta }{\int }}{T}^T\left(\zeta \right)·K_2·T\left(s\right)\hspace{1em}dW\left(s\right),\hfill \end{array}$$

where

$$T^T(s-\tau )={\left(\stackrel{⌢}{Q}·T\left(s\right))\right.}^T=T^T\left(s\right)·{\stackrel{⌢}{Q}}^T,$$

where ${\stackrel{⌢}{Q}}_{2m\times 2m}$ is defined as the Kronecker delta (p. 5), and

$$\begin{array}{c}\hfill \begin{array}{c}{T}^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+T^T\left(\zeta \right)·K_1\left(\underset{0}{\overset{\tau }{\int }}T\left(s\right)·\mu (\zeta -\tau )ds\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+T^T\left(\zeta \right)·K_1·\left(\underset{\tau }{\overset{\zeta }{\int }}T\left(s\right)·T^T\left(s\right)·{\widehat{Q}}^T·Jds\right)\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+T^T\left(\zeta \right)·K_2\left(\underset{0}{\overset{\zeta }{\int }}T\left(s\right)dW\left(s\right)\right).\hfill \end{array}\end{array}$$

Using previous relations:

$$T^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+T^T\left(\zeta \right)·K_1·I_1+T^T\left(\zeta \right)·K_1·\left(\underset{\tau }{\overset{\zeta }{\int }}\hspace{1em}\tilde{Z}·T\left(s\right)ds\right)+T^T\left(\zeta \right)·K_2·I_2,$$

where $\tilde{Z}=diag({\stackrel{⌢}{Q}}^T·J)$ is a 2m × 2m diagonal matrix. Notice that

$$\begin{array}{c}{I}_1=\underset{0}{\overset{\tau }{\int }}T\left(s\right)\mu (s-\tau )ds=\left[\begin{array}{c}\underset{0}{\overset{\tau }{\int }}{T}^1\left(s\right)·\mu (s-\tau )ds\hfill \\ \underset{0}{\overset{\tau }{\int }}{T}^2\left(s\right)·\mu (s-\tau )ds\hfill \end{array}\right]\hfill \\ \hspace{1em}=\left[\begin{array}{c}\underset{0}{\overset{\tau }{\int }}\left(\left(1+i\right)-\frac{s}{h}\right)·\mu (s-\tau )ds\hfill \\ \underset{0}{\overset{\tau }{\int }}\left(\frac{s}{h}-i\right)·\mu (s-\tau )ds\hfill \end{array}\right]=\left[\begin{array}{c}\left(1+i\right)\underset{0}{\overset{\tau }{\int }}\mu (s-\tau )ds-\frac{1}{h}\underset{0}{\overset{\tau }{\int }}s·\mu (s-\tau )ds\hfill \\ \hspace{1em}\frac{1}{h}\underset{0}{\overset{\tau }{\int }}s·\mu (s-\tau )ds-i\underset{0}{\overset{\tau }{\int }}\mu (s-\tau )ds\hfill \end{array}\right],\hfill \end{array}$$

and

$$\begin{array}{c}{I}_2=\underset{0}{\overset{\tau }{\int }}T\left(s\right)·\mu (s-\tau )dW\left(s\right)=\left[\begin{array}{c}\underset{0}{\overset{\tau }{\int }}{T}^1\left(s\right)·\mu (s-\tau )dW\left(s\right)\hfill \\ \underset{0}{\overset{\tau }{\int }}{T}^2\left(s\right)·\mu (s-\tau )dW\left(s\right)\hfill \end{array}\right]\hfill \\ =\left[\begin{array}{c}\underset{0}{\overset{\tau }{\int }}\left((1+i)-\frac{s}{h}\right)·\mu (s-\tau )dW\left(s\right)\hfill \\ \underset{0}{\overset{\tau }{\int }}\left(\frac{s}{h}-i\right)·\mu (s-\tau )dW\left(s\right)\hfill \end{array}\right]\hfill \\ =\left[\begin{array}{c}(1+i)\underset{0}{\overset{\tau }{\int }}\mu (s-\tau )dW\left(s\right)-\frac{1}{h}\underset{0}{\overset{\tau }{\int }}s·\mu (s-\tau )dW\left(s\right)\hfill \\ \frac{1}{h}\underset{0}{\overset{\tau }{\int }}s·\mu (s-\tau )dW\left(s\right)-i\underset{0}{\overset{\tau }{\int }}\mu (s-\tau )dW\left(s\right)\hfill \end{array}\right].\hfill \end{array}$$

Then,

$$T^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+T^T\left(\zeta \right)·K_1·I_1+T^T\left(\zeta \right)·K_1·{\stackrel{⌢}{Q}}^T·J·P·T^T\left(\zeta \right)·K_2·I_2,$$

in which $K_1·{\stackrel{⌢}{Q}}^T·J·P$ are 2m × 2m matrices, and

$$\begin{array}{c}\hfill T^T\left(\zeta \right)·\left(K_1·\tilde{Z}·P\right)·T\left(\zeta \right)=B^T·T\left(\zeta \right)=T^T\left(\zeta \right)·B,\hspace{1em}\end{array}$$

where B is a 2m-D vector with components equal to the diagonal entries of the matrix $\left(K_1·{\stackrel{⌢}{Q}}^T·J·P\right)$, B can be written as B = Π · J and $\Pi =K_1·P^T·{\stackrel{⌢}{Q}}^T$. Finally, we arrive at:

$$\begin{array}{c}\hfill T^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+T^T\left(\zeta \right)·K_1·I_1+T^T\left(\zeta \right)·\Pi ·J+T^T\left(\zeta \right)·K_2·I_2.\end{array}$$

By multiplying both sides in T(ζ) and replacing ≃ with =, we have:

$$\begin{array}{c}\hfill \left(I-\Pi \right)·J\simeq \left(J_0+K_1·I_1+K_2·I_2\right),\end{array}$$

which we denote by J₀₀ = J₀ + K₁ · I₁ + K₂ · I₂, R = I − Π. The solution J is obtained by solving the algebraic equation:

$$\begin{array}{c}\hfill R·J=J_{00}.\end{array}$$

(10)

3.2. Solving SDDEs Based on Triangular Function Operational Matrices with Multiplicative Noise

In the following linear stochastic Volterra integral equations (SVIEs) with single time delay, our purpose is to obtain the TF coefficients of j(ζ) in the interval $\zeta \in [0,Z)$:

$$\left\{\begin{array}{c}j\left(\zeta \right)=j_0\left(\zeta \right)+\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta )·j(s-\tau )ds+\underset{0}{\overset{\zeta }{\int }}{k}_2(s,\zeta )·j(s-\tau )dW\left(s\right),\hspace{1em}\hspace{1em}\tau \in \left[0,Z),\right.\hfill \\ j\left(0\right)=j_0,\hfill \\ \\ j\left(\zeta \right)=\mu \left(\zeta \right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-\tau \le \zeta <0.\hspace{1em}\hfill \end{array}\right.$$

(11)

Substituting function approximations (6)–(9) into (11), we arrive at:

$$\begin{array}{c}{T}^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+\underset{0}{\overset{\tau }{\int }}{T}^T\left(\zeta \right)·K_1·T\left(s\right)·J(s-\tau )ds+\underset{\tau }{\overset{\zeta }{\int }}{T}^T\left(\zeta \right)·K_1·T\left(s\right)·T^T(s-\tau )·Jds\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\underset{0}{\overset{\tau }{\int }}{T}^T\left(\zeta \right)·K_2·T\left(s\right)·J(s-\tau )dW\left(s\right)+\underset{\tau }{\overset{\zeta }{\int }}{T}^T\left(\zeta \right)·K_2·T\left(s\right)·T^T(s-\tau )·JdW\left(s\right),\hfill \end{array}$$

(12)

and

$$\begin{array}{c}{T}^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+T^T\left(\zeta \right)·K_1·\left(\underset{0}{\overset{\tau }{\int }}T\left(s\right)·\mu (\zeta -\tau )ds\right)\hfill \\ \hspace{1em}\hspace{1em}+T^T\left(\zeta \right)·K_1·\left(\underset{\tau }{\overset{\zeta }{\int }}T\left(s\right)·T^T\left(s\right){\widehat{Q}}^T·Jds\right)\hfill \\ \hspace{1em}\hspace{1em}+T^T\left(\zeta \right)·K_2·\left(\underset{0}{\overset{\tau }{\int }}T\left(s\right)·\mu (\zeta -\tau )dW\left(s\right)\right)\hfill \\ \hspace{1em}\hspace{1em}+T^T\left(\zeta \right)·K_2·\left(\underset{\tau }{\overset{\zeta }{\int }}T\left(s\right)·T^T\left(s\right)·{\widehat{Q}}^T·JdW\left(s\right)\right)·\hfill \end{array}$$

Therefore,

$$\begin{array}{c}{T}^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+T^T\left(\zeta \right)·K_1·I_1+T^T\left(\zeta \right)·K_1·\tilde{Z}·\left(\underset{\tau }{\overset{\zeta }{\int }}T\left(s\right)ds\right)\hfill \\ \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+T^T\left(\zeta \right)·K_2·I_2+T^T\left(\zeta \right)·K_2·\tilde{Z}·\left(\underset{\tau }{\overset{\zeta }{\int }}\hspace{1em}T\left(s\right)dW\left(s\right)\right),\hfill \end{array}$$

where $I_1$, $I_2$ and both integrals are as considered in the previous section. Then,

$$\begin{array}{cc}\hfill & T^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+T^T\left(\zeta \right)·K_1·I_1+\hfill \\ & T^T\left(\zeta \right)·K_1·{\widehat{Q}}^T·J·P·T\left(\zeta \right)+T^T\left(\zeta \right)·K_2·I_2+T^T\left(\zeta \right)·K_2·{\widehat{Q}}^T·J·P_S·T\left(\zeta \right).\hfill \end{array}$$

(13)

Thus, $\widehat{Q},\Pi ,{\Pi }_S$ are $2m\times 2m$ matrices, and

$$T^T\left(\zeta \right)·J\simeq T^T\left(\zeta \right)·J_0+T^T\left(\zeta \right)·K_1·I_1+T^T\left(\zeta \right)·\Pi ·J+T^T\left(\zeta \right)·K_2·I_2+T^T\left(\zeta \right)·{\Pi }_S·J,$$

(14)

in which

$$\begin{array}{c}{B}_1=K_1·P^T·{\tilde{Z}}^T=\left[K_1·P^T·{\widehat{Q}}^T\right]J=\Pi ·J\Pi =K_1·P^T·{\widehat{Q}}^T,\hfill \\ \\ B_2=K_2·{P_S}^T·{\tilde{Z}}^T=\left[K_2·{P_S}^T·{\widehat{Q}}^T\right]J={\Pi }_S·J{\Pi }_S=K_2{P_S}^T·{\widehat{Q}}^T.\hfill \end{array}$$

(15)

Lastly, by multiplying both sides in $T\left(\zeta \right)$, we arrive at:

$$\left(I-\Pi -{\Pi }_S\right)·J\simeq \left(J_0+K_1·I_1+K_2·I_2\right),$$

(16)

and replacing ≃ with =, we have

$$R·J=J_{00}.$$

(17)

Equations (10) and (12) are linear systems of equations with lower triangular coefficient matrices that give the approximate triangular coefficients of the unknown stochastic processes $j\left(\zeta \right)$.

This fulfills the assertions.

4. Error Analysis

This part addresses the rate of convergence. The results show a good level of accuracy of order $O\left(h^2\right)$. We utilize Theorem 1 [16] and assert Theorem 2.

Theorem 1

([16]). Suppose that $j_0\left(\zeta \right)$ is an arbitrary real bounded function, which is twice differentiable in the interval $\left[0,1\right)$, ${\overline{j}}_{0,m}\left(\zeta \right)$ correspond to TFs, and $\left|{j^{″}}_0\left(\zeta \right)\right|<M$ for every $\zeta \in [0,1].$ Then, $∥j\left(\zeta \right)-{j^{″}}_0\left(\zeta \right)∥<O\left(h^2\right).$

Theorem 2.

Let $j\left(\zeta \right)$ and $\overline{j}\left(\zeta \right)$ be solutions of Equations (1) and (6), respectively, and let $∥j\left(\zeta \right)∥<C$ and $∥k_i∥<C$ for $i=1,2.$ Then,

$$E\left({∥j\left(\zeta \right)-\overline{j}\left(\zeta \right)∥}^2\right)\le O\left(h^2\right)\zeta \in \left[0,1\right),$$

and

$${sup}_{0\le \zeta <Z}{\left(E{\left(∥\left(j\left(\zeta \right)-\overline{j}\left(\zeta \right)\right)∥\right)}^2\right)}^{\frac{1}{2}}=O\left(h^4\right)\hspace{1em}\hspace{1em}\zeta \in \left[0,1\right].$$

Proof. 

We follow with an error analysis in two manners:

• A: Based on Equation $\left(1\right)$,

$$\begin{array}{c}\hfill \begin{array}{c}j\left(\zeta \right)-\overline{j}\left(\zeta \right)=\left\{j_0\left(\zeta \right)-{\overline{j}}_0\left(\zeta \right)\right\}+\hfill \\ \left\{\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta )·j(s-\tau )ds-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_1(s,\zeta )·\overline{j}(s-\tau )ds\right\}+\hfill \\ \left\{\underset{0}{\overset{\zeta }{\int }}{k}_2(s,\zeta )dW\left(s\right)-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_2(s,\zeta )dW\left(s\right)\right\},\hfill \end{array}\end{array}$$

(18)

and considering the Euclidean norm:

$$\begin{array}{c}\hfill \begin{array}{c}E{∥j\left(\zeta \right)-\overline{j}\left(\zeta \right)∥}^2=E∥\left\{j_0\left(\zeta \right)-{\overline{j}}_0\left(\zeta \right)\right\}+\left\{\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta )·j(s-\tau )ds-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_1(s,\zeta )·\overline{j}(s-\tau )ds\right\}+\hfill \\ {\left\{\underset{0}{\overset{\zeta }{\int }}{k}_2(s,\zeta )dW\left(s\right)-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_2(s,\zeta )dW\left(s\right)\right\}∥}^2,\hfill \end{array}\end{array}$$

(19)

and ${(a+b+c)}^2\le 3(a^2+b^2+c^2),$ we obtain:

$$\begin{array}{c}\le 3\left(E{∥j_0\left(\zeta \right)-{\overline{j}}_0\left(\zeta \right)∥}^2+\right.\hfill \\ E{∥\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta )·j(s-\tau )ds-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_1(s,\zeta )·\overline{j}(s-\tau )ds∥}^2+E{∥\underset{0}{\overset{\zeta }{\int }}{k}_2(s,\zeta )dW\left(s\right)-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_2(s,\zeta )dW\left(s\right)∥}^2.\hfill \end{array}$$

Then, we obtain the last two parts one by one:

$$\begin{array}{c}\hfill \begin{array}{c}{I}_1=E{∥\underset{0}{\overset{\zeta }{\int }}\left(k_1(s,\zeta )·j(s-\tau )-{\overline{k}}_1(s,\zeta )·\overline{j}(s-\tau )\right)ds∥}^2\hfill \\ \\ =E∥\left(\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta )·j(s-\tau )-{\overline{k}}_1(s,\zeta )·j(s-\tau )ds\right)+\hfill \\ {\left(\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_1(s,\zeta )·j(s-\tau )-{\overline{k}}_1(s,\zeta )·\overline{j}(s-\tau )ds\right)∥}^2\hfill \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le 2E{∥\underset{0}{\overset{\zeta }{\int }}\left\{k_1(s,\zeta )-{\overline{k}}_1(s,\zeta )\right\}·j(s-\tau )ds∥}^2+2E{∥\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_1(s,\zeta )·\left\{j(s-\tau )-\overline{j}(s-\tau )\right\}ds∥}^2\hfill \\ \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le C\underset{0}{\overset{\zeta }{\int }}E{∥k_1(s,\zeta )-{\overline{k}}_1(s,\zeta )∥}^{2}ds+C\underset{0}{\overset{\zeta }{\int }}\left(E{∥j(s-\tau )-\overline{j}(s-\tau )∥}^2\right)ds\hfill \\ \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\le C·O\left(h^4\right)+C·O\left(h^4\right)=O\left(h^4\right),\hfill \end{array}\end{array}$$

(20)

and $I_2=E{∥\underset{0}{\overset{\zeta }{\int }}{k}_2(s,\zeta )-{\overline{k}}_2(s,\zeta )dW\left(s\right)∥}^2=\underset{0}{\overset{\zeta }{\int }}E\left({∥k_2(s,\zeta )-{\overline{k}}_2(s,\zeta )∥}^2\right)ds\le O\left(h^4\right),(It\stackrel{⌢}{o}isometry$ ).

Therefore

$$E{∥j\left(\zeta \right)-\overline{j}\left(\zeta \right)∥}^2\le 3(O\left(h^4\right)+O\left(h^4\right)+O\left(h^4\right)),$$

$$E∥j\left(\zeta \right)-\overline{j}\left(\zeta \right)∥\le C.O\left(h^2\right).$$

• B: Based on Equation (1)

$$\begin{array}{c}j\left(\zeta \right)-\overline{j}\left(\zeta \right)=\left\{j_0\left(\zeta \right)-{\overline{j}}_0\left(\zeta \right)\right\}+\left\{\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta )·z(s-\tau )ds-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_1(s,\zeta )·\overline{j}(s-\tau )ds\right\}\hfill \\ \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left\{\underset{0}{\overset{\zeta }{\int }}{k}_2(s,\zeta )·j(s-\tau )dW\left(s\right)-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_2(s,\zeta )·\overline{j}(s-\tau )dW\left(s\right)\right\},\hfill \end{array}$$

and we obtain the Euclidean norm:

$$\begin{array}{c}{∥j\left(\zeta \right)-\overline{j}\left(\zeta \right)∥}^2={∥j_0\left(\zeta \right)-{\overline{j}}_0\left(\zeta \right)∥}^2+{∥\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta )·j(s-\tau )ds-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_1(s,\zeta )·\overline{j}(s-\tau )ds∥}^2\hfill \\ \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{∥\underset{0}{\overset{\zeta }{\int }}{k}_2(s,\zeta )·j(s-\tau )dW\left(s\right)-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_2(s,\zeta )·\overline{j}(s-\tau )dW\left(s\right)∥}^2\hfill \end{array}$$

and

$$\begin{array}{c}E\left({∥j\left(\zeta \right)-\overline{j}\left(\zeta \right)∥}^2\right)\le 3\left(E{∥j_0\left(\zeta \right)-{\overline{j}}_0\left(\zeta \right)∥}^2+E{∥\underset{0}{\overset{\zeta }{\int }}{k}_1(s,\zeta ).j(s-\tau )ds-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_1(s,\zeta ).j(s-\tau )ds∥}^2\right.\hfill \\ \\ \left.+E{∥\underset{0}{\overset{\zeta }{\int }}{k}_2(s,\zeta ).j(s-\tau )dW\left(s\right)-\underset{0}{\overset{\zeta }{\int }}{\overline{k}}_2(s,\zeta ).\overline{j}(s-\tau )dW\left(s\right)∥}^2\right).\hfill \end{array}$$

By using the Itô isometry property [17] and Theorem 1, we have:

$$\begin{array}{c}\le CE\left({∥j_0\left(\zeta \right)-{\overline{j}}_0\left(\zeta \right)∥}^2\right)+6\underset{0}{\overset{\zeta }{\int }}E\left({∥\left(k_1(s,\zeta )-{\overline{k}}_1(s,\zeta )\right)·j(s-\tau )∥}^2\right)ds\hfill \\ \\ +6\underset{0}{\overset{\zeta }{\int }}E\left({∥{\overline{k}}_1(s,\zeta )·(j(s-\tau )-\overline{j}(s-\tau ))∥}^2\right)ds+6\underset{0}{\overset{\zeta }{\int }}E\left({∥\left(k_2(s,\zeta )-{\overline{k}}_2(s,\zeta )\right)·j(s-\tau )∥}^2\right)ds\hfill \\ \\ +6\underset{0}{\overset{\zeta }{\int }}E\left({∥{\overline{k}}_2(s,\zeta )·(j(s-\tau )-\overline{j}(s-\tau ))∥}^2\right)ds.\hfill \end{array}$$

By the problem’s assumption,

$$\begin{array}{c}\le CE\left({∥j_0\left(\zeta \right)-{\overline{j}}_0\left(\zeta \right)∥}^2\right)+6C^2\underset{0}{\overset{\zeta }{\int }}E\left({∥\left(k_1(s,\zeta )-{\overline{k}}_1(s,\zeta )\right)∥}^2\right)ds+6C^2\underset{0}{\overset{\zeta }{\int }}E\left({∥(j\left(s\right)-\overline{j}\left(s\right))∥}^2\right)ds\hfill \\ \\ +6C^2\underset{0}{\overset{\zeta }{\int }}E\left({∥\left(k_2(s,\zeta )-{\overline{k}}_2(s,\zeta )\right)∥}^2\right)ds+6C^2\underset{0}{\overset{\zeta }{\int }}E\left({∥(j(s-\tau )-\overline{j}(s-\tau ))∥}^2\right)ds,\hfill \end{array}$$

and

$$E\left({∥j\left(\zeta \right)-\overline{j}\left(\zeta \right)∥}^2\right)\le O\left(h^2\right),$$

We complete the proof.  □

5. Illustrative Examples

Examples are used to demonstrate the theoretical outcomes stated in Section 3 and Section 4. The related computations are performed in $Matlab$$\left(2015a\right).$

Let us consider the assumptions:

• h: $h=\frac{Z}{m}$ where $Z=1$ in $\zeta \in [0,Z)$.
• ${\zeta }_i$: The number of interval division. $t_i=i·h$ node points for $i=0,1,\dots ,m$.
• $\tau$: Delay parameter $\tau =(k+\alpha )·h$ with nonnegative integer $\alpha =0$ ($0\le \alpha <1$).
• n: Number of iterations.
• p: Dividing of interval $[0,Z)$ in $p=m.z$ equal units. z is a coefficient of $m=k.N$.
• $\sigma$: Standard deviation.
• $j_i$: TF coefficient of the analytic solution.
• ${\overline{j}}_i$: TF coefficient of the approximated method.
• e: Absolute error computed by ${∥e∥}_{\infty }={max}_{1\le i\le m}\left|j_i-{\overline{j}}_i\right|.$
• ${\overline{z}}_e$: Mean of error.
• $s_e$: Standard deviation of error.
• $UB$: Upper bound.
• $LB$: Lower bound.

Time delayed equations are mostly solved analytically in a step-wise procedure, called method of steps, with initial condition given [11].

Example 1.

Consider a linear stochastic Volterra integral equation with additive noise input and constant time delay:

$$\left\{\begin{array}{c}j\left(\zeta \right)=j\left(0\right)+{\int }_0^{\zeta }4\hspace{1em}j(s-0.25)ds+4b{\int }_0^{\zeta }\hspace{1em}dW\left(s\right),\hspace{1em}\hspace{1em}\zeta \in \left[0,1\right),\hfill \\ j\left(\zeta \right)=1,\hfill \\ j\left(0\right)=1,\hfill \end{array}\right.$$

(21)

where $j\left(\zeta \right)$ is an unknown stochastic process on the probability space $(\Omega ,F,P)$, and $W\left(\zeta \right)$ is a Wiener process:

$$j\left(\zeta \right)=\left\{\begin{array}{l}1+4bW\left(\zeta \right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}0\le \zeta <\frac{1}{4}\\ 1+4(\zeta -\frac{1}{4})+4b\left(4\left(\underset{\frac{1}{4}}{\overset{\zeta }{\int }}W(s-\frac{1}{4})ds\right)+W\left(\zeta \right)\right),\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{1}{4}\le \zeta <\frac{2}{4}\\ 1+4\left(\zeta -\frac{1}{4}\right)+8{(\zeta -\frac{1}{2})}^2+\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{2}{4}\le \zeta <\frac{3}{4}\\ 4b\left[16(\zeta -\frac{1}{2})+\left(\underset{\frac{1}{4}}{\overset{\zeta -\frac{1}{4}}{\int }}W(s-\frac{1}{4})ds\right)+4\left(\underset{\frac{1}{2}}{\overset{\zeta }{\int }}W(u-\frac{1}{4})du\right)+4\left(\underset{\frac{1}{4}}{\overset{\frac{1}{2}}{\int }}W(s-\frac{1}{4})ds\right)+W\left(\zeta \right)\right],\\ 1+4\left(\zeta -\frac{1}{4}\right)+8{(\zeta -\frac{1}{2})}^2+\frac{32}{3}{(\xi -\frac{3}{4})}^3+\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{3}{4}\le \zeta <1\\ 4b\left[\begin{array}{c}\frac{4^3}{2}(\zeta -\frac{3}{4})\left(\underset{\frac{1}{4}}{\overset{\zeta -\frac{1}{2}}{\int }}W(s-\frac{1}{4})ds\right)+4^2(\zeta -\frac{3}{4})\left(\underset{\frac{1}{2}}{\overset{\zeta -\frac{1}{4}}{\int }}W(u-\frac{1}{4})du\right)+W\left(\zeta \right)+\hfill \\ 4\left(\underset{\frac{3}{4}}{\overset{\zeta }{\int }}W(v-\frac{1}{4})dv\right)+8(1+2(\zeta -\frac{3}{4}))\left(\underset{\frac{1}{4}}{\overset{\frac{2}{4}}{\int }}W(s-\frac{1}{4})ds\right)+4\left(\underset{\frac{1}{2}}{\overset{\frac{3}{4}}{\int }}W(u-\frac{1}{4})du\right)\hfill \end{array}\right]\end{array}\right.$$

Example 2.

Assume a linear stochastic Volterra integral equation with multiplicative noise input and time-varying delay as follows:

$$\left\{\begin{array}{c}j\left(\zeta \right)=j\left(0\right)+{\int }_0^{\zeta }4j(s-0.25)ds+4b{\int }_0^{\zeta }j(s-0.25)dW\left(s\right),\hspace{1em}\zeta \in \left[0,1\right),\hfill \\ j\left(\zeta \right)=1,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\zeta \in \left[-0.25,0\right),\hfill \\ j\left(0\right)=1,\hfill \end{array}\right.$$

where $j\left(\zeta \right)$ is an unknown stochastic process on the probability space $(\Omega ,F,P)$, and $W\left(\zeta \right)$ is a Wiener process:

$$j\left(\zeta \right)=\left\{\begin{array}{c}1,\hspace{1em}\hspace{1em}0\le \zeta <\frac{1}{4}\hfill \\ \\ 1+4(\zeta -\frac{1}{4})+b\left(W\left(\zeta \right)-W\left(\frac{1}{4}\right)\right),\frac{1}{4}\le \zeta <\frac{2}{4}\hfill \\ \\ 1+4\left(\zeta -\frac{1}{4}\right)+8{(\zeta -\frac{1}{2})}^2+4b\left[\left(\underset{\frac{1}{2}}{\overset{\zeta }{\int }}W(s-\frac{1}{4})ds\right)-\left(\underset{\frac{1}{2}}{\overset{\zeta }{\int }}W\left(s\right)ds\right)+\frac{2}{4}\le \zeta <\frac{3}{4}\right.\hfill \\ \left.+\zeta \left(W\left(\zeta \right)-W\left(\frac{1}{4}\right)\right)\right]+b\left[b\left(\underset{\frac{1}{2}}{\overset{\zeta }{\int }}W(s-\frac{1}{4})dW\left(s\right)\right)-(1+bW\left(\frac{1}{4}\right))W\left(\zeta \right)+W\left(\frac{1}{4}\right)(1+bW\left(\frac{1}{2}\right))\right],\hfill \\ \hfill \\ \hfill \\ 1+4\left(\zeta -\frac{1}{4}\right)+8{\left(\zeta -\frac{1}{2}\right)}^2+\frac{32}{3}{(\zeta -\frac{3}{4})}^3+\frac{3}{4}\le \zeta <1\hfill \\ \hfill \\ 4b\left[\left(3+bW\left(\frac{1}{4}\right)\right)\left(\underset{\frac{3}{4}}{\overset{\zeta }{\int }}W\left(u\right)du\right)-\left(b\left(W\left(\zeta \right)-W\left(\frac{3}{4}\right)\right)+4(\zeta -\frac{3}{4})\right)\left(\underset{\frac{1}{2}}{\overset{\zeta -\frac{1}{4}}{\int }}W\left(s\right)ds\right)-\right.\hfill \\ \left(\underset{\frac{1}{2}}{\overset{\frac{3}{4}}{\int }}W\left(s\right)ds\right)+\left(3+bW\left(\frac{1}{4}\right)\right)\left(\underset{\frac{3}{4}}{\overset{\zeta }{\int }}W\left(u\right)du\right)-\left(\underset{\frac{1}{2}}{\overset{\frac{3}{4}}{\int }}W\left(s\right)ds\right)\hfill \\ -(2+bW\left(\frac{1}{4}\right))\left(\underset{\frac{3}{4}}{\overset{\zeta }{\int }}W(u-\frac{1}{4})du\right)+\left(\underset{\frac{1}{2}}{\overset{\frac{3}{4}}{\int }}W(s-\frac{1}{4})ds\right)+2{\zeta }^2\left(W\left(\zeta \right)-W\left(\frac{1}{4}\right)\right)\hfill \\ +\left(b\left(W\left(\zeta \right)-W\left(\frac{3}{4}\right)\right)+4(\zeta -\frac{3}{4})\right)\left(\underset{\frac{1}{2}}{\overset{\zeta -\frac{1}{4}}{\int }}W(s-\frac{1}{4})ds\right)+b\left(\underset{\frac{3}{4}}{\overset{\zeta }{\int }}uW(u-\frac{1}{4})dW\left(u\right)\right)-\frac{7}{8}W\left(\frac{1}{4}\right)+\hfill \\ \\ \left.\zeta \left(\left(2+bW\left(\frac{1}{2}\right)\right)W\left(\frac{1}{4}\right)-\left(2+bW\left(\frac{1}{4}\right)\right)W\left(\zeta \right)\right)-\left(b\left(W\left(\zeta \right)-W\left(\frac{3}{4}\right)\right)+4(\zeta -\frac{3}{4})\right)\left(\underset{\frac{1}{2}}{\overset{\zeta -\frac{1}{4}}{\int }}W\left(s\right)ds\right)\right]\hfill \\ \hfill \\ +b^2\left(\left(b\left(W\left(\zeta \right)-W\left(\frac{3}{4}\right)\right)+4(\zeta -\frac{3}{4})\right)\left(\underset{\frac{1}{2}}{\overset{\zeta -\frac{1}{4}}{\int }}W(s-\frac{1}{4})dW\left(s\right)\right)+\left(\underset{\frac{1}{2}}{\overset{\frac{3}{4}}{\int }}W(s-\frac{1}{4})dW\left(s\right)\right)-\right.\hfill \\ \left.(2+bW\left(\frac{1}{4}\right))\left(\underset{\frac{3}{4}}{\overset{\zeta }{\int }}W(u-\frac{1}{4})dW\left(u\right)\right)+\left(2W\left(\frac{1}{4}\right)+bW\left(\frac{1}{4}\right)W\left(\frac{1}{2}\right)+\frac{7}{2}b\right)W\left(\zeta \right)-W\left(\frac{1}{2}\right)W\left(\frac{1}{4}\right)\left(2+bW\left(\frac{3}{4}\right)\right)\right).\hfill \end{array}\right.$$

(22)

The numerical results of the mean and standard deviation with 95% confidence intervals for some different values of ${\zeta }_i$ in the points $h=\frac{1}{m}$ are shown in

Table 1. SDDE with additive noise, n = 40, N = 4, k = 2, p = 1040, σ = 0.99, b = 0.85.

${\mathit{\zeta }}_{\mathit{i}}$	${\overline{\mathit{z}}}_{\mathit{e}}$	${\mathit{s}}_{\mathit{e}}$	LB	UB
0	0	0	0	0
$0.125$	0	0	0	0
$0.255$	$-0.00013233$	$0.013495$	$-0.00039169$	$0.00012703$
$0.375$	$-0.00012545$	$0.012793$	$-0.00037133$	$0.00012042$
$0.5$	$-0.00025845$	$0.026356$	$-0.00076499$	$0.00024809$
$0.625$	$-0.00031402$	$0.032023$	$-0.00092948$	$0.00030143$
$0.75$	$-0.00053561$	$0.054619$	$-0.0015854$	$0.00051413$
$0.875$	$-0.0006232$	$0.063551$	$-0.0018446$	$0.00059821$
1	$-0.00071978$	$0.0734$	$-0.0021305$	$0.00069092$

and

Table 2. SDDE with multiplicative noise, n = 40, N = 4, k = 2, p = 1040, σ = 0.99, b = 0.85.

${\mathit{\zeta }}_{\mathit{i}}$	${\overline{\mathit{z}}}_{\mathit{e}}$	${\mathit{s}}_{\mathit{e}}$	LB	UB
0	0	0	0	0
$0.125$	0	0	0	0
$0.25$	$-0.00011832$	$0.012066$	$-0.00035023$	$0.00011358$
$0.375$	$-0.00015062$	$0.01536$	$-0.00044583$	$0.00014459$
$0.5$	$-0.00017632$	$0.01798$	$-0.00052188$	$0.00016925$
$0.625$	$-0.00026521$	$0.027044$	$-0.00078498$	$0.00025457$
$0.75$	$-0.00019342$	$0.019724$	$-0.00057251$	$0.00018567$
$0.875$	$-0.00026007$	$0.026521$	$-0.00076978$	$0.00024964$
1	$-0.00051652$	$0.052672$	$-0.0015289$	$0.00049581$

. Figure 1 and Figure 2 express a trajectory of the approximated solution calculated by the proposed method, along with a trajectory of the analytic solution. It should be noted that due to the high accuracy of the method, in both graphs, the approximated and numerical solutions practically coincide.>

6. Conclusions

In this paper, a new method based on operational matrices of TFs was applied to solve SDDEs with constant time delay. The numerical scheme is computationally simple and revealed good accuracy. An $O\left(h^2\right)$ convergence rate is ensured, contrasting with the $O\left(h\right)$ of other methods. Further research will address numerical solutions to diverse orthogonal functions.