# Numerical Solution of Multidimensional Stochastic Itô-Volterra Integral Equation Based on the Least Squares Method and Block Pulse Function

Authors: Ting Ke; Guo Jiang; Mengting Deng
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, a method based on the least squares method and block pulse function is proposed to solve the multidimensional stochastic Itô-Volterra integral equation. The Itô-Volterra integral equation is transformed into a linear algebraic equation. Furthermore, the error analysis is given by the isometry property and Doob’s inequality. Numerical examples verify the effectiveness and precision of this method.


---

## Body

## 1. Introduction

Stochastic Volterra integral equations have been generally used in many fields such as mechanics, medicine, materials, finance, and physical science, especially in the process of engineering practice, because many material coefficients are uncertain, it is necessary to introduce stochastic integral equations to establish models. Meanwhile, the stochastic integral equation is used to solve a variety of engineering problems, for example, material science modeling [1], automated systems science [2], dynamics [3], and biology [4]. However, most of the stochastic Volterra integral equations cannot be solved explicitly. Therefore, it is essential to provide numerical solutions to these equations. Many scholars have developed various methods to study the stochastic Volterra integral equation. Different Volterra integral equations are handled by Bernstein polynomials, block pulse functions, least squares method, Haar wavelets, Walsh functions, Chebyshev and Legendre polynomials, and other functions. We only mentioned the referenced such as [5–21] and other relevant literatures. On the other hand, some authors obtained the numerical solution of stochastic Volterra integral equation by Euler–Maruyama approximation or iterative algorithm, for example [22–27].Recently, there have been many articles on the numerical solution of the stochastic Itô-Volterra integral equation. For example, in [13], Maleknejad et al. considered linear stochastic Itô-Volterra integral equations (SIVIEs) by block pulse functions (BPFs). In [15], Wu et al. used Haar wavelets (HWs) to solve nonlinear SIVIEs. Mirzaee et al. studied nonlinear SIVIEs of fractional order based on the hybrid of block pulse and parabolic functions [16]. In [28], Jiang et al. used BPFs to obtain the numerical solution of two-dimensional nonlinear SIVIEs. Ahmadinia et al. presented a new method that applies the least squares method to study SIVIEs [29]. Moreover, Maleknejad et al. simplified m-dimensional SIVIEs into linear lower trigonometric equations by using BPFs and stochastic integral operation matrix, and then the equations were solved [30].Motivated by the abovementioned literatures, we consider the following multidimensional linear SIVIEs which are studied relatively little and propose a method that combines BPFs and least squares method.(1)Xt=ft+∫0tS˜u,tXudu+∑k=1q∫0tS^ku,tXudBku,t∈0,T,where Xt,S˜u,t, and S^ku,t,k=1,2,…,q, for 0<u<t, are stochastic processes on the same probability space of Ω,ℱ,P, ft is an initial function, and Xt is unknown. Bku are Brownian motions and ∫0tS^ku,tXudBku are Itô integrals.In contrast to the articles [13, 29], this paper’s difference is to study linear SIVIEs driven by multiple independent Brownian motions. In addition, the numerical solution based on the least squares method and BPFs is more accurate than other methods in the reference [30]. Finally, this method can be applied to the models of engineering and material science, which makes the practical data more meaningful.In Section2, BPFs are introduced. In Section 3, stochastic integral operation matrix is given. In Section 4, the least squares method is recommended. In Section 5, the error analysis is obtained. In Section 6, the accuracy of the method is verified by two numerical examples. In Section 7, the conclusion is given.


## 2. Preliminaries

Definition 1.
(see [31]) (BPFs).
BPFs have been considered and used to solve various equations by numerous scholars, such as [8, 12, 13, 28, 30]. BPFs are denoted as(2)ψit=1,i−1h≤t<ih,0,otherwise,where t∈0,T,i=1,…,m, and h=T/m.
The basic properties of BPFs are described as follows:(1)
Disjointness:(3)ψitψjt=δijψit,i,j=1,2,…,m,whereδij is the Kronecker delta.(2)
Orthogonality:(4)∫0Tψitψjtdt=hδij.(3)
Completeness: Parseval’s identity is true for everyf∈L20,T and ft=∑i=1∞fiψit,(5)∫0Tf2tdt=∑i=1∞fi2ψit,where fi=1/h∫0Tftψitdt.
BPFs are expressed as a vector(6)Ψt=ψ1t,ψ2t,…,ψmtT.
It can be concluded from the above description that(7)ΨmtΨmTt=ψ1t0⋯00ψ2t⋯0⋮⋮⋱000⋯ψmtm×m,ΨmTtΨmt=1,ΨmtΨmTtFm=DFmΨmt.where Fm=f1,f2,…,fmT and DFm=diagFm.
Any functionXt∈L20,T can be approximately expressed as(8)Xt≃Xmt=∑i=1mciψit=ΨmTtCm,Cm=c1,c2,…,cmT,where Xmt is a linear combination of BPFs.
LetSu,t∈L20,T1×0,T2; it can be extended as(9)Su,t=Ψm1TuSΨm2t=Ψm2TtSTΨm1u,where S=sijm1×m2,(10)sij=1h1h2∫0T1∫0T2Su,tψiuϕjtdudt,and h1=T1/m1,h2=T2/m2. In order to facilitate, we put m1=m2=m in the following sections.


## 3. Stochastic Integration Operational Matrix

The section gives the relevant lemmas.Lemma 1.
(see [8, 13]). Supposing Ψmt is given in (6), we obtain(11)∫0tΨmudu≃QΨmt,where(12)Q=h2122⋯2012⋯2001⋯2⋮⋮⋮⋱⋮000⋯1m×m.Therefore, each functionXt∈L20,T can be approximately expressed as(13)∫0tXudu≃CmTQΨmt.Lemma 2 (see [8, 13]).
SupposingΨmt is given in (6), we obtain(14)∫0tΨmudBu≃QBΨmt,where QB can be expressed as(15)QB=β1α1α1⋯α10β2α2⋯α200β3⋯α3⋮⋮⋮⋱⋮000⋯βmm×m,where αi=Bih−Bi−1h,i=1,2,…,m−1; βj=Bih/2−Bi−1h/2,j=1,2,…,m.Therefore, each functionXt∈L20,T can be approximately expressed as(16)∫0tXudBu≃CmTQBΨmt.


## 4. Method Description

In this section, the multidimensional linear SIVIEs are transformed into linear algebraic equations by the least square method and BPFs.The linear operatorL is taken into account,(17)LXt≔Xt−∫0tS˜u,tXudu−∑k=1q∫0tS^ku,tXudBku.IfXt is the exact solution to (1), then the residual norm disappears:(18)LX−fL2=0.Forε>0, Xεt is the approximate solution, so that the residual norm of Xεt is less than ε,(19)Xεt−fL2<ε.An approximate solution can be obtained by (19), and Xm is a linear combination by BPFs,(20)Xm=∑i=1mciψit,where ci is the unknown coefficient.The minimization problem is studied below,(21)minc1,…,cmLXm−fL2.Minimizing (21) to get a set of values of c^1,…,c^m, then ∑i=1mc^iψit is an approximate solution to (1), and we’ll prove it.(22)Eminc1,…,cmLXm−fL2=EL∑i=1mc^iψi−fL2⟶0,where m⟶∞ (or h⟶0).To obtain the minimum value of the approximate solution in (21), we have to take the partial derivative in ci, and let(23)∂∂ci∫0T∑j=1mcjLψjt−ft2dt=0,where i=1,2,…,m.Then,(24)∫0T∑j=1mcjLψjtLψitdt=∫0TftLψitdt,that is,(25)∑j=1mcjLψit,Lψjt=Lψit,ft,where Lψit,Lψjt:=∫0TLψjtLψitdt.The following is expressed as a matrix:(26)c=G−1a,where c=c1,c2,…,cmT∈ℝm,a=a1,a2,…,amT∈ℝm, and G=gij∈ℝm×m such that(27)ai=Lψit,ft,gij=Lψit,Lψjt,i,j=1,2,…,m.For solving equation (1), Xt,ft,S˜u,t, and S^u,t can be approximated by BPFs in the following forms:(28)Xt≃Xmt=CmTΨmt=ΨmTtCm,(29)ft≃f0t=FmTΨmt=ΨmTtFm,(30)S˜u,t=ΨmTuS1Ψmt=ΨmTtS1TΨmu,(31)S^u,t=ΨmTuS2Ψmt=ΨmTtS2TΨmu,where Cm and Fm are BPF coefficient vectors; S1 and S2 are BPF coefficient matrices by the equation of (10).Next, for obtainingai and gij, by (17) and (28)–(31), we have(32)Lψit=ψit−∫0tS˜u,tψiudu−∑k=1q∫0tS^ku,tψiudBku=IiΨmt−∫0tIiΨmuΨmTuS1Ψmtdu−∑k=1q∫0tIiΨmuΨmTuS2ΨmtdBku=IiΨmt−Ii∫0tΨmuΨmTuduS1Ψmt−Ii∑k=1q∫0tΨmuΨmTudBkuS2Ψmt,where Ii is the ith row of an m×m identity matrix.LetS1i and S2i be the ith row of the matrix S1 and S2, respectively, Ri and RBi be the ith row of the integral operation matrix Q and the stochastic integral matrix QB, respectively, and DS1i and DS2i be a diagonal matrix with S1i and S2i as the diagonal elements, respectively. By (11) and (14), we get(33)∫0tΨmuΨmuTduS1Ψmt,=R1ΨmtS11ΨmtR2ΨmtS12Ψmt⋮RmΨmtS1mΨmt=R1DS11R2DS12⋮RmDS1mΨmt=S˜Ψmt,where(34)S˜=h2S˜112S˜122S˜13⋯2S˜1m0S˜222S˜23⋯2S˜2m00S˜33⋯2S˜3m⋮⋮⋮⋱⋮000⋯S˜mmm×m,∫0tΨmuΨmuTdBuS2Ψmt=RB1ΨmtS21ΨmtRB2ΨmtS22Ψmt⋮RBmΨmtS2mΨmt=RB1DS21RB2DS22⋮RBmDS2mΨmt=S^Ψmt,where(35)S^k=S^11kβ1S^12kα1S^13kα1⋯S^1mkα10S^22kβ2S^23kα2⋯S^2mkα200S^33kβ3⋯S^3mkα3⋮⋮⋮⋱⋮000⋯S^mmkβmm×m.According to Lemmas1 and 2,(36)Lψit=IiΨmt−IiS˜Ψmt−Ii∑k=1qS^kΨmt=Ii−IiS˜−Ii∑k=1qS^kΨmt.Then,(37)ai=Lψit,ft=∫0TLψitftdt=∫0TIi−IiS˜−Ii∑k=1qS^kΨmtFmTΨmtdt=∫0TIi−IiS˜−Ii∑k=1qS^kFmΨmTtΨmtdt=Ii−IiS˜−Ii∑k=1qS^kFm∫0TΨmTtΨmtdt=Ii−IiS˜−Ii∑k=1qS^kFmT,gij=Lψit,Lψjt=∫0TLψitLψjtdt=∫0TIi−IiS˜−Ii∑k=1qS^kΨmtIj−IjS˜−Ij∑k=1qS^kΨmtdt=∫0TIi−IiS˜−Ii∑k=1qS^kIj−IjS˜−Ij∑k=1qS^kTΨmTtΨmtdt=Ii−IiS˜−Ii∑k=1qS^kIj−IjS˜−Ij∑k=1qS^kT∫0TΨmTtΨmtdt=Ii−IiS˜−Ii∑k=1qS^kIj−IjS˜−Ij∑k=1qS^kTT.


## 5. Error Analysis

We will give the error analysis in this section.Lemma 3 (see [29]) (continuous module).
f with respect to ϱ in 0,T is the definition of continuous modulus ωf,ϱ(38)ωf,ϱ=supfx−fyx,y∈0,T,x−y≤ϱ.As shown in the reference [32], if and only if limρ⟶0ωf,ϱ=0, on 0,T, ft is uniformly continuous.Lemma 4 (see [29]).
Supposef∈C0,T, Qh=∑i=1meiψit, where ei=fti−1+ti/2, ti=ih, h=T/m, i=0,1,…,m. Then,(39)f−Qh∞≤ωf,h.Theorem 1.
SupposeXt be the exact solution of (1) and X^mt=∑i=1meiψit, where ei=Xti−1+ti/2, S˜u,t, and S^ku,t are deterministic functions, S˜∞≤M,S^k∞≤M, where M is a positive constant. Then, we get the following conclusions, when h⟶0:(40)iX^mt−Xt∞,E=Esupt∈0,TX^mt−Xt⟶0,iiEminc1,…,cmXmt−ft−∫0tS˜u,tXmudu−∑k=1q∫0tS^ku,tXmudBkuL22⟶0.Proof.
(i) By Lemma4, we have(41)‖X^mt−Xt‖∞,E=Esupt∈0,TX^mt−Xt≤ωX,h⟶0.
(ii)(42)Eminc1,…,cmXmt−ft−∫0tS˜u,tXmudu−∑k=1q∫0tS^ku,tXmudBkuL22≤EX^mt−ft−∫0tS˜u,tX^mudu−∑k=1q∫0tS^u,tX^mudBkuL22≤T2X^mt−ft−∫0tS˜u,tX^mudu−∑k=1q∫0tS^u,tX^mudBku∞,E2=T2X^mt−Xt+∫0tS˜u,tXu−X^mudu+∑k=1q∫0tS^u,tXu−X^mudBku∞,E2≤3T2X^mt−Xt∞,E2+∫0tS˜u,tXu−X^mudu∞,E2+∑k=1q∫0tS^u,tXu−X^mudBku∞,E2≤3T2X^mt−Xt∞,E2+3T2M2Xt−X^mt∞,E2+3T2∑k=1q∫0tS^u,tXu−X^mudBku∞,E2.
According to isometry property and Doob’s inequality,(43)∑k=1q∫0tS^u,tX−X^mudBku∞,E2=Esup0≤τ≤t∑k=1q∫0τS^u,tXu−X^mudBku2≤4E∑k=1q∫0tS^u,tXu−X^mudBku2≤4Eq∑k=1q∫0tS^u,tXu−X^mudBku2=4Eq2∫0tS^u,t2Xu−X^mu2du≤4TM2q2Xt−X^mt∞,E2.
By (41)–(43), we get(44)Eminc1,…,cmXmt−ft−∫0tS˜u,tXmudu−∑k=1q∫0tS^ku,tXmudBkuL22≤3T2+3T4M2+12T3M2q2X^mt−Xt∞,E2≤3T2+3T4M2+12T3M2q2ωX,h⟶0,and the proof of the theorem is completed.


## 6. Numerical Examples

The following two examples show the availability and accuracy of the method mentioned in Section4. The algorithms were implemented by MATLAB 2016a.Example 1.
The linear SIVIE is studied as follows:(45)Xt=X0+∫0trXudu+∑k=1q∫0takXudBku,u,t∈0,1,where Xt=X0er−1/2∑k=1qak2t+∑k=1qakBkt, 0≤t<1, and Bku are q-dimensional Brownian motions. When X0=1/200,r=1/20,a1=1/50,a2=2/50,a3=4/50, and a4=9/50, the numerical results are shown in Tables 1 and 2. The trajectories are shown in Figures 1 and 2, where n=100 is the number of iterations.
We can see from Figures1 and 2 and Tables 1 and 2 that the approximate solution and the exact solution are very close. This method is more accurate than the method in [30], and this method is feasible and effective for solving the q-dimensional SIVIEs.Table 1
Whenm=24, the table shows error means Em, error standard deviations Es, and confidence intervals for different time t.
tEmEs95% confidence interval for error meanLowerUpper1/161.36075252×10−46.80376264×10−52.7215050×10−62.6942900×10−43/161.13092940×10−45.65464702×10−52.2618588×10−62.2392402×10−45/165.33574672×10−52.66787336×10−51.0671493×10−61.0564778×10−47/162.67128325×10−51.33564162×10−55.3425665×10−75.2891408×10−59/162.53074880×10−61.26537440×10−65.0614976×10−85.0108826×10−6Table 2
Whenm=25,n=100, this table shows error means Em, error standard deviations Es, and confidence intervals for different time t.
tEmEs95% confidence interval for error meanLowerUpper1/321.47378183×10−47.36890915×10−52.9475636×10−62.9180880×10−47/329.11708782×10−54.55854391×10−51.8234175×10−61.8051833×10−413/325.04388689×10−52.52194344×10−51.0087773×10−69.9868960×10−519/329.47544810×10−64.73772405×10−61.8950896×10−71.8761387×10−525/324.47443350×10−52.23721675×10−58.9488670×10−78.8593783×10−5Figure 1
m=24, the approximate solution and exact solution for Example 1.Figure 2
m=25, the approximate solution and exact solution for Example 1.Example 2.
The linear SIVIE is studied as follows:(46)Xt=X0+∫0te−t−uXudu+∑k=1q∫0take−t−uXudBku,u,t∈0,1,where Bku are q-dimensional Brownian motions. When X0=1,a1=1/50,a2=2/50,a3=4/50, and a4=9/50, the trajectories are shown in Figures 3 and 4.Figure 3
m=24, the approximate solution and mean solution for Example 2.Figure 4
m=25, the approximate solution and mean solution for Example 2.


## 7. Conclusion

In this paper, a numerical method based on the least squares method and BPFs is proposed to solve multidimensional linear SIVIEs different from that of Maleknejad et al., who only used BPFs to solve equations [30]. The least squares method is a mathematical optimization technique that seeks the best function match of data by minimizing errors and quickly obtaining unknown coefficients. Combining it with BPFs can make the error of the numerical solution of multidimensional linear SIVIEs smaller. In Section 6, the numerical simulations are carried out and compared with the literature of [30]; this paper concludes that this method is better than applying only the block pulse functions directly.



---

*Source: 6662604-2021-02-24.xml*