Numerical Solution of Two-Dimensional Nonlinear Stochastic Itô-Volterra Integral Equations by Applying Block Pulse Functions

Guo Jiang, Xiaoyan Sang^*, Jieheng Wu, Biwen Li

School of Mathematics and Statistics, Hubei Normal University, Huangshi, China.

DOI: 10.4236/apm.2019.92004                

Abstract

This paper investigates the numerical solution of two-dimensional nonlinear stochastic Itô-Volterra integral equations based on block pulse functions. The nonlinear stochastic integral equation is transformed into a set of algebraic equations by operational matrix of block pulse functions. Then, we give error analysis and prove that the rate of convergence of this method is efficient. Lastly, a numerical example is given to confirm the method.

Keywords

Block Pulse Functions, Integration Operational Matrix, Stochastic Itô-Volterra Integral Equations

1. Introduction

Two-dimensional stochastic Itô-Volterra integral equations arise from many phenomena in physics and engineering fields [1] . Some different orthogonal basis functions, polynomials and wavelets are used to approximate the solution of two-dimensional Volterra integral equations. For example, block pulse functions, triangular functions, modification of hat functions, Legender polynomials and Haar wavelet and the like (see [2] [3] [4] [5] [6] ).

Especially, Fallahpour et al. [3] introduced the following two-dimensional linear stochastic Volterra integral equation by Haar wavelet

$\begin{array}{c}x\left(t_1,t_2\right)=f\left(t_1,t_2\right)+{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\tilde{k}\left(t_1,t_2,s_1,s_2\right)x\left(s_1,s_2\right)\text{d}{s}_1\text{d}{s}_2\\ +{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\hat{k}\left(t_1,t_2,s_1,s_2\right)x\left(s_1,s_2\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right),\end{array}$ (1)

where $x\left(t_1\mathrm{,}{t}_2\right)$ is unknown and called the solution of the Equation (1), $\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ , $\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ and $f\left(t_1\mathrm{,}{t}_2\right)$ are known functions $\left(t_1\mathrm{,}{t}_2\right)\in \left[\mathrm{0,}{T}_1\right)\times \left[\mathrm{0,}{T}_2\right)$ , $s_1\le t_1\mathrm{,}{s}_2\le t_2$ . $B\left(s_1\right)$ and $B\left(s_2\right)$ are two independent Brownian motions and ${\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)x\left(s_1\mathrm{,}{s}_2\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)$ is the double Itô integral. The authors transformed stochastic Volterra integral equations to algebra equations by Haar wavelet and gave the numerical solutions to the equations. Similarly, Fallahpour et al. [7] obtained a numerical method for two-dimensional linear stochastic Volterra integral equations by block pulse functions.

For nonlinear determinate Volterra integral equations, Maleknejad et al. [8] and Nemati et al. [6] used two-dimensional block pulse functions and Legendre polynomials to solve those respectively. Both Babolian et al. [2] and Maleknejad et al. [9] employed triangular functions to get the numerical solutions. Mirzaee et al. [5] [10] applied modified two-dimensional block pulse functions to approximate the following determinate equation

$f\left(t_1,t_2\right)={\displaystyle {\int }_0^{t_1}}{\displaystyle {\int }_0^{t_2}}\tilde{k}\left(t_1,t_2,s_1,s_2\right){\left[x\left(s_1,s_2\right)\right]}^n\text{d}{s}_2\text{d}{s}_1,\left(t_1,t_2\right)\in \left[0,T_1\right)\times \left[0,T_2\right),$ (2)

where nonlinear term ${\left[x\left(s_1\mathrm{,}{s}_2\right)\right]}^n$ is power function and $x\left(s_1\mathrm{,}{s}_2\right)$ is unknown, n is a positive integer. $\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ is determinate kernel function $0\le s_1\le t_1\le T_1\mathrm{,0}\le s_2\le t_2\le T_2$ . The authors revealed the accuracy and efficiency of the proposed method by some examples and gave the rate of convergence to the numerical solution.

However, as far as we known, there are hardly any papers about the numerical solution of two-dimensional nonlinear stochastic Itô-Volterra integral equations. Inspired by the above literatures, we introduce an efficient numerical method for the following nonlinear stochastic integral equation based on block pulse functions.

$\begin{array}{c}x\left(t_1,t_2\right)=x_0\left(t_1,t_2\right)+{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\tilde{k}\left(t_1,t_2,s_1,s_2\right)\sigma \left(x\left(s_1,s_2\right)\right)\text{d}{s}_1\text{d}{s}_2\\ +{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\hat{k}\left(t_1,t_2,s_1,s_2\right)g\left(x\left(s_1,s_2\right)\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right),\end{array}$ (3)

where $x\left(t_1\mathrm{,}{t}_2\right)$ is unknown function and is called the solution of the Equation (3) defined on district $D=\left[0,1\right)\times \left[0,1\right)$ . $x_0\left(t_1\mathrm{,}{t}_2\right)$ is known determinate function. $\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ and $\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ are determinate kernel functions. ${\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\hat{k}\left(t_1,t_2,s_1,s_2\right)g\left(x\left(s_1,s_2\right)\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)$ is the double Itô integral. $B\left(s_1\right)$ and $B\left(s_2\right)$ are two independent Brownian motions. $\sigma$ and g are analytical functions.

In Section 2, we recall the definition and properties of block pulse function. In Section 3 and 4, we show the integration operational matrix about two-dimensional block pulse functions. In Section 5, an efficient numerical method to nonlinear stochastic Itô-Volterra integral equation is obtained. In Section 6, the error and the rate of convergence of this method are given. It’s important to emphasize that the error is analyzed by Gronwall’s inequality and the interchangeability of integral and expectation. However, the norm was used in the literature [11] , it is a pity that the interchangeability of norm and integral wasn’t proved. In Section 7, we give a numerical example to illustrate the validity of the method. In the final Section 8, we make some conclusions and look ahead to further work.

2. Two-Dimensional Block Pulse Functions

One dimensional block pulse functions (BPFs) have been widely studied and applied to solve different problems. For example, the article [12] and their relative references give a detailed description. A $m_1{m}_2$ -set of two-dimensional block pulse functions (2D-BPFs) ${\varphi }_{a_1\mathrm{,}{a}_2}\left(t_1\mathrm{,}{t}_2\right)$ in the region of $D=\left[0,1\right)\times \left[0,1\right)$ are defined as:

${\varphi }_{a_1\mathrm{,}{a}_2}\left(t_1\mathrm{,}{t}_2\right)=\{\begin{array}{l}1\left(a_1-1\right)h_1\le t_1<a_1{h}_1,\left(a_2-1\right)h_2\le t_2<a_2{h}_2\\ 0\text{otherwise}\mathrm{,}\end{array}$

where $a_i=1,2,\cdots ,m_i$ , $h_i=\frac{1}{m_i}$ , $m_i=2^n$ , $m_i$ and n are arbitrary positive integers and $i=1,2$ .

Similar to the one-dimensional case [12] . There are some elementary properties for 2D-BPFs as follows:

1) Disjointness:

${\varphi }_{a_1,a_2}\left(t_1,t_2\right){\varphi }_{b_1,b_2}\left(t_1,t_2\right)=\{\begin{array}{l}{\varphi }_{a_1,a_2}\left(t_1,t_2\right)\text{if}{a}_1=b_1,a_2=b_2\\ 0\text{otherwise},\end{array}$ (4)

where $a_i,b_i=1,2,\cdots ,m_i$ , $i=1,2$ .

2) Orthogonality:

${\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}{\varphi }_{a_1,a_2}\left(t_1,t_2\right){\varphi }_{b_1,b_2}\left(t_1,t_2\right)\text{d}{t}_1\text{d}{t}_2=\{\begin{array}{l}{h}_1{h}_2\text{if}{a}_1=b_1,a_2=b_2\\ 0\text{otherwise}.\end{array}$ (5)

3) Completeness: for every $f\in \left(L^2\left(D\right)\right)$ , when $m_1$ and $m_2$ approach to the infinity, Parseval’s identity holds:

${\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}{f}^2\left(t_1,t_2\right)\text{d}{t}_1\text{d}{t}_2={\displaystyle \underset{a_1=1}{\overset{\infty }{\sum }}}{\displaystyle \underset{a_2=1}{\overset{\infty }{\sum }}}{f}_{a_1,a_2}^2{\Vert {\varphi }_{a_1,a_2}\left(t_1,t_2\right)\Vert }^2,$ (6)

where

$f_{a_1,a_2}=\frac{1}{h_1{h}_2}{\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}f\left(t_1,t_2\right){\varphi }_{a_1,a_2}\left(t_1,t_2\right)\text{d}{t}_1\text{d}{t}_2.$

The set of 2D-BPFs may be written as a vector $\Phi \left(t_1\mathrm{,}{t}_2\right)$ of dimension $\left(m_1{m}_2\right)$ :

${\Phi }_{m_1{m}_2}\left(t_1,t_2\right)={\left({\varphi }_{1,1}\left(t_1,t_2\right),\cdots ,{\varphi }_{1,m_2}\left(t_1,t_2\right),\cdots ,{\varphi }_{m_1,1}\left(t_1,t_2\right),\cdots ,{\varphi }_{m_1,m_2}\left(t_1,t_2\right)\right)}^{\text{T}},$ (7)

where $\left(t_1\mathrm{,}{t}_2\right)\in D$ .

From the above representation and disjointness property, it follows that:

${\Phi }_{m_1{m}_2}\left(t_1,t_2\right){\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1,t_2\right)={\left(\begin{array}{cccc}\text{}\text{}\text{}{\varphi }_{1,1}\left(t_1,t_2\right)\text{}\text{}\text{}& \text{}\text{}\text{}0\text{}\text{}\text{}& \text{}\text{}\text{}\cdots \text{}\text{}\text{}& \text{}\text{}\text{}0\text{}\text{}\text{}\\ \text{}\text{}\text{}0\text{}\text{}\text{}& \text{}\text{}\text{}{\varphi }_{1,2}\left(t_1,t_2\right)\text{}\text{}\text{}& \text{}\text{}\text{}\cdots \text{}\text{}\text{}& \text{}\text{}\text{}0\text{}\text{}\text{}\\ \text{}\text{}\text{}⋮\text{}\text{}\text{}& \text{}\text{}\text{}⋮\text{}\text{}\text{}& \text{}\text{}\text{}\ddots \text{}\text{}\text{}& \text{}\text{}\text{}⋮\text{}\text{}\text{}\\ \text{}\text{}\text{}0\text{}\text{}\text{}& \text{}\text{}\text{}0\text{}\text{}\text{}& \text{}\text{}\text{}\cdots \text{}\text{}\text{}& \text{}\text{}\text{}{\varphi }_{m_1,m_2}\left(t_1,t_2\right)\text{}\text{}\text{}\end{array}\right)}_{m_1{m}_2\times m_1{m}_2},$ (8)

${\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1,t_2\right){\Phi }_{m_1{m}_2}\left(t_1,t_2\right)=1,$

${\Phi }_{m_1{m}_2}\left(t_1,t_2\right){\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1,t_2\right)G=\tilde{G}{\Phi }_{m_1{m}_2}\left(t_1,t_2\right),$ (9)

where G is a $\left(m_1{m}_2\right)$ -vector and the matrix $\tilde{G}=diag\left(G\right)$ . Moreover, it is easy to conclude that for every $\left(m_1{m}_2\right)\times \left(m_1{m}_2\right)$ matrix $A$

${\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)A{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)={\hat{A}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\mathrm{,}$ (10)

where $\hat{A}$ is a $\left(m_1{m}_2\right)$ -vector with elements equal to the diagonal entries of matrix $A$ .

Any function $x\left(t_1\mathrm{,}{t}_2\right)$ which is square integrable in the interval D can be expanded in terms of BPFs as

$x\left(t_1\mathrm{,}{t}_2\right)\simeq x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)={\displaystyle \underset{a_1=1}{\overset{m_1}{\sum }}}{\displaystyle \underset{a_2=1}{\overset{m_2}{\sum }}}{x}_{a_1\mathrm{,}{a}_2}{\varphi }_{a_1\mathrm{,}{a}_2}\left(t_1\mathrm{,}{t}_2\right)=X_{m_1{m}_2}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t\right)\mathrm{,}$ (11)

where $x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)$ is $m_1{m}_2$ approximations of 2D-BPFs of $x\left(t_1\mathrm{,}{t}_2\right)$ , $x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)$ is a coefficient $\left(m_1{m}_2\right)$ -vector, i.e.

$X_{m_1{m}_2}={\left(x_{1,1},\cdots ,x_{1,m_2},\cdots ,x_{m_1,1},\cdots ,x_{m_1,m_2}\right)}^{\text{T}},$ (12)

where the block pulse coefficients $x_{a_1\mathrm{,}{a}_2}$ are obtained as

$x_{a_1,a_2}=\frac{1}{h_1{h}_2}{\displaystyle {\int }_{\left(a_2-1\right)h_2}^{a_2{h}_2}}{\displaystyle {\int }_{\left(a_1-1\right)h_1}^{a_1{h}_1}}x\left(t_1,t_2\right)\text{d}{t}_1\text{d}{t}_2.$

Similarly, a function of four variables $k\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ on $L^2\left(D\times D\right)$ may be approximated with respect to 2D-BPFs such as

$k\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\simeq {\Phi }_{m_1{m}_2}{\left(t_1\mathrm{,}{t}_2\right)}^{\text{T}}K{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\mathrm{,}$

where ${\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)$ is a 2D-BPFs vector of dimension $\left(m_1{m}_2\right)$ , $K$ is the $\left(m_1{m}_2\right)\times \left(m_1{m}_2\right)$ two-dimensional block pulse coefficient matrix in the following form

$K={\left(K_{a_1{b}_1}\right)}_{m_1\times m_1}\mathrm{,}{K}_{a_1{b}_1}={\left(k_{a_1{a}_2{b}_1{b}_2}\right)}_{m_2\times m_2}\mathrm{,}$

$a_i,b_i=1,\cdots ,m_i,i=1,2$ and two-dimensional block pulse coefficients $k_{a_1{a}_2{b}_1{b}_2}$ are given by

$k_{a_1{a}_2{b}_1{b}_2}=\frac{1}{h_1^2{h}_2^2}\left[{\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}k\left(t_1,t_2,s_1,s_2\right){\varphi }_{b_1,b_2}\left(t_1,t_2\right){\varphi }_{a_1,a_2}\left(s_1,s_2\right)\text{d}{t}_1\text{d}{t}_2\text{d}{s}_1\text{d}{s}_2\right].$ (13)

The more details can also reference to [7] .

3. Operational Matrix of Integration

Let $M={\left({\xi }_{ij}\right)}_{M_1\times M_2}$ and $N={\left({\eta }_{ij}\right)}_{N_1\times N_2}$ be matrices. $M_l\mathrm{,}{N}_l$ are positive

integers, $l=1,2$ . We have

$M\otimes N=\left({\xi }_{ij}N\right)={\left(\begin{array}{cccc}{\xi }_{11}N& {\xi }_{12}N& \cdots & {\xi }_{1M_2}N\\ {\xi }_{21}N& {\xi }_{22}N& \cdots & {\xi }_{2M_2}N\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_{M_{1}1}N& {\xi }_{M_{1}2}N& \cdots & {\xi }_{M_1{M}_2}N\end{array}\right)}_{M_1{N}_1\times M_2{N}_2}\mathrm{,}$

where $\otimes$ denotes the Kronecker product defined as [13] . Each ${\xi }_{ij}N$ is a block of size $N_1\times N_2$ , $M\otimes N$ is of size $M_1{N}_1\times M_2{N}_2$ .

Then the vector ${\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)$ can be showed as following

$\begin{array}{l}{\Phi }_{m_1{m}_2}\left(t_1,t_2\right)\\ ={\Phi }_{m_1}\left(t_1\right)\otimes {\Phi }_{m_2}\left(t_2\right)\\ ={\left({\varphi }_1\left(t_1\right),{\varphi }_2\left(t_1\right),\cdots ,{\varphi }_{m_1}\left(t_1\right)\right)}^{\text{T}}\otimes {\left({\varphi }_1\left(t_2\right),{\varphi }_2\left(t_2\right),\cdots ,{\varphi }_{m_2}\left(t_2\right)\right)}^{\text{T}}\\ ={\left({\varphi }_1\left(t_1\right){\varphi }_1\left(t_2\right),\cdots ,{\varphi }_1\left(t_1\right){\varphi }_{m_2}\left(t_2\right),\cdots ,{\varphi }_{m_1}\left(t_1\right){\varphi }_1\left(t_2\right),\cdots ,{\varphi }_{m_1}\left(t_1\right){\varphi }_{m_2}\left(t_2\right)\right)}^{\text{T}}\end{array}$

where ${\varphi }_{a_i}\left(t_i\right)$ are one dimensional BPFs, ${\Phi }_{m_i}\left(t_i\right)$ are vectors of one dimensional BPFs, $a_i=1,2,\cdots ,m_i,i=1,2$ .

The integration of the vector ${\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)$ defined in (7) can be approximately obtained as following

$\begin{array}{c}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\text{d}{s}_1\text{d}{s}_2={\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1}\left(s_1\right)\otimes {\Phi }_{m_2}\left(s_2\right)\text{d}{s}_1\text{d}{s}_2\\ ={\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1}\left(s_1\right)\text{d}{s}_1\otimes {\displaystyle {\int }_0^{t_2}}{\Phi }_{m_2}\left(s_2\right)\text{d}{s}_2\\ \simeq P_1{\Phi }_{m_1}\left(t_1\right)\otimes P_2{\Phi }_{m_2}\left(t_2\right)\\ ={\left(P_1\otimes P\right)}_2{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\\ =P{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\mathrm{,}\end{array}$ (14)

where $t_1\in \left[\mathrm{0,1}\right)\mathrm{,}{t}_2\in \left[\mathrm{0,1}\right)$ , $P$ is the $\left(m_1{m}_2\right)\times \left(m_1{m}_2\right)$ operational matrix of integration for 2D-BPFs and $P_i$ , $\left(i=1,2\right)$ are the operational matrix of one-dimensional BPFs [12] defined over $\left[\mathrm{0,1}\right)$ as following.

$P_i=\frac{h}{2}{\left(\begin{array}{ccccc}1& 2& 2& \cdots & 2\\ 0& 1& 2& \cdots & 2\\ 0& 0& 1& \cdots & 2\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & 1\end{array}\right)}_{\left(m_i\times m_i\right)}\mathrm{.}$

For details, see [7] , so

$\begin{array}{c}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}x\left(s_1\mathrm{,}{s}_2\right)\text{d}{s}_1\text{d}{s}_2\simeq {\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{X}_{m_1{m}_2}^T{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\text{d}{s}_1\text{d}{s}_2\\ =X_{m_1{m}_2}^{T}P{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\mathrm{.}\end{array}$ (15)

4. Stochastic Integration Operational Matrix

Similarly, we obtain the stochastic integration of the vector ${\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)$ defined in (7) as following

$\begin{array}{l}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)\\ ={\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1}\left(s_1\right)\otimes {\Phi }_{m_2}\left(s_2\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)\\ ={\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1}\left(s_1\right)\text{d}B\left(s_1\right)\otimes {\displaystyle {\int }_0^{t_2}}{\Phi }_{m_2}\left(s_2\right)\text{d}B\left(s_2\right)\\ \simeq P_{s_1}{\Phi }_{m_1}\left(t_1\right)\otimes P_{s_2}{\Phi }_{m_2}\left(t_2\right)\\ =\left(P_{s_1}\otimes P_{s_2}\right){\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)=P_s{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\mathrm{,}\end{array}$ (16)

where $t_1\in \left[\mathrm{0,1}\right)\mathrm{,}{t}_2\in \left[\mathrm{0,1}\right)$ , $P_s$ is the $\left(m_1{m}_2\right)\times \left(m_1{m}_2\right)$ stochastic operational matrix of integration for 2D-BPFs and $P_{s_i}$ , $\left(i=1,2\right)$ are the stochastic operational matrix of one-dimensional BPFs [12] defined over $\left[\mathrm{0,1}\right)$ as following.

$P_{s_i}={\left(\begin{array}{ccccc}B\left(\frac{h_i}{2}\right)\text{}\text{}\text{}& B\left(h_i\right)\text{}\text{}& B\left(h_i\right)\text{}\text{}& \cdots & B\left(h_i\right)\text{}\text{}\\ 0& B\left(\frac{3h_i}{2}\right)-B\left(h_i\right)\text{}\text{}& B\left(2h_i\right)-B\left(h_i\right)\text{}\text{}& \cdots & B\left(2h_i\right)-B\left(h_i\right)\text{}\text{}\\ 0& 0& B\left(\frac{5h_i}{2}\right)-B\left(5h_i\right)\text{}\text{}& \cdots & B\left(3h_i\right)-B\left(2h_i\right)\text{}\text{}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \cdots & B\left(\frac{\left(2m_i-1\right)h_i}{2}\right)-B\left(\left(m_i-1\right)h_i\right)\end{array}\right)}_{m_i\times m_i}\mathrm{.}$ (17)

For details, see [7] . Therefore,

$\begin{array}{c}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}x\left(s_1\mathrm{,}{s}_2\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)\simeq {\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{X}_{m_1{m}_2}^{\text{T}}{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)\\ =X_{m_1{m}_2}^{\text{T}}{P}_s{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\mathrm{.}\end{array}$ (18)

5. Numerical Method

In this section, we first provide a useful result for solving two-dimensional nonlinear stochastic Itô-Volterra integral Equation (3).

Lemma 1. Let $\sigma \left(t\right)={\displaystyle \sum }{a}_j{t}^j$ , $g\left(t\right)={\displaystyle \sum }{b}_j{t}^j$ be the analytic functions for positive integer $j\in \left(\mathrm{0,}\infty \right)$ , then

$\sigma \left(x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\right)={\sigma }^{\text{T}}\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right),$

$g\left(x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\right)=g^{\text{T}}\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right),$

where ${\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)$ and $X_{m_1{m}_2}$ are derived in (7) and (12),

${\sigma }^{\text{T}}\left(X_{m_1{m}_2}\right)=\left(\sigma \left(x_{1,1}\right),\cdots ,\sigma \left(x_{1,m_2}\right),\cdots ,\sigma \left(x_{m_1,1}\right),\cdots ,\sigma \left(x_{m_1,m_2}\right)\right),$

$g^{\text{T}}\left(X_{m_1{m}_2}\right)=\left(g\left(x_{1,1}\right),\cdots ,g\left(x_{1,m_2}\right),\cdots ,g\left(x_{m_1,1}\right),\cdots ,g\left(x_{m_1,m_2}\right)\right).$

Proof. By virtue of the known conditions and the disjointness properties of 2D-BPFs defined in (4), we can get

$\begin{array}{l}\sigma \left(x_{m_1{m}_2}\left(t_1,t_2\right)\right)={\displaystyle \sum }{a}_j{\left(x_{m_1{m}_2}\left(t_1,t_2\right)\right)}^j={\displaystyle \sum }{a}_j{\left[{\displaystyle \underset{a_1=1}{\overset{m_1}{\sum }}}{\displaystyle \underset{a_2=1}{\overset{m_2}{\sum }}}{x}_{a_1,a_2}{\varphi }_{a_1,a_2}\left(t_1,t_2\right)\right]}^j\\ ={\displaystyle \sum }{a}_j[x_{1,1}{\varphi }_{1,1}\left(t_1,t_2\right)+\cdots +x_{1,m_2}{\varphi }_{1,m_2}\left(t_1,t_2\right)+\cdots +x_{m_1,1}{\varphi }_{m_1,1}\left(t_1,t_2\right)\\ +{\cdots +x_{m_1,m_2}{\varphi }_{m_1,m_2}\left(t_1,t_2\right)]}^j\\ ={\displaystyle \sum }{a}_j\left(x_{1,1}^j,\cdots ,x_{1,m_2}^j,\cdots ,x_{m_1,1}^j,\cdots ,x_{m_1,m_2}^j\right){\Phi }_{m_1{m}_2}\left(t_1,t_2\right)\\ ={\sigma }^{\text{T}}\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(t_1,t_2\right),\end{array}$

thus,

$\sigma \left(x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\right)={\sigma }^{\text{T}}\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)={\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)\sigma \left(X_{m_1{m}_2}\right)\mathrm{,}$ (19)

$g\left(x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\right)=g^{\text{T}}\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)={\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)g\left(X_{m_1{m}_2}\right)\mathrm{.}$ (20)

The proof is completed. □

Now we suppose $x\left(t_1\mathrm{,}{t}_2\right)$ , $x_0\left(t_1\mathrm{,}{t}_2\right)$ , $\sigma \left(x\left(t_1\mathrm{,}{t}_2\right)\right)$ , $g\left(x\left(t_1\mathrm{,}{t}_2\right)\right)$ , $\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ and $\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ can be approximated in terms of 2D-BPFs.

$x\left(t_1\mathrm{,}{t}_2\right)\simeq x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)=X_{m_1{m}_2}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)={\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)X_{m_1{m}_2}\mathrm{,}$ (21)

$x_0\left(t_1\mathrm{,}{t}_2\right)\simeq x_{0_{m_1{m}_2}}\left(t_1\mathrm{,}{t}_2\right)=X_{0_{m_1{m}_2}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)={\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)X_{0_{m_1{m}_2}}\mathrm{,}$ (22)

$\sigma \left(x\left(t_1\mathrm{,}{t}_2\right)\right)\simeq \sigma \left(x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\right)={\sigma }^{\text{T}}\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)={\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)\sigma \left(X_{m_1{m}_2}\right)\mathrm{,}$ (23)

$g\left(x\left(t_1\mathrm{,}{t}_2\right)\right)\simeq g\left(x_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\right)=g^{\text{T}}\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)={\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)g\left(X_{m_1{m}_2}\right)\mathrm{,}$ (24)

$\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\simeq {\tilde{k}}_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)={\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_1{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\mathrm{,}$ (25)

$\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\simeq {\hat{k}}_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)={\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_2{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\mathrm{,}$ (26)

where $X_{m_1{m}_2}$ , $X_{0_{m_1{m}_2}}$ , $\sigma \left(X_{m_1{m}_2}\right)$ and $g\left(X_{m_1{m}_2}\right)$ are two-dimensional block

pulse coefficient vectors. $K_1$ and $K_2$ are two-dimensional block pulse coefficient matrices.

Now, by (21)-(26), we approximate the Equation (3)

$\begin{array}{l}{X}_{m_1{m}_2}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\\ =X_{0_{m_1{m}_2}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\\ +{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_1{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right){\Phi }_{m_1{m}_2}^{\text{T}}\left(s_1\mathrm{,}{s}_2\right)\sigma \left(X_{m_1{m}_2}\right)\text{d}{s}_1\text{d}{s}_2\\ +{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_2{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right){\Phi }_{m_1{m}_2}^{\text{T}}\left(s_1\mathrm{,}{s}_2\right)g\left(X_{m_1{m}_2}\right)\text{d}B\left(s_1\right)\text{d}B(s2)\end{array}$

$\begin{array}{l}=X_{0_{m_1{m}_2}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\\ +{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_1{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right){\Phi }_{m_1{m}_2}^{\text{T}}\left(s_1\mathrm{,}{s}_2\right)\sigma \left(X_{m_1{m}_2}\right)\text{d}{s}_1\text{d}{s}_2\\ +{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right){\Phi }_{m_1{m}_2}^{\text{T}}\left(s_1\mathrm{,}{s}_2\right)g\left(X_{m_1{m}_2}\right)\text{d}B\left(s_1\right)\text{d}B(s2)\end{array}$

$\begin{array}{l}=X_{0_{m_1{m}_2}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\\ +{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_1{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\tilde{\sigma }\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\text{d}{s}_1\text{d}{s}_2\\ +{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\tilde{g}\left(X_{m_1{m}_2}\right){\Phi }_{m_1{m}_2}\left(s_1\mathrm{,}{s}_2\right)\text{d}B\left(s_1\right)\text{d}B(s2)\end{array}$

$\begin{array}{l}=X_{0_{m_1{m}_2}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\\ +{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_1\tilde{\sigma }\left(X_{m_1{m}_2}\right){\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\text{d}{s}_1\text{d}{s}_2\\ +{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_2\tilde{g}\left(X_{m_1{m}_2}\right){\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)\mathrm{,}\end{array}$

by (15) and (18), we have

$\begin{array}{l}{X}_{m_1{m}_2}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\\ =X_{0_{m_1{m}_2}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)+{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_1\tilde{\sigma }\left(X_{m_1{m}_2}\right)P{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\\ \text{}\text{}+{\Phi }_{m_1{m}_2}^{\text{T}}\left(t_1\mathrm{,}{t}_2\right)K_2\tilde{g}\left(X_{m_1{m}_2}\right)P_s{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\mathrm{,}\end{array}$

let $Q=K_1\tilde{\sigma }\left(X_{m_1{m}_2}\right)P$ and $Q_s=K_2\tilde{g}\left(X_{m_1{m}_2}\right)P_s$ , they both are $\left(m_1{m}_2\right)\times \left(m_1{m}_2\right)$ matrices. By (10), we have

$X_{m_1{m}_2}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)=X_{0_{m_1{m}_2}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)+{\hat{Q}}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)+{\hat{Q}}_s^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1\mathrm{,}{t}_2\right)\mathrm{,}$

where $\hat{Q}$ and ${\hat{Q}}_s$ are $\left(m_1{m}_2\right)$ -vectors with elements equal to the diagonal entries of matrices $Q$ and $Q_s$ . Then

$X_{m_1{m}_2}^{\text{T}}=X_{0_{m_1{m}_2}}^{\text{T}}+{\hat{Q}}^{\text{T}}+{\hat{Q}}_s^{\text{T}}.$ (27)

There are various methods to solve the nonlinear system of Equation (27) of $X_{m_1{m}_2}$ . In this paper, we will use the int () function provided by Matlab 2015b [14] to solve it. According to the coefficient vector $X_{m_1{m}_2}$ , we obtain that the approximation solution of Equation (3) $x_{m_1{m}_2}\left(t_1,t_2\right)=X_{m_1{m}_2}^{\text{T}}{\Phi }_{m_1{m}_2}\left(t_1,t_2\right)$ .

6. Error Analysis

In this section, for convenience, we assume $m_1=m_2=m$ and prove that the approximation solution is convergent of order $O\left(h\right),h=\frac{1}{m}$ .

Lemma 2. Let $v\left(s_1\mathrm{,}{s}_2\right)$ be an arbitrary bounded function on $D=\left[0,1\right)\times \left[0,1\right)$ and ${\tilde{e}}_{mm}\left(s_1,s_2\right)=v\left(s_1,s_2\right)-v_{mm}\left(s_1,s_2\right)$ , which $v_{mm}\left(s_1\mathrm{,}{s}_2\right)$ is m² approximations of 2D-BPFs of $v\left(s_1\mathrm{,}{s}_2\right)$ , then

${\Vert \tilde{e}\Vert }_{L^2\left(D\right)}^2={\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}{\tilde{e}}_{mm}^2\left(s_1\mathrm{,}{s}_2\right)\text{d}{s}_1\text{d}{s}_2\le O\left(h^2\right)\mathrm{.}$ (28)

Proof. Similar to [15] [16] . □

Lemma 3. Let $v\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ be an arbitrary bounded function on $D\times D$ and ${\hat{e}}_{mm}\left(t_1,t_2,s_1,s_2\right)=v\left(t_1,t_2,s_1,s_2\right)-v_{mm}\left(t_1,t_2,s_1,s_2\right)$ , which $v_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ is m² approximations of 2D-BPFs of $v\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ , then

${\Vert \hat{e}\Vert }_{L^2\left(D\times D\right)}^2={\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}{\displaystyle {\int }_0^1}{\hat{e}}_{mm}^2\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\text{d}{s}_1\text{d}{s}_2\text{d}{t}_1\text{d}{t}_2\le O\left(h^2\right)\mathrm{.}$ (29)

Proof. Similar to [15] [16] . □

Next, let

$\begin{array}{l}{e}_{mm}\left(t_1\mathrm{,}{t}_2\right)=x\left(t_1\mathrm{,}{t}_2\right)-x_{mm}\left(t_1\mathrm{,}{t}_2\right)\\ =x_0\left(t_1\mathrm{,}{t}_2\right)-x_{0_{mm}}\left(t_1\mathrm{,}{t}_2\right)+{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}[\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\sigma \left(x\left(s_1\mathrm{,}{s}_2\right)\right)\\ -{\tilde{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\sigma \left(x_{mm}\left(s_1\mathrm{,}{s}_2\right)\right)]\text{d}{s}_1\text{d}{s}_2\\ +{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}[\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)g\left(x\left(s_1\mathrm{,}{s}_2\right)\right)\\ -{\hat{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)g\left(x_{mm}\left(s_1\mathrm{,}{s}_2\right)\right)]\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)\mathrm{.}\end{array}$ (30)

where $x_{mm}\left(t_1\mathrm{,}{t}_2\right)$ is the approximation solution of $x\left(t_1\mathrm{,}{t}_2\right)$ defined in (3), $x_{0_{mm}}\left(t_1\mathrm{,}{t}_2\right)$ , ${\tilde{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ and ${\hat{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ are m² approximations of 2D-BPFs of $x_0\left(t_1\mathrm{,}{t}_2\right)\mathrm{,}\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ and $\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)$ , respectively.

Theorem 1. For analytic functions $\sigma$ and g, there are constant numbers satisfy the following conditions:

1) $\left|\sigma \left(x\right)-\sigma \left(y\right)\right|\le l_1\left|x-y\right|$ , $\left|g\left(x\right)-g\left(y\right)\right|\le l_3\left|x-y\right|$ ,

2) $\left|\sigma \left(x\right)\right|\le l_2\mathrm{,}\left|g\left(y\right)\right|\le l_4$ ,

where $x\mathrm{,}y\in R$ and let $\left|\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right|\le l_5\mathrm{,}\left|\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right|\le l_6$ be determinate bounded kernel functions, where $l_i$ , $i=1,2,\cdots ,6$ are constant numbers. Then,

$\begin{array}{l}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}\mathbb{E}\left({\left|e_{mm}\left(t_1\mathrm{,}{t}_2\right)\right|}^2\right)\text{d}{t}_1\text{d}{t}_2\\ ={\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}\mathbb{E}\left({\left|x\left(t_1\mathrm{,}{t}_2\right)-x_{mm}\left(t_1\mathrm{,}{t}_2\right)\right|}^2\right)\text{d}{t}_1\text{d}{t}_2\le O\left(h^2\right)\mathrm{,}T\in \left[\mathrm{0,1}\right)\mathrm{.}\end{array}$

Proof. For (30), we have

$\begin{array}{c}\mathbb{E}\left({\left|e_{mm}\left(t_1,t_2\right)\right|}^2\right)\le 3[\mathbb{E}\left({\left|x_0\left(t_1,t_2\right)-x_{0_{mm}}\left(t_1,t_2\right)\right|}^2\right)\\ +\mathbb{E}(|{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}[\tilde{k}\left(t_1,t_2,s_1,s_2\right)\sigma \left(x\left(s_1,s_2\right)\right)\\ -{{\tilde{k}}_{mm}\left(t_1,t_2,s_1,s_2\right)\sigma \left(x_{mm}\left(s_1,s_2\right)\right)]\text{d}{s}_1\text{d}{s}_2|}^2)\\ +\mathbb{E}(|{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}[\hat{k}\left(t_1,t_2,s_1,s_2\right)g\left(x\left(s_1,s_2\right)\right)\\ -{{\hat{k}}_{mm}\left(t_1,t_2,s_1,s_2\right)g\left(x_{mm}\left(s_1,s_2\right)\right)]\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)|}^2)].\end{array}$

According to Itô isometry, Cauchy-Schwartz inequality and Lipschitz conditions, we can write

$\begin{array}{l}\mathbb{E}\left({\left|e_{mm}\left(t_1,t_2\right)\right|}^2\right)\\ \le 3[\mathbb{E}\left({\left|x_0\left(t_1,t_2\right)-x_{0_{mm}}\left(t_1,t_2\right)\right|}^2\right)\\ +\mathbb{E}\left({\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\sigma \left(x\left(s_1\mathrm{,}{s}_2\right)\right)-{\tilde{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\sigma \left(x_{mm}\left(s_1\mathrm{,}{s}_2\right)\right)\right|}^2\text{d}{s}_1\text{d}{s}_2\right)\\ +\mathbb{E}\left({\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)g\left(x\left(s_1\mathrm{,}{s}_2\right)\right)-{\hat{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)g\left(x_{mm}\left(s_1\mathrm{,}{s}_2\right)\right)\right|}^2\text{d}{s}_1\text{d}{s}_2\right)]\end{array}$

$\begin{array}{l}=3[\mathbb{E}\left({\left|x_0\left(t_1,t_2\right)-x_{0_{mm}}\left(t_1,t_2\right)\right|}^2\right)\\ +{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\mathbb{E}(|\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\left[\sigma \left(x\left(s_1\mathrm{,}{s}_2\right)\right)-\sigma \left(x_{mm}\left(s_1\mathrm{,}{s}_2\right)\right)\right]\\ +{\sigma \left(x_{mm}\left(s_1\mathrm{,}{s}_2\right)\right)\left[\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)-{\tilde{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right]|}^2)\text{d}{s}_1\text{d}{s}_2\\ +{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\mathbb{E}(|\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\left[g\left(x\left(s_1\mathrm{,}{s}_2\right)\right)-g\left(x_{mm}\left(s_1\mathrm{,}{s}_2\right)\right)\right]\\ +{g\left(x_{mm}\left(s_1\mathrm{,}{s}_2\right)\right)\left[\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)-{\hat{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right]|}^2)\text{d}{s}_1\text{d}{s}_2]\end{array}$

$\begin{array}{l}\le 3[{\left|x_0\left(t_1,t_2\right)-x_{0_{mm}}\left(t_1,t_2\right)\right|}^2\\ +2l_1^2{l}_5^2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\mathbb{E}\left({\left|e_{mm}\left(s_1,s_2\right)\right|}^2\right)\text{d}{s}_1\text{d}{s}_2\\ +2l_2^2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|\tilde{k}\left(t_1,t_2,s_1,s_2\right)-{\tilde{k}}_{mm}\left(t_1,t_2,s_1,s_2\right)\right|}^2\text{d}{s}_1\text{d}{s}_2\\ +2l_3^2{l}_6^2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\mathbb{E}\left({\left|e_{mm}\left(s_1,s_2\right)\right|}^2\right)\text{d}{s}_1\text{d}{s}_2\\ +2l_4^2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|\hat{k}\left(t_1,t_2,s_1,s_2\right)-{\hat{k}}_{mm}\left(t_1,t_2,s_1,s_2\right)\right|}^2\text{d}{s}_1\text{d}{s}_2].\end{array}$

Then, we can get

$\mathbb{E}\left({\left|e_{mm}\left(t_1,t_2\right)\right|}^2\right)\le \beta \left(t_1\mathrm{,}{t}_2\right)+\alpha {\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\mathbb{E}\left({\left|e_{mm}\left(s_1,s_2\right)\right|}^2\right)\text{d}{s}_1\text{d}{s}_2\mathrm{,}$

where,

$\alpha =6\left(l_1^2{l}_5^2+l_3^2{l}_6^2\right)\mathrm{.}$

$\begin{array}{c}\beta \left(t_1\mathrm{,}{t}_2\right)=3[{\left|x_0\left(t_1\mathrm{,}{t}_2\right)-x_{0_{mm}}\left(t_1\mathrm{,}{t}_2\right)\right|}^2\\ +2l_2^2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)-{\tilde{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right|}^2\text{d}{s}_1\text{d}{s}_2\\ +2l_4^2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)-{\hat{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right|}^2\text{d}{s}_1\text{d}{s}_2].\end{array}$

Let $f\left(t_1\mathrm{,}{t}_2\right)=\mathbb{E}\left({\left|e_{mm}\left(t_1,t_2\right)\right|}^2\right)$ , we get

$f\left(t_1\mathrm{,}{t}_2\right)\le \beta \left(t_1\mathrm{,}{t}_2\right)+\alpha {\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}f\left({\tau }_1\mathrm{,}{\tau }_2\right)\text{d}{\tau }_1\text{d}{\tau }_2\mathrm{,}{\tau }_1\in \left[\mathrm{0,}{t}_1\right)\mathrm{,}{\tau }_2\in \left[\mathrm{0,}{t}_2\right)\mathrm{.}$

By Gronwall’s inequality, we have

$f\left(t_1\mathrm{,}{t}_2\right)\le \beta \left(t_1\mathrm{,}{t}_2\right)+\alpha {\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\text{e}}^{{\displaystyle {\int }_0^{{\tau }_2}}{\displaystyle {\int }_0^{{\tau }_1}}\alpha \text{d}{s}_1\text{d}{s}_2}\beta \left({\tau }_1\mathrm{,}{\tau }_2\right)\text{d}{\tau }_1\text{d}{\tau }_2\mathrm{,}{t}_1\mathrm{,}{t}_2\in \left[\mathrm{0,1}\right)\mathrm{.}$

Then, for $T\in \left[\mathrm{0,1}\right)$

$\begin{array}{l}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}f\left(t_1\mathrm{,}{t}_2\right)\text{d}{t}_1\text{d}{t}_2\\ ={\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}\mathbb{E}\left({\left|e_{mm}\left(t_1,t_2\right)\right|}^2\right)\text{d}{t}_1\text{d}{t}_2\\ \le {\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}\left(\beta \left(t_1\mathrm{,}{t}_2\right)+\alpha {\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\text{e}}^{{\displaystyle {\int }_0^{{\tau }_2}}{\displaystyle {\int }_0^{{\tau }_1}}\alpha \text{d}{s}_1\text{d}{s}_2}\beta \left({\tau }_1\mathrm{,}{\tau }_2\right)\text{d}{\tau }_1\text{d}{\tau }_2\right)\text{d}{t}_1\text{d}{t}_2\\ ={\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}\beta \left(t_1\mathrm{,}{t}_2\right)\text{d}{t}_1\text{d}{t}_2+\alpha {\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\text{e}}^{{\displaystyle {\int }_0^{{\tau }_2}}{\displaystyle {\int }_0^{{\tau }_1}}\alpha \text{d}{s}_1\text{d}{s}_2}\beta \left({\tau }_1\mathrm{,}{\tau }_2\right)\text{d}{\tau }_1\text{d}{\tau }_2\text{d}{t}_1\text{d}{t}_2\\ \le {\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}\beta \left(t_1\mathrm{,}{t}_2\right)\text{d}{t}_1\text{d}{t}_2+\alpha {\text{e}}^{\alpha T^2}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\beta \left({\tau }_1\mathrm{,}{\tau }_2\right)\text{d}{\tau }_1\text{d}{\tau }_2\text{d}{t}_1\text{d}{t}_2\end{array}$

$\begin{array}{l}=3{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}{\left|x_0\left(t_1\mathrm{,}{t}_2\right)-x_{0_{mm}}\left(t_1\mathrm{,}{t}_2\right)\right|}^2\text{d}{t}_1\text{d}{t}_2\\ +6l_2^2{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|\tilde{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)-{\tilde{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right|}^2\text{d}{s}_1\text{d}{s}_2\text{d}{t}_1\text{d}{t}_2\\ +6l_4^2{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|\hat{k}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)-{\hat{k}}_{mm}\left(t_1\mathrm{,}{t}_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right|}^2\text{d}{s}_1\text{d}{s}_2\text{d}{t}_1\text{d}{t}_2\\ +\alpha {\text{e}}^{\alpha T^2}[3{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\left|x_0\left({\tau }_1\mathrm{,}{\tau }_2\right)-x_{0_{mm}}\left({\tau }_1\mathrm{,}{\tau }_2\right)\right|}^2\text{d}{\tau }_1\text{d}{\tau }_2\text{d}{t}_1\text{d}{t}_2\end{array}$

$\begin{array}{l}+6l_2^2{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\displaystyle {\int }_0^{{\tau }_2}}{\displaystyle {\int }_0^{{\tau }_1}}{\left|\tilde{k}\left({\tau }_1\mathrm{,}{\tau }_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)-{\tilde{k}}_{mm}\left({\tau }_1\mathrm{,}{\tau }_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right|}^2\text{d}{s}_1\text{d}{s}_2\text{d}{\tau }_1\text{d}{\tau }_2\text{d}{t}_1\text{d}{t}_2\\ +6l_4^2{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}{\displaystyle {\int }_0^{{\tau }_2}}{\displaystyle {\int }_0^{{\tau }_1}}{\left|\hat{k}\left({\tau }_1\mathrm{,}{\tau }_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)-{\hat{k}}_{mm}\left({\tau }_1\mathrm{,}{\tau }_2\mathrm{,}{s}_1\mathrm{,}{s}_2\right)\right|}^2\text{d}{s}_1\text{d}{s}_2\text{d}{\tau }_1\text{d}{\tau }_2\text{d}{t}_1\text{d}{t}_2]\\ \mathrm{<mo>\; =\; </mo>3}{I}_1+6l_2^2{I}_2+6l_4^2{I}_3+\alpha {\text{e}}^{\alpha T^2}\left[3I_4+6l_2^2{I}_5+6l_4^2{I}_6\right]\mathrm{,}\end{array}$

by using (28) (29), the integrals

$I_i\le c_i{h}^2,i=1,2,\cdots ,6,$

the last equation can be converted into

$\begin{array}{l}{\displaystyle {\int }_0^T}{\displaystyle {\int }_0^T}\mathbb{E}{\left|e_{mm}\left(t_1,t_2\right)\right|}^2\text{d}{t}_1\text{d}{t}_2\\ \le \left[\left(3c_1+6l_2^2{c}_2+6l_4^2{c}_3\right)+\alpha {\text{e}}^{\alpha T^2}\left(3c_4+6l_2^2{c}_5+6l_4^2{c}_6\right)\right]h^2\le O\left(h^2\right).\end{array}$

where $c_i,i=1,2,\cdots ,6$ are independent nonnegative constants.

The proof is completed. □

7. Numerical Examples

In the last section, we give a numerical example which illustrates the feasibility of the above method. The approximation solutions and mean solutions of the equations are shown in Figures 1-4.

Example 1. Consider the following two-dimensional nonlinear stochastic Itô-Volterra integral equation (one-dimensional case can reference to Example 1 in [17] ).

$\begin{array}{c}x\left(t_1\mathrm{,}{t}_2\right)=\frac{1}{10}-{\left(\frac{1}{30}\right)}^2{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}x\left(s_1\mathrm{,}{s}_2\right)\left(1-x^2\left(s_1\mathrm{,}{s}_2\right)\right)\text{d}{s}_1\text{d}{s}_2\\ +\frac{1}{30}{\displaystyle {\int }_0^{t_2}}{\displaystyle {\int }_0^{t_1}}\left(1-x^2\left(s_1\mathrm{,}{s}_2\right)\right)\text{d}B\left(s_1\right)\text{d}B\left(s_2\right)\mathrm{.}\end{array}$

The front view and the top view of the approximation solutions of the Example 1 for m = 8 are given in Figure 1.

The front view and the top view of the mean solutions of the Example 1 for m = 8 are given in Figure 2.

The front view and the top view of the approximation solutions of the Example 1 for m = 16 are given in Figure 3.

The front view and the top view of the mean solutions of the Example 1 for m = 16 are given in Figure 4.

From these figures, we find the general trends of the solutions are similar for different m, and the absolute error of mean solution is very small. This method is efficient and the accuracy is credible.

8. Conclusion

For some stochastic Volterra integral equations, exact solutions cannot be expressed. But, the numerical solution can be conveniently obtained based on different stochastic numerical methods. As the complexity of the system, we use

Figure 1. The front view and top view of the approximation solutions for m = 8.

Figure 2. The front view and top view of the mean solutions for m = 8.

Figure 3. The front view and top view of the approximation solutions for m = 16.

Figure 4. The front view and top view of the mean solutions for m = 16.

BPFs as the basis function to solve the two-dimensional nonlinear stochastic Volterra integral equation. This numerical method is simple and effective. In the future, we will try to extend it to n-dimensional space and solve more problems.

Acknowledgements

We thank the Editors and the Reviewers for their helps and comments. This article is funded by NSF Grants 11471105 of China, NSF Grants 2016CFB526 of Hubei Province, Innovation Team of the Educational Department of Hubei Province T201412, and Innovation Items of Hubei Normal University 2018032 and 2018105. These supports are greatly appreciated.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hanson, R. and Phillips, J. (1978) Numerical Solution of Two-Dimensional Integral Equations Using Linear Elements. SIAM Journal on Numerical Analysis, 15, 113-121. https://doi.org/10.1137/0715007
[2]	Babolian, E., Maleknejad, K. and Roodaki, M. (2010) Two-Dimensional Triangular Functions and Their Applications to Nonlinear 2D Volterra Fredholm Integral Equations. Computers & Mathematics with Applications, 60, 1711-1722. https://doi.org/10.1016/j.camwa.2010.07.002
[3]	Fallahpour, M., Khodabin, M. and Maleknejad, K. (2015) Approximation Solution of Two-Dimensional Linear Stochastic Volterra Integral Equations by Applying the Haar Wavelets. International Journal of Mathematical Modelling and Computations (IJM2C), 5, 361-372.
[4]	Jiang, Z.H. and Schaufelberger, W. (1992) Block Pulse Functions and Their Applications in Control Systems. Spriger-Verlag, Berlin.
[5]	Mirzaee, F. and Hadadiyan, E. (2014) Using Modified Two-Dimensional Block-Pulse Functions for the Numerical Solution of Nonlinear Two-Dimensional Volterra Integral Equations. Journal of Hyperstructures, 3, 68-80.
[6]	Nemati, S., Lima, P.M. and Ordokhani, Y. (2013) Numerical Solution of a Class of Two-Dimensional Nonlinear Volterra Integral Equations Using Legendre Polynomials. Journal of Computational and Applied Mathematics, 242, 53-69. https://doi.org/10.1016/j.cam.2012.10.021
[7]	Fallahpour, M., Khodabin, M. and Maleknejad, K. (2016) Approximation Solution of Two-Dimensional Linear Stochastic Volterra-Fredholm Integral Equation via Two-Dimensional Block-Pulse Functions. International Journal of Industrial Mathematics, 8, Article ID: IJIM-00774.
[8]	Maleknejad, K., Sohrabi, S. and Baranji, B. (2010) Application of 2D-BPFs to Nonlinear Integral Equations. Communications in Nonlinear Science and Numerical Simulation, 15, 527-535. https://doi.org/10.1016/j.cnsns.2009.04.011
[9]	Maleknejad, K. and Jafaribehbahani, Z. (2012) Applications of Two-Dimensional Triangular Functions for Solving Nonlinear Class of Mixed Volterra-Fredholm Integral Equations. Mathematical and Computer Modelling, 55, 1833-1844. https://doi.org/10.1016/j.mcm.2011.11.041
[10]	Mirzaee, F. and Hadadiyan, E. (2012) Approximate Solutions for Mixed Nonlinear Volterra-Fredholm Type Integral Equations via Modified Block-Pulse Functions. Journal of the Association of Arab Universities for Basic and Applied Sciences, 12, 65-73. https://doi.org/10.1016/j.jaubas.2012.05.001
[11]	Aleknejad, K., Khodabin, M. and Shekarabi, F.H. (2014) Modified Block Pulse Functions for Numerical Solution of Stochastic Volterra Integral Equations. Journal of Applied Mathematics, 2014, Article ID: 469308.
[12]	Maleknejad, K., Khodabin, M. and Rostami, M. (2012) Numerical Solution of Stochastic Volterra Integral Equation by a Stochastic Operational Matrix Based on Block Pulse Function. Mathematical and Computer Modelling, 55, 791-800. https://doi.org/10.1016/j.mcm.2011.08.053
[13]	Langville, A.N. and Stewart, W.J. (2004) The Kronecker product and Stochastic Automata Networks. Journal of Computational and Applied Mathematics, 167, 429-447. https://doi.org/10.1016/j.cam.2003.10.010
[14]	Moler, C.B. (2006) Numerical Computing with MATLAB. China Machine Press, Beijing.
[15]	Ezzati, R., Khodabin, M. and Sadati, Z. (2014) Numerical Implementation of Stochastic Operational Matrix Driven by a Fractional Brownian Motion for Solving a Stochastic Differential Equation. Abstract and Applied Analysis, 2014, Article ID: 523163.
[16]	Maleknejad, K., Khodabin, M. and Rostami, M. (2012) A Numerical Method for Solving m-Dimensional Stochastic Itô-Volterra Integral Equations by Stochastic Operational Matrix. Computers & Mathematics with Applications, 63, 133-143. https://doi.org/10.1016/j.camwa.2011.10.079
[17]	Mirzaee, F. and Samadyar, N. (2018) Numerical Solution of Nonlinear Stochastic Itô-Volterra Integral Equations Driven by Fractional Brownian Motion. Mathematical Methods in the Applied Sciences, 14, 1410-1423. https://doi.org/10.1002/mma.4671