Numerical approximation of nonlinear stochastic Volterra integral equation based on Walsh function

Abstract

This paper adopts a highly effective numerical approach for approximating non-linear stochastic Volterra integral equations (NLSVIEs) based on the operational matrices of the Walsh function and the collocation method. The method transforms the integral equation into a system of algebraic equations, which allows for the derivation of an approximate solution. Error analysis is performed, confirming the effectiveness of the proposed method, which results in a linear order of convergence. Numerical examples are provided to illustrate the precision and effectiveness of the proposed method.

Introduction

Stochastic integral equations have found extensive use across various disciplines, including financial mathematics, mechanics, biology, engineering, and many other sectors [1, 2]. Since it is not always possible to obtain an exact solution to a problem, numerical approximation to the solution of the integral equation becomes essential. In order to approximate the solution of an integral equation, orthogonal functions such as the block pulse function, Hermite wavelet [3], Legendre polynomial [4], Chebyshev wavelet [5], meshless method [6, 7], finite difference [8, 9] and certain hybrid functions [10] have been applied to approximate the solutions of such integral equations. In recent decades, non-linear stochastic Volterra integral equations (NLSVIEs) have found various applications in the biological sciences, particularly in modelling and simulating biological systems. NLSVIEs have been used for modelling population growth and extinction in various species. For example, they have been used to study the dynamics of predator–prey systems, where the population of predators and prey interact with each other in a non-linear way. NLSVIEs have been used to investigate the spread of infectious diseases in populations and model the dynamics of disease transmission, including the rate of infection and recovery and the effect of different control measures. Also, NLSVIEs have been used to model the dynamics of neuronal networks and synaptic plasticity in the brain and the non-linear interactions between neurons. However, it is impossible to derive exact solutions for all NLSVIEs, various numerical techniques are adopted to obtain approximate solutions. In particular, numerical techniques, including orthogonal functions, are often applied to solve these problems. Orthogonal functions such as block pulse function (BPF) [11], Haar wavelet [12], and Legendre polynomials [13] have been implemented to simulate the solution of NLSVIEs. The Walsh function has also been employed for solving the stochastic Volterra integral equation [14].

In this work, we investigate the approximated solution x(t) of the NLSVIE by using the Walsh function [15] as

1

$$\begin{array}{c}\hfill x(t)=x_0+{\int }_0^t{k}_1(s,t)\beta (x(s))\text{d}s+{\int }_0^t{k}_2(s,t)\sigma (x(s))\text{d}B(s)\end{array}$$where x(t) s an unknown function to be determined and $k_1(s,t)$, $k_2(s,t)$ for $s,t\in [0,T)$, represent the stochastic processes based on the same probability space $(\mathrm{\Omega },F,P)$. Here, the notations B(t) and ${\int }_0^t{k}_2(s,t)x(s)\text{d}B(s)$ denote the Brownian motion and It$\widehat{o}$ integral, respectively [1, 2].

The rest of the current research paper is organised as follows: Sect. 2 introduces the Walsh function and its properties and describes the relationship between the Walsh function and block pulse functions (BPFs). Section 3 presents a numerical technique based on the operational matrices of the Walsh function and the collocation method to discretize the NLSVIEs. Section 4 discusses the convergence and error analysis of the method to demonstrate the method’s validity and precision. Section 5 gives two numerical examples by using the proposed method to demonstrate the efficacy of the method. Finally, Sect. 6 contains some concluding remarks and summarises the main findings.

Walsh function and its properties

Definition 1

Rademacher function $r_i(t)$, $i=1,2,\dots$, for $t\in [0,1)$ is defined by [15]$$\begin{array}{c}\hfill r_i(t)=\left\{\begin{array}{c}1i=0,\hfill \\ \text{sgn}(sin(2^i\pi t))\text{otherwise}\hfill \end{array}\right)\end{array}$$where,$$\begin{array}{c}\hfill \text{sgn}(x)=\left\{\begin{array}{cc}1\hfill & x>0,\hfill \\ 0\hfill & x=0,\hfill \\ -1\hfill & x<0.\hfill \end{array}\right)\end{array}$$

Definition 2

The $n^{\mathit{th}}$ Walsh function for $n=0,1,2,\dots ,$ denoted by $w_n(t)$, $t\in [0,1)$ is defined [15] as$$\begin{array}{c}\hfill w_n(t)={(r_q(t))}^{b_q}.{(r_{q-1}(t))}^{b_{q-1}}.{(r_{q-2}(t))}^{b_{q-2}}\dots {(r_1(t))}^{b_1}\end{array}$$in which $n=b_q{2}^{q-1}+b_{q-1}{2}^{q-2}+b_{q-2}{2}^{q-3}+\dots +b_1{2}^0$ is the binary expression of n. So, q, the number of digits present in the binary expression of n can be obtained by $q=[{log}_{2}n]+1$ in which $[·]$ is the greatest integer less than or equal’s to ${}^{\prime }{·}^{\prime }$.

The first m Walsh functions for $m\in \mathbb{N}$ can be restated as an m-vector by $W(t)={\left[\begin{array}{ccc}{w}_0(t)& w_1(t)& w_2(t)\dots w_{m-1}(t)\end{array}\right]}^T;$$t\in [0,1)$. The Walsh functions satisfy the following properties.

Orthonormality

The set of Walsh functions is orthonormal. i.e.,$$\begin{array}{c}\hfill {\int }_0^1{w}_i(t)w_j(t)\text{d}t=\left\{\begin{array}{cc}1\hfill & i=j,\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}\right)\end{array}$$

Completeness

For every $f\in L^2([0,1))$$$\begin{array}{c}\hfill {\int }_0^1{f}^2(t)\text{d}t=\sum _{i=0}^{\infty }{f}_i^2||w_i(t){||}^2,\end{array}$$in which $f_i={\int }_0^{1}f(t)w_i(t)\text{d}t$.

Walsh function approximation

Any real-valued function $f(t)\in L^2([0,1))$ can be approximated as$$\begin{array}{c}\hfill f_m(t)\simeq \sum _{i=0}^{m-1}{c}_i{w}_i(t),\end{array}$$where, $c_i={\int }_0^{1}f(t)w_i(t)\text{d}t$ and $m=2^k$, $k=0,1,2\dots$

The matrix form can be represented by

2

$$\begin{array}{c}\hfill f(t)\simeq F^T{T}_{W}W(t),\end{array}$$where $F={\left[\begin{array}{ccc}{f}_0& f_1& f_2\dots f_{m-1},\end{array}\right]}^T$, $f_i={\int }_{\mathit{ih}}^{(i+1)h}f(s)\text{d}s$ and $h=\frac{1}{m}$.

Here, $T_W=[w_i({\eta }_j)]$ is named as the Walsh operational matrix where ${\eta }_j\in (jh,(j+1)h)$.

Similarly, function $k(s,t)\in L^2([0,1)\times [0,1))$ can be approximated by$$\begin{array}{c}\hfill k_m(s,t)\simeq \sum _{i=0}^{m-1}\sum _{j=0}^{m-1}{c}_{\mathit{ij}}{w}_i(s)w_j(t)\end{array}$$in which $c_{\mathit{ij}}={\int }_0^1{\int }_0^{1}k(s,t)w_i(s)w_j(t)\text{d}t\text{d}s$ with the matrix form represented by

3

$$\begin{array}{c}\hfill k(s,t)\simeq W^T(s)T_{W}K{T}_{W}W(t)=W^T(t)T_W{K}^T{T}_{W}W(s),\end{array}$$in which $K={[k_{\mathit{ij}}]}_{m\times m},k_{\mathit{ij}}={\int }_{\mathit{ih}}^{(i+1)h}{\int }_{\mathit{jh}}^{(j+1)h}k(s,t)\text{d}t\text{d}s$.

Definition 3

For a fixed positive integer m, an m-set of BPFs ${\varphi }_i(t),t\in [0,1)$ for $i=0,1,...,m-1$ can be defined as$$\begin{array}{c}\hfill {\varphi }_i(t)=\left\{\begin{array}{cc}1\hfill & \text{if}\frac{i}{m}\le t<\frac{(i+1)}{m},\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}\right)\end{array}$$in which ${\varphi }_i$ is known as the ith BPF.

The set of all m BPFs are stated as an m-vector, $\mathrm{\Phi }(t)={\left[\begin{array}{ccc}{\varphi }_0(t)& {\varphi }_1(t)& {\varphi }_2(t)\dots {\varphi }_{m-1}(t)\end{array}\right]}^T$, $t\in [0,1)$. The BPFs are disjoint, complete, and orthogonal [16]. The BPFs in vector form satisfy$$\begin{array}{c}\hfill \mathrm{\Phi }(t)\mathrm{\Phi }{(t)}^{T}X=\tilde{X}\mathrm{\Phi }(t)\text{and}{\mathrm{\Phi }}^T(t)A\mathrm{\Phi }(t)=\widehat{A}\mathrm{\Phi }(t)\end{array}$$in which $X\in {\mathbb{R}}^{m\times 1},\tilde{X}$ is the $m\times m$ diagonal matrix with $\tilde{X}(i,i)=X(i)\text{for}i=1,2,3,\dots m,A\in {\mathbb{R}}^{m\times m}$ and $\widehat{A}={\left[\begin{array}{cccc}{a}_{11}& a_{22}& \dots & a_{\mathit{mm}}\end{array}\right]}^T$ is the m-vector with elements equal to the diagonal entries of A. The integration of BPF vector $\mathrm{\Phi }(t)$, $t\in [0,1)$ can be performed by [17]

4

$$\begin{array}{c}\hfill {\int }_0^t\mathrm{\Phi }(\tau )\text{d}\tau =P\mathrm{\Phi }(t),t\in [0,1),\end{array}$$Consequently, the integral of every function $f(t)\in L^2[0,1)$ can be estimated as$$\begin{array}{c}\hfill {\int }_0^{t}f(s)\text{d}s=F^{T}P\mathrm{\Phi }(t).\end{array}$$The integration of the BPF vector $\mathrm{\Phi }(t)$, with $t\in [0,1)$, via the It$\widehat{o}$ integral can be executed by [18]

5

$$\begin{array}{c}\hfill {\int }_0^t\mathrm{\Phi }(\tau )\text{d}B(\tau )=P_S\mathrm{\Phi }(t),t\in [0,1).\end{array}$$Hence, the It$\widehat{o}$ integral of every function $f(t)\in L^2[0,1)$ can be represented as$$\begin{array}{c}\hfill {\int }_0^{t}f(s)\text{d}B(s)=F^T{P}_S\mathrm{\Phi }(t).\end{array}$$The next theorem elucidates a correlation between the Walsh function and the BPF.

Theorem 1

(See [14].) Let the m-set of Walsh function and BPF vectors be W(t) and $\mathrm{\Phi }(t)$, respectively. Then the BPF vectors $\mathrm{\Phi }(t)$ can be used to approximate W(t) as $W(t)=T_W\mathrm{\Phi }(t)$, $m=2^k$, and $k=0,1,\dots$, where $T_W=[c_{\mathit{ij}}{]}_{m\times m}$, $c_{\mathit{ij}}=w_i({\eta }_j)$, for some ${\eta }_j\in (\frac{j}{m},\frac{j+1}{m})$ and $i,j=0,1,2,\dots m-1$.

One can see that [19]$$\begin{array}{c}\hfill T_W{T}_W^T=mI\text{and}{T}_W^T=T_W\end{array}$$which implies that $\mathrm{\Phi }(t)=\frac{1}{m}{T}_{W}W(t).$

Lemma 1

(See [14].) Assume that W(t) represent a Walsh function vector, then the integral of W(t) w.r.t. t is defined by ${\int }_0^{t}W(s)\text{d}s=\wedge W(t)$, where $\wedge =\frac{1}{m}{T}_{W}P{T}_W$ and$$\begin{array}{c}\hfill P=\frac{1}{h}\left[\begin{array}{ccccc}1& 2& 2& \dots & 2\\ 0& 1& 2& \dots & 2\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 0& 0& 0& \dots & 1\end{array}\right].\end{array}$$

Lemma 2

(See [14].) Assume that W(t) denote a Walsh function vector, then the It$\widehat{o}$ integral of W(t) is defined by ${\int }_0^{t}W(s)\text{d}B(s)={\wedge }_{S}W(t)$, where ${\wedge }_S=\frac{1}{m}{T}_W{P}_S{T}_W$ and$$\begin{array}{c}\hfill P_S=\left[\begin{array}{cccc}B(\frac{h}{2})& B(h)& \dots & B(h)\\ 0& B(\frac{3h}{2})-B(h)& \dots & B(2h)-B(h)\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & B(\frac{(2m-1)h}{2})-B((m-1)h)\end{array}\right].\end{array}$$

Numerical method of the NLSVIE

Let us consider the NLSVIE as

6

$$\begin{array}{c}\hfill x(t)=x_0+{\int }_0^t{k}_1(s,t)\beta (x(s))\text{d}s+{\int }_0^t{k}_2(s,t)\sigma (x(s))\text{d}B(s),\end{array}$$in which x(t), $k_1(s,t)$, $k_2(s,t)$ for $s,t\in [0,T)$, represent the stochastic processes based on the same probability space $(\mathrm{\Omega },F,P)$ and x(t) is unknown. Here B(t) is Brownian motion [1, 2] and ${\int }_0^t{k}_2(s,t)x(s)\text{d}B(s)$ is the It$\widehat{o}$ integral.

Let $z_1(t)=\beta (x(s))$ and $z_2(t)=\sigma (x(s))$ which implies,

7

$$\begin{array}{c}\hfill x(t)=x_0+{\int }_0^t{k}_1(s,t)z_1(s)\text{d}s+{\int }_0^t{k}_2(s,t)z_2(s)\text{d}B(s).\end{array}$$The function $z_1(t)$, $z_2(t)$, $k_1(s,t)$, $k_2(s,t)$ can be approximated for $s,t\in [0,T)$ as

8

$$\begin{array}{c}\hfill z_{\eta }(t)\simeq Z_{\eta }^T{T}_{W}W(t),\eta =1,2,\end{array}$$where $Z_{\eta }={\left[\begin{array}{ccc}{c}_0& c_1& c_2\dots c_{m-1}\end{array}\right]}^T$ and $c_i={\int }_{\mathit{ih}}^{(i+1)h}{z}_{\eta }(s)\text{d}s$.

Similarly, for $\gamma =1,2$

9

$$\begin{array}{c}\hfill k_{\gamma }(s,t)\simeq W^T(s)T_W{K}_{\gamma }{T}_{W}W(t)=W^T(t)T_W{K}_{\gamma }^T{T}_{W}W(s)\end{array}$$in which $K_{\gamma }={[k_{\mathit{ij}}]}_{m\times m},k_{\mathit{ij}}={\int }_{\mathit{ih}}^{(i+1)h}{\int }_{\mathit{jh}}^{(j+1)h}{k}_{\gamma }(s,t)\text{d}t\text{d}s$.

Assume that

10

$$\begin{array}{c}\hfill x(t)\simeq X^T{T}_{W}W(t),\end{array}$$in which $X={\left[\begin{array}{ccc}{x}_0& x_1& x_2\dots x_{m-1}\end{array}\right]}^T$ and $x_i={\int }_{\mathit{ih}}^{(i+1)h}x(s)\text{d}s$. Replacing (8), (10) and (9) in (7) leads to

11

$$\begin{array}{cc}\hfill X^T{T}_{W}W(t)\simeq & x_0+{\int }_0^t{W}^T(t)T_W{K}_1^T{T}_{W}W(s)W^T(s)T_W{Z}_1\text{d}s\hfill \\ \hfill & +{\int }_0^t{W}^T(t)T_W{K}_2^T{T}_{W}W(s)W^T(s)T_W{Z}_2\text{d}B(s)\hfill \\ \hfill =& x_0+W^T(t)T_W{K}_1^T{T}_W{\int }_0^{t}W(s)W^T(s)T_W{Z}_1\text{d}s\hfill \\ \hfill & +W^T(t)T_W{K}_2^T{T}_W{\int }_0^{t}W(s)W^T(s)T_W{Z}_2\text{d}B(s).\hfill \end{array}$$Now

12

$$\begin{array}{c}\hfill {\int }_0^{t}W(s)W^T(s)T_W{Z}_1\text{d}s=T_W\widetilde{Z_1}P{T}_{W}W(t).\end{array}$$Similarly,

13

$$\begin{array}{c}\hfill {\int }_0^{t}W(s)W^T(s)T_W{Z}_2\text{d}B(s)=T_W\widetilde{Z_2}{P}_S{T}_{W}W(t).\end{array}$$Inserting (12) and (13) in (11), we obtain$$\begin{array}{cc}\hfill X^T{T}_{W}W(t)=& x_0+m{W}^T(t)T_W{K}_1^T\widetilde{Z_1}P{T}_{W}W(t)\hfill \\ \hfill & +m{W}^T(t)T_W{K}_2^T\widetilde{Z_2}{P}_S{T}_{W}W(t)\hfill \\ \hfill =& x_0+W^T(t)T_W{H}_1{T}_{W}W(t)+W^T(t)T_W{H}_2{T}_{W}W(t)\hfill \\ \hfill =& x_0+m{\widehat{H_1}}^T{T}_{W}W(t)+m{\widehat{H_2}}^T{T}_{W}W(t)\hfill \end{array}$$i.e.,

14

$$\begin{array}{c}\hfill (X^T-m{\widehat{H_1}}^T-m{\widehat{H_2}}^T)T_{W}W(t)=x_0\end{array}$$and

15

$$\begin{array}{c}\hfill Z_1^T{T}_{W}W(t)=\beta (x_0+m{\widehat{H_1}}^T{T}_{W}W(t)+m{\widehat{H_2}}^T{T}_{W}W(t)),\\ \hfill Z_2^T{T}_{W}W(t)=\sigma (x_0+m{\widehat{H_1}}^T{T}_{W}W(t)+m{\widehat{H_2}}^T{T}_{W}W(t)),\end{array}$$in which $H_1=m{K}_1^T\widetilde{Z_1}P$, $H_2=m{K}_2^T\widetilde{Z_2}{P}_S$ and ${\widehat{H}}_i$ denotes the m-vector with elements equal to the diagonal entries of $H_i$.

To calculate $Z_1$ and $Z_2$, we collocate the aforementioned Eq. (15) at $t_j=\frac{2j+1}{2m}$ for $j=0,1,\dots ,m-1$ and solve the following system

16

$$\begin{array}{c}\hfill Z_1^T{T}_{W}W(t_j)=\beta (x_0+m{\widehat{H_1}}^T{T}_{W}W(t_j)+m{\widehat{H_2}}^T{T}_{W}W(t_j)),\\ \hfill Z_2^T{T}_{W}W(t_j)=\sigma (x_0+m{\widehat{H_1}}^T{T}_{W}W(t_j)+m{\widehat{H_2}}^T{T}_{W}W(t_j)).\end{array}$$

Error analysis of the method

This section focuses on analysing the discrepancy between the approximate and exact solutions of the NLSVIE. Prior to commencing the analysis, we define the notation ${E(|X|}^2{)}^{\frac{1}{2}}={\Vert X\Vert }_2$.

Theorem 2

(See [14].) If $f\in L^2[0,1)$ and fulfills the Lipschitz condition with a Lipschitz constant C, then the 2-norm of $e_m(t)$ is $\mathcal{O}(h)$, where $e_m(t)=|f(t)-\sum _{i=0}^{m-1}{c}_i{w}_i(t)|$ and $c_i={\int }_0^{1}f(s)w_i(s)\text{d}s$.

Theorem 3

(See [14].) Assume that $k\in L^2([0,1)\times [0,1))$ fulfills the Lipschitz condition with Lipschitz constant L. If $k_m(x,y)={\sum }_{i=0}^{m-1}{\sum }_{j=0}^{m-1}{c}_{\mathit{ij}}{w}_i(x)w_j(y)$, $c_{\mathit{ij}}={\int }_0^1{\int }_0^{1}k(s,t)w_i(s)w_j(t)\text{d}t\text{d}s$, then $\Vert e_m{(x,y)\Vert }_2=\mathcal{O}(h)$, where $|e_m(x,y)|=|k(x,y)-k_m(x,y)|$.

Theorem 4

Assume that $x_m(t)$ be the approximated solution of the NLSIE (7). If

$f\in L^2[0,1)$, $k_1(s,t)$ and $k_2(s,t)\in L^2([0,1)\times [0,1))$ fulfills the Lipschitz condition with Lipschitz constants with Lipschitz constants C, $L_1$ and $L_2$ respectively

$|k_1(s,t)|\le {\rho }_1$, $|k_2(s,t)|\le {\rho }_2$, $|\beta (x(t))|\le {\zeta }_1$ and $|\sigma (x(t))|\le {\zeta }_2$, and

for ${\eta }_1,{\eta }_2\ge 0$, $|\beta (x(t))-\beta (x_m(t))|\le {\eta }_1|x(t)-x_m(t)|$, $|\sigma (x(t))-\sigma (x_m(t))|\le {\eta }_2|x(t)-x_m(t)|$.

Then, we have $\Vert x(t)-x_m{(t)\Vert }_2^2=O(h^2)$.

Proof

Consider the NLSIE (7) and let $x_m(t)$ be the approximation of the solution obtained using the Walsh function. Then, we have$$\begin{array}{cc}\hfill x(t)-x_m(t)=& f(t)-f_m(t)\hfill \\ \hfill & +{\int }_0^t(k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s)))\text{d}s\hfill \\ \hfill & +{\int }_0^t(k_2(s,t)\sigma (x(s))-k_{2m}(s,t)\sigma (x_m(s)))\text{d}B(s)\hfill \end{array}$$or,$$\begin{array}{cc}\hfill |x(t)-x_m(t)|\le & |f(t)-f_m(t)|\hfill \\ \hfill & +|{\int }_0^t(k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s)))\text{d}s|\hfill \\ \hfill & +|{\int }_0^t(k_2(s,t)\sigma (x(s))-k_{2m}(s,t)\sigma (x_m(s)))\text{d}B(s)|.\hfill \end{array}$$We know that, ${(a+b+c)}^2\le 3a^2+3b^2+3c^2$$$\begin{array}{cc}\hfill |x(t)-x_m{(t)|}^2\le & 3|f(t)-f_m{(t)|}^2\hfill \\ \hfill & +3|{\int }_0^t(k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s)))\text{d}s{|}^2\hfill \\ \hfill & +3|{\int }_0^t(k_2(s,t)\sigma (x(s))-k_{2m}(s,t)\sigma (x_m(s)))\text{d}B(s){|}^2.\hfill \end{array}$$which implies that

17

$$\begin{array}{cc}\hfill E(|x(t)-x_m{(t)|}^2)\le & E(3|f(t)-f_m{(t)|}^2)\hfill \\ \hfill & +E(3|{\int }_0^t(k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s)))\text{d}s{|}^2)\hfill \\ \hfill & +E(3|{\int }_0^t(k_2(s,t)\sigma (x(s))-k_{2m}(s,t)\sigma (x_m(s)))\text{d}B(s){|}^2).\hfill \\ \hfill \end{array}$$Suppose,$$\begin{array}{c}\hfill I_1=E(3|{\int }_0^t(k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s)))\text{d}s{|}^2)\end{array}$$and$$\begin{array}{c}\hfill I_2=E(3|{\int }_0^t(k_2(s,t)\sigma (x(s))-k_{2m}(s,t)\sigma (x_m(s)))\text{d}B(s){|}^2).\end{array}$$Applying the Theorem 2 in inequality (17) results

18

$$\begin{array}{c}\hfill E(|x(t)-x_m{(t)|}^2)\le C^2{h}^2+I_1+I_2\end{array}$$Now,$$\begin{array}{cc}\hfill |k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s))|=& |k_1(s,t)\beta (x(s))-k_1(s,t)\beta (x_m(s))\hfill \\ \hfill & +k_1(s,t)\beta (x_m(s))-k_{1m}(s,t)\beta (x_m(s))\Vert \hfill \\ \hfill \le & |k_1(s,t)||\beta (x(s))-\beta (x_m(s))|\hfill \\ \hfill & +|k_1(s,t)-k_{1m}(s,t)||\beta (x_m(s))|\hfill \\ \hfill =& |k_1(s,t)||\beta (x(s))-\beta (x_m(s))|\hfill \\ \hfill & +|k_1(s,t)-k_{1m}(s,t)||\beta (x(s))-\beta (x(s))+\beta (x_m(s))|\hfill \\ \hfill \le & |k_1(s,t)||\beta (x(s))-\beta (x_m(s))|\hfill \\ \hfill & +|k_1(s,t)-k_{1m}(s,t)||\beta (x(s))|\hfill \\ \hfill & +|k_1(s,t)-k_{1m}(s,t)||\beta (x(s))-\beta (x_m(s))|\hfill \\ \hfill \end{array}$$Let $|k_1(s,t)|\le {\rho }_1$, $|\beta (x(t))|\le {\zeta }_1$, $|\beta (x(t))-\beta (x_m(t))|\le {\eta }_1|x(t)-x_m(t)|$ and using Theorem 3, we get$$\begin{array}{cc}\hfill |k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s))|\le & {\rho }_1{\eta }_1|x(s)-x_m(s)|\hfill \\ \hfill & +\sqrt{2}{L}_{1}h{\zeta }_1+\sqrt{2}{L}_{1}h{\eta }_1|x(s)-x_m(s)|\hfill \end{array}$$which arrives at

19

$$\begin{array}{c}\hfill |k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s))|\le \sqrt{2}{L}_{1}h{\zeta }_1+({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)|x(s)-x_m(s)|\\ \hfill \end{array}$$which gives,$$\begin{array}{cc}\hfill I_1\le & E(3({\int }_0^t|k_1(s,t)\beta (x(s))-k_{1m}(s,t)\beta (x_m(s))|\text{d}s{)}^2)\hfill \\ \hfill \le & E(3({\int }_0^t(\sqrt{2}{L}_{1}h{\zeta }_1\hfill \\ \hfill & +({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)|x(s)-x_m(s)|)\text{d}s{)}^2)\hfill \end{array}$$In virtue of the Cauchy–Schwarz inequality, for $t>0$ and $f\in L^2[0,t)$$$\begin{array}{c}\hfill |{\int }_0^{t}f(s)\text{d}s{|}^2\le t{\int }_0^t{|f|}^2\text{d}s\end{array}$$then, we can write$$\begin{array}{cc}\hfill I_1\le & E(3{\int }_0^t(\sqrt{2}{L}_{1}h\zeta +({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)|x(s)-x_m(s)|{)}^2\text{d}s)\hfill \\ \hfill \le & E(6{\int }_0^t((\sqrt{2}{L}_{1}h{\zeta }_1{)}^2+(({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)|x(s)-x_m(s)|{)}^2)\text{d}s)\hfill \\ \hfill \le & E(6(\sqrt{2}{L}_{1}h{\zeta }_1{)}^2+6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2{\int }_0^t{|x(s)-x_m(s)|}^2\text{d}s)\hfill \end{array}$$Hence,

20

$$\begin{array}{cc}\hfill I_1\le & 6(\sqrt{2}{L}_{1}h{\zeta }_1{)}^2\hfill \\ \hfill & +6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2{\int }_0^{t}E({|x(s)-x_m(s)|}^2)\text{d}s.\hfill \end{array}$$Now,$$\begin{array}{cc}\hfill I_2\le & E(3{\int }_0^t|k_2(s,t)\sigma (x(s))-k_{2m}(s,t)\sigma (x_m(s)){|}^2\text{d}s)\hfill \end{array}$$Let $|k_2(s,t)|\le {\rho }_2$, $|\sigma (x(t))|\le {\zeta }_2$, $|\sigma (x(t))-\sigma (x_m(t))|\le {\eta }_2|x(t)-x_m(t)|$ and using Theorem 3, we get$$\begin{array}{cc}\hfill & |k_2(s,t)\sigma (x(s))-k_{2m}(s,t)\sigma (x_m(s))|\hfill \\ \hfill & \le \sqrt{2}{L}_{2}h{\zeta }_2+({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)|x(s)-x_m(s)|\hfill \\ \hfill & I_2\le E(3{\int }_0^t(\sqrt{2}{L}_{2}h{\zeta }_2+({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)|x(s)-x_m(s)|{)}^2\text{d}s)\hfill \\ \hfill & \le E(6{\int }_0^t({(\sqrt{2}{L}_{2}h{\zeta }_2)}^2+{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2{|x(s)-x_m(s)|}^2)\text{d}s)\hfill \\ \hfill & =E(6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2+6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2{\int }_0^t{|x(s)-x_m(s)|}^2\text{d}s).\hfill \end{array}$$Hence,

21

$$\begin{array}{cc}\hfill I_2\le & 6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2\hfill \\ \hfill & +6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2{\int }_0^{t}E({|x(s)-x_m(s)|}^2)\text{d}s\hfill \end{array}$$Employing (20) and (21) in (18), we have$$\begin{array}{cc}\hfill E(|x(t)-x_m{(t)|}^2)\le & C^2{h}^2+6(\sqrt{2}{L}_{1}h{\zeta }_1{)}^2\hfill \\ \hfill & +6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2{\int }_0^{t}E({|x(s)-x_m(s)|}^2)\text{d}s\hfill \\ \hfill & +6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2+6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2{\int }_0^{t}E({|x(s)-x_m(s)|}^2)\text{d}s\hfill \end{array}$$which concludes that

22

$$\begin{array}{cc}\hfill E(|x(t)-x_m{(t)|}^2)& \le (C^2{h}^2+6(\sqrt{2}{L}_{1}h{\zeta }_1{)}^2+6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2)\hfill \\ \hfill & +(6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2+6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2)\hfill \\ \hfill & \times {\int }_0^{t}E({|x(s)-x_m(s)|}^2)\text{d}s.\hfill \end{array}$$Let$$\begin{array}{cc}\hfill R_1=& (C^2{h}^2+6(\sqrt{2}{L}_{1}h{\zeta }_1{)}^2+6{(\sqrt{2}{L}_{2}h{\zeta }_2)}^2),\hfill \\ \hfill R_2=& (6{({\rho }_1{\eta }_1+\sqrt{2}{L}_{1}h{\eta }_1)}^2+6{({\rho }_2{\eta }_2+\sqrt{2}{L}_{2}h{\eta }_2)}^2).\hfill \end{array}$$By using Gronwall’s inequality, we have

23

$$\begin{array}{cc}\hfill E(|x(t)-x_m{(t)|}^2)\le & R_{1}exp({\int }_0^t{R}_2\text{d}s),\hfill \end{array}$$which holds that,

24

$$\begin{array}{c}\hfill \Vert x(t)-x_m{(t)\Vert }_2^2=E({|x(t)-x_m(t)|}^2)\le R_1{e}^{R_2}=\mathcal{O}(h^2).\end{array}$$$\square$

Numerical examples

This section employs the proposed method to solve the NLSVIE. Two numerical examples are given to demonstrate the convergence of the method by comparing the approximate and analytical results. To measure the error between the two solutions, the infinity norm of the error E is defined as ${\Vert E\Vert }_{\infty }=\underset{1\le i\le m}{max}|X_i-Y_i|$, where $X_i$ and $Y_i$ are the Walsh coefficients of the exact and approximate solutions, respectively. The number of iterations for each example is denoted as n, while the mean and standard deviation of the error E are represented as $\overline{x_E}$ and $s_E$, respectively. All computations are performed using Matlab 2013(a).

Example 1

Let us consider the NLSVIE [1] as$$\begin{array}{c}\hfill x(t)=x_0-\frac{a^2}{2}{\int }_0^{t}tanh(x(s)){\text{sech}}^2(x(s))\text{d}s+a{\int }_0^t\text{sech}(x(s))\text{d}B(s),\end{array}$$where, $x_0=0$ and $a=\frac{1}{30}$ and x(t) represents the unknown stochastic process based on the same probability space $(\mathrm{\Omega },F,P)$, and B(t) is a Brownian motion process. The exact solution is given by $x(t)={sinh}^{-1}(aB(t)+sinh(x_0))$.

Table 1 reports ${\overline{x}}_E$ and $s_E$ errors as well as the interval of confidence for the mean error with $m=16$ and 50 iterations. Figure 1 displays the numerical and exact solutions with $m=32$ and $m=64$. Figure 2 shows the behaviour of error solutions with $m=16$, $m=32$, and $n=50$. Table  2 represents ${\overline{x}}_E$ and $s_E$ errors as well as the interval of confidence for the mean error with $m=32$ and 50 iterations.

Table 1. ${\overline{x}}_E$ and $s_E$ errors as well as interval of confidence for mean error in Example 1 with $m=16$ and 50 iterations

t	${\overline{x}}_E$	$s_E$	95% interval of confidence for error mean
			Lower	Upper
0.1	$0.0073\times {10}^{-5}$	$0.0137\times {10}^{-5}$	$0.0348\times {10}^{-6}$	$0.0111\times {10}^{-5}$
0.3	$0.0330\times {10}^{-5}$	$0.1151\times {10}^{-5}$	$0.0109\times {10}^{-6}$	$0.0649\times {10}^{-5}$
0.5	$0.0900\times {10}^{-5}$	$0.3352\times {10}^{-5}$	$0.0290\times {10}^{-6}$	$0.1829\times {10}^{-5}$
0.7	$0.1402\times {10}^{-5}$	$0.4698\times {10}^{-5}$	$0.0995\times {10}^{-6}$	$0.2704\times {10}^{-5}$
0.9	$0.1939\times {10}^{-5}$	$0.6806\times {10}^{-5}$	$0.0519\times {10}^{-6}$	$0.3825\times {10}^{-5}$

[See PDF for image]

Fig. 1

The numerical and exact solutions of Example 1 with $m=32$ and $m=64$

[See PDF for image]

Fig. 2

The behavior of error solutions in Example 1 with $m=16$, $m=32$, and $n=50$

Table 2. ${\overline{x}}_E$ and $s_E$ errors as well as interval of confidence for mean error in Example 1 with $m=32$ and 50 iterations

t	${\overline{x}}_E$	$s_E$	95% interval of confidence for error mean
			Lower	Upper
0.1	$0.0192\times {10}^{-5}$	$0.0223\times {10}^{-5}$	$0.0130\times {10}^{-5}$	$0.0254\times {10}^{-5}$
0.3	$0.0809\times {10}^{-5}$	$0.1477\times {10}^{-5}$	$0.0400\times {10}^{-5}$	$0.1219\times {10}^{-5}$
0.5	$0.1577\times {10}^{-5}$	$0.3446\times {10}^{-5}$	$0.0621\times {10}^{-5}$	$0.2532\times {10}^{-5}$
0.7	$0.2262\times {10}^{-5}$	$0.4336\times {10}^{-5}$	$0.1061\times {10}^{-5}$	$0.3464\times {10}^{-5}$
0.9	$0.2992\times {10}^{-5}$	$0.5677\times {10}^{-5}$	$0.1418\times {10}^{-5}$	$0.4565\times {10}^{-5}$

Example 2

Consider the NLSVIE [1] as

25

$$\begin{array}{c}\hfill x(t)=x_0-a^2{\int }_0^{t}x(s)(1-x^2(s))\text{d}s+a{\int }_0^t(1-x^2(s))\text{d}B(s),\end{array}$$where, $x_0=\frac{1}{10}$ and $a=\frac{1}{30}$ and x(t) represents the unknown stochastic process based on the same probability space $(\mathrm{\Omega },F,P)$, and B(t) is a Brownian motion process. The exact solution is given by $x(t)=tanh(aB(t)+{tanh}^{-1}(x_0))$.

Table 3 displays ${\overline{x}}_E$ and $s_E$ errors as well as interval of confidence for mean error with $m=16$ and 50 iterations. Table 4 presents ${\overline{x}}_E$ and $s_E$ errors as well as interval of confidence for mean error with $m=32$ and 50 iterations. Figure 3 shows the numerical and exact solutions with $m=16$ and $m=632$.

Table 3. ${\overline{x}}_E$ and $s_E$ error as well as interval of confidence for mean error in Example 2 with $m=16$ and 50 iterations

t	${\overline{x}}_E$	$s_E$	95% interval of confidence for error mean.
			Lower	Upper
0.1	$0.1029\times {10}^{-4}$	$0.0045\times {10}^{-4}$	$0.1016\times {10}^{-4}$	$0.0104\times {10}^{-3}$
0.3	$0.3070\times {10}^{-4}$	$0.0284\times {10}^{-4}$	$0.2992\times {10}^{-4}$	$0.0315\times {10}^{-3}$
0.5	$0.5373\times {10}^{-4}$	$0.0761\times {10}^{-4}$	$0.5163\times {10}^{-4}$	$0.0558\times {10}^{-3}$
0.7	$0.7689\times {10}^{-4}$	$0.1193\times {10}^{-4}$	$0.7358\times {10}^{-4}$	$0.0802\times {10}^{-3}$
0.9	$0.9666\times {10}^{-4}$	$0.1706\times {10}^{-4}$	$0.9193\times {10}^{-4}$	$0.1014\times {10}^{-3}$

Table 4. ${\overline{x}}_E$ and $s_E$ errors as well as interval of confidence for mean error in Example 2 with $m=32$ and 50 iterations

t	${\overline{x}}_E$	$s_E$	95% interval of confidence for error mean
			Lower	Upper
0.1	$0.1184\times {10}^{-4}$	$0.0064\times {10}^{-4}$	$0.1156\times {10}^{-4}$	$0.0121\times {10}^{-3}$
0.3	$0.3161\times {10}^{-4}$	$0.0281\times {10}^{-4}$	$0.3038\times {10}^{-4}$	$0.0328\times {10}^{-3}$
0.5	$0.5215\times {10}^{-4}$	$0.0832\times {10}^{-4}$	$0.4850\times {10}^{-4}$	$0.0557\times {10}^{-3}$
0.7	$0.7284\times {10}^{-4}$	$0.1288\times {10}^{-4}$	$0.6719\times {10}^{-4}$	$0.0784\times {10}^{-3}$
0.9	$0.9216\times {10}^{-4}$	$0.1893\times {10}^{-4}$	$0.8387\times {10}^{-4}$	$0.1004\times {10}^{-3}$

[See PDF for image]

Fig. 3

The numerical and exact solutions of Example 2 with $m=16$ and $m=32$

Conclusion

This paper implemented the proposed numerical approach based on the Walsh function to solve the integral equation, which transforms the problem into a system of algebraic equations. This system was then solved to obtain an approximation of the solution. Through convergence and error analysis, the method’s order of convergence was found to be $\mathcal{O}(h)$. The method was applied to numerical examples with known exact solutions in the final section, and the results are presented in tables and figures. These results demonstrated that the method is highly accurate and efficient, with high precision achievable using a limited number of basis functions. Increasing the number of basis functions reduces absolute errors. The method can be adapted to solve non-linear stochastic integral equations induced by fractional Brownian motion. Additionally, the approximated solution to the stochastic integral equation with a singular kernel can be obtained using the above-mentioned method. The Walsh function can be modified further to enhance the method and attain a higher order of convergence.

Acknowledgements

The authors would like to thank the reviewers for their helpful suggestions and comments, which improved the quality of this paper.

Data availibility

No data was used for the research described in the article.

Declarations

Conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this article.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Kloeden, P.E., Platen, E., Kloeden, P.E., Platen, E.: Stochastic differential equations (Springer, 1992)

2. Øksendal, B., Øksendal, B.: Stochastic differential equations (Springer, 2003)

3. Rostami, Y. Two approximated techniques for solving of system of two-dimensional partial integral differential equations with weakly singular kernels. Comput. Appl. Math.; 2021; 40, 6 217.4301837 [DOI: https://dx.doi.org/10.1007/s40314-021-01608-1]

4. Rostami, Y; Maleknejad, K. Comparison of two hybrid functions for numerical solution of nonlinear mixed partial integro-differential equations. Iran. J. Sci. Technol. Trans. A: Sci.; 2022; 46, 2 645.4410094 [DOI: https://dx.doi.org/10.1007/s40995-022-01277-7]

5. Rostami, Y. Operational matrix of two dimensional Chebyshev wavelets and its applications in solving nonlinear partial integro-differential equations. Eng. Comput.; 2020; 38, 2 745. [DOI: https://dx.doi.org/10.1108/EC-03-2020-0162]

6. Golbabai, A; Nikan, O; Tousi, JR. Note on using radial basis functions method for solving nonlinear integral equations. Commun. Numer. Anal.; 2016; 2016, 2 81.3506409 [DOI: https://dx.doi.org/10.5899/2016/cna-00257]

7. Nikan, O; Molavi-Arabshai, SM; Jafari, H. Numerical simulation of the nonlinear fractional regularized long-wave model arising in ion acoustic plasma waves, Discrete Contin. Dyn. Syst.-S.; 2021; 14, 10 3685.4296169

8. Zheng, X; Qiu, W; Chen, H. Three semi-implicit compact finite difference schemes for the nonlinear partial integro-differential equation arising from viscoelasticity. Int. J. Model. Simul.; 2021; 41, 3 234. [DOI: https://dx.doi.org/10.1080/02286203.2020.1720566]

9. Qiao, L; Qiu, W; Tang, B. A fast numerical solution of the 3D nonlinear tempered fractional integrodifferential equation. Numer. Methods Partial Differ. Equ.; 2023; 39, 2 1333.4563172 [DOI: https://dx.doi.org/10.1002/num.22936]

10. Rostami, Y; Maleknejad, K. The solution of the nonlinear mixed partial integro-differential equation via two-dimensional hybrid functions. Mediterr. J. Math.; 2022; 19, 2 89.4396821 [DOI: https://dx.doi.org/10.1007/s00009-022-01998-4]

11. Asgari, M; Shekarabi, FH. Numerical solution of nonlinear stochastic differential equations using the block pulse operational matrices. Math. Sci.; 2013; 7, 1.3877778 [DOI: https://dx.doi.org/10.1186/2251-7456-7-47]

12. Wu, J; Jiang, G; Sang, X. Numerical solution of nonlinear stochastic itô-Volterra integral equations based on Haar wavelets. Adv. Differ. Equ.; 2019; 2019, 1 1. [DOI: https://dx.doi.org/10.1186/s13662-019-2440-6]

13. Zeghdane, R. Numerical approach for solving nonlinear stochastic Itô-Volterra integral equations using shifted Legendre polynomials. Int. J. Dyn. Syst. Differ. Equ.; 2021; 11, 1 69.4249576

14. Paikaray, P.P., Parida, N.C., Beuria, S.: Numerical approximation of $p$-dimensional stochastic Volterra integral equation using walsh function. J. Math. Comput. Sci. 31(4), 448 (2023)

15. Walsh, JL. A closed set of normal orthogonal functions. Am. J. Math.; 1923; 45, 1 5.1506485 [DOI: https://dx.doi.org/10.2307/2387224]

16. Rao, G.P.: Piecewise constant orthogonal functions and their application to systems and control (Springer, 1983)

17. Hatamzadeh-Varmazyar, S; Masouri, Z; Babolian, E. Numerical method for solving arbitrary linear differential equations using a set of orthogonal basis functions and operational matrix. Appl. Math. Model.; 2016; 40, 1 233.3432084 [DOI: https://dx.doi.org/10.1016/j.apm.2015.04.048]

18. Maleknejad, K; Khodabin, M; Rostami, M. Numerical solution of stochastic Volterra integral equations by a stochastic operational matrix based on block pulse functions. Math. Comput. Model.; 2012; 55, 3–4 791.2887420 [DOI: https://dx.doi.org/10.1016/j.mcm.2011.08.053]

19. Cheng, C; Tsay, Y; Wu, T. Walsh operational matrices for fractional calculus and their application to distributed systems. J. Frankl. Inst.; 1977; 303, 3 267. [DOI: https://dx.doi.org/10.1016/0016-0032(77)90029-1]