Galerkin Approximation for Stochastic Volterra Integral Equations with Doubly Singular Kernels

Abstract

This paper is concerned with the more general nonlinear stochastic Volterra integral equations with doubly singular kernels, whose singular points include both $s=t$ and $s=0$. We propose a Galerkin approximate scheme to solve the equation numerically, and we obtain the strong convergence rate for the Galerkin method in the mean square sense. The rate is $min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}$ (where ${\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2$ are positive numbers satisfying $0<{\alpha }_1+{\beta }_1<1$, $0<{\alpha }_2+{\beta }_2<\frac{1}{2}$), which improves the results of some numerical schemes for the stochastic Volterra integral equations with regular or weakly singular kernels. Moreover, numerical examples are given to support the theoretical result and explain the priority of the Galerkin method.

1. Introduction

Stochastic Volterra integral equations (SVIEs) were firstly introduced by Berger and Mizel [1,2] and can be regarded as generalizations of stochastic differential equations or deterministic Volterra integral equations. Due to the advantages of the SVIEs in describing memory, heritability and ubiquitous noise perturbations, these kinds of equations have extensive applications in many fields (e.g., finance [3], control science [4] and mathematical biology [5]) and have spurred great research enthusiasm. An area of particular interest in the study of SVIEs has been numerical analysis because the analytic solution to SVIEs is rarely known and the numerical approximations provide a powerful tool for understanding the behavior of the solution. Up to now, most numerical methods are developed to deal with SVIEs with regular kernels. For example, Tudor [6] first proposed a one-step numerical approximation for It$\widehat{o}$–Volterra equations and obtained a basic convergence theorem. Wen and Zhang [7,8] studied the rectangular method. Methods based on the operational matrix were introduced in [9,10,11]. Xiao et al. [12] introduced a split-step collocation method for SVIEs. The Euler–Maruyama (EM) methods were discussed in [13,14,15]. Recently, Conte et al. [16] introduced improved stochastic $\theta$-methods, having better stability properties, for the numerical integration of stochastic Volterra integral equations.

However, stochastic Volterra integral equations with singular kernels arise while dealing with some problems in stochastic partial differential equations [17,18,19], such as certain heat conduction problems with mixed boundary conditions, and in the analysis of the fractional Brownian motion [20]. It is therefore of great value and significance to consider SVIEs involving the kernels with singularities. Wang [21] established the existence and uniqueness theorem for stochastic Volterra equations with singular kernels. Zhang [22] investigated the convergence of the EM method for this kind of equation. In particular, Li et al. [23] recently proposed a $\theta$-Euler–Maruyama scheme and a Milstein scheme to solve SVIEs with weakly singular kernels numerically. It is observed that the singular SVIEs have been relatively less studied. It is probably because the singularity of the integrand kernel brings us more difficulties. In this case, the key It$\widehat{o}$ formula, which is a powerful and necessary tool in the study of stochastic differential equations (SDEs), is not available and we have to seek other tools and techniques.

In this paper, we will consider a class of more general stochastic Volterra integral equations with doubly singular kernels of the following form

$$x\left(t\right)=x_0+{\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}f\left(x\left(s\right)\right)ds+{\int }_0^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)dB\left(s\right),0\le t\le T,$$

(1)

where $T>0$ and ${\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2$ are positive numbers satisfying $0<{\alpha }_1+{\beta }_1<1$, $0<{\alpha }_2+{\beta }_2<\frac{1}{2}$; $f,g:[0,T]\to \mathbb{R}$ are given measurable functions; $B\left(t\right)$ is a standard Brownian motion defined on a complete probability space $(\mathrm{\Omega },\mathcal{F},\mathbb{P})$ and adapted to the filtration ${\left\{{\mathcal{F}}_t\right\}}_{0\le t\le T}$; and $x_0$ is ${\mathcal{F}}_0$-measurable random variable such that $\mathbb{E}|x_0{|}^2<+\infty$.

As we can see, SVIE (1) involves both diagonal and boundary singular kernels, that is, the kernel functions are singular when $s=t$ or $s=0$, which are more general than the case considered in [23,24,25], whose singular point only includes $s=t$. The existence and uniqueness theorem of the true solution as well as the strong convergence of the Euler–Maruyama method for Equation (1) have been developed in [26]. Li et al. [27] obtained the asymptotic separation for Equation (1).

Our aim here is to establish the Galerkin approximations for Equation (1). This method is high-order accurate, easy-to-handle complicated geometries and side conditions, highly parallelizable, nonlinear stable and has the ability to capture discontinuities without spurious oscillations [28]. Hence, it has been extensively used for stochastic partial differential equations [29,30,31] and stochastic fractional differential equations [32,33], while it has not been analyzed to singular SVIEs.

The rest of this paper is organized as follows.

In Section 2, we present some preliminaries. Section 3 provides the main result and Section 4 covers numerical experiments to illustrate the efficiency of the numerical method. We conclude this work in Section 5.

2. Preliminaries

2.1. Lemmas and Assumption

Let us first introduce two important lemmas which will play critical roles in proving our main result.

Denote $l_{n+1}=l_n·\frac{\mathrm{\Gamma }\left(n\right(1-\alpha -\beta )+1-\beta )}{\mathrm{\Gamma }\left(n\right(1-\alpha -\beta )+2-\alpha -\beta )}$, $n=0,1,2,\cdots$ with $l_0=1$, where $\mathrm{\Gamma }\left(z\right)={\int }_0^{\infty }{x}^{z-1}exp(-x)dx$, $z>0$ is the Gamma function. Define

$$\begin{array}{c}\hfill E_{1-\alpha ,1-\beta }\left(s\right)=\sum _{n=0}^{\infty }{l}_n{s}^{n(1-\alpha -\beta )}.\end{array}$$

(2)

This function $E_{1-\alpha ,1-\beta }$ is closely related to the Mittag-Leffler function and has the following asymptotic property.

Lemma 1

([27]). The function $E_{1-\alpha ,1-\beta }$ defined as in (2) satisfies

$$\begin{array}{c}\hfill E_{1-\alpha ,1-\beta }\left(s\right)=O\left(s^{\frac{1}{2}(\frac{1-\alpha -\beta }{1-\alpha }-1+\beta )}exp\left(\frac{1-\alpha }{1-\alpha -\beta }·s^{\frac{1-\alpha -\beta }{1-\alpha }}\right)\right)\mathrm{as}s\to \infty .\end{array}$$

Lemma 2

([27]). Assume that $\alpha >0,\beta >0,\alpha +\beta <1,a\ge 0,b\ge 0$ and u is a nonnegative function satisfying that $t^{-\beta }u\left(t\right)$ is locally integrable on ${\mathbb{R}}_{+}$.

(1) If u satisfies

$$u\left(t\right)\le a+b{\int }_0^t{(t-s)}^{-\alpha }{s}^{-\beta }u\left(s\right)ds\forall t\in {\mathbb{R}}_{+},$$

then

$$u\left(t\right)\le a{E}_{1-\alpha ,1-\beta }\left({\left(b\mathrm{\Gamma }(1-\alpha )\right)}^{1/(1-\alpha -\beta )}t\right)\forall t\in {\mathbb{R}}_{+}.$$

(2) If u satisfies

$$u\left(t\right)\ge a+b{\int }_0^t{(t-s)}^{-\alpha }{s}^{-\beta }u\left(s\right)ds\forall t\in {\mathbb{R}}_{+},$$

then

$$u\left(t\right)\ge a{E}_{1-\alpha ,1-\beta }\left({\left(b\mathrm{\Gamma }(1-\alpha )\right)}^{1/(1-\alpha -\beta )}t\right)\forall t\in {\mathbb{R}}_{+}.$$

Throughout this paper, we impose the following conditions on the coefficients f and g.

Assumption 1.

Assume that there exists a constant $L>0$ such that for any $x,y\in \mathbb{R}$, we have

$$\left|f\right(x)-f(y\left)\right|+\left|g\right(x)-g(y\left)\right|\le L|x-y|,\left|f\right(x\left)\right|+\left|g\right(x\left)\right|\le L(1+|x\left|\right).$$

(3)

2.2. A Numerical Approximation of White Noise for Stochastic Integral

Next, we introduce the approximation of the noise. Following the method discussed in [29], let

$$0=t_0<t_1<\cdots <t_N=T$$

with $t_i=i\mathrm{\Delta }t$ and $\mathrm{\Delta }t=\frac{T}{N}$ for $i=0,...,N$ being a partition of $[0,T]$. Define

$$\frac{d\widehat{B}\left(t\right)}{dt}=\frac{1}{\sqrt{\mathrm{\Delta }t}}\sum _{i=1}^N{\eta }_i{\xi }_i\left(t\right),$$

(4)

where

$${\xi }_i\left(t\right)=\left\{\begin{array}{c}\hfill 1,t_i\le t<t_{i+1},\\ \hfill 0,otherwise,\end{array}\right.$$

and ${\eta }_i\sim N(0,1)$ is defined by

$${\eta }_i=\frac{1}{\sqrt{\mathrm{\Delta }t}}{\int }_{t_i}^{t_{i+1}}dB\left(t\right),i=0,\cdots ,N-1.$$

We substitute $dB\left(t\right)$ with $d\widehat{B}\left(t\right)$ in Equation (1) to obtain the following equation

$$\widehat{x}\left(t\right)=x_0+{\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}f\left(\widehat{x}\left(s\right)\right)ds+{\int }_0^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(\widehat{x}\left(s\right)\right)d\widehat{B}\left(s\right),0\le t\le T.$$

(5)

Let $\widehat{x}\left(t\right)$ be the solution of the approximate SVIEs

$$\frac{d\widehat{x}\left(t\right)}{dt}=\sum _{i=1}^N{(t-t_i)}^{-{\alpha }_1}{t}_i^{-{\beta }_1}f\left(\widehat{x}\left(t\right)\right)+\sum _{i=1}^N{(t-t_i)}^{-{\alpha }_2}{t}_i^{-{\beta }_2}g\left(\widehat{x}\left(t\right)\right)\frac{1}{\sqrt{\mathrm{\Delta }t}}{\eta }_i{\xi }_i\left(t\right).$$

(6)

2.3. The Discontinuous Galerkin Method

Because (6) is a system of ODEs, we apply the standard discontinuous Galerkin method for deterministic ODEs [28]. Multiplying (6) by v, integrating over $I_n$ and using integration by parts, we obtain the discontinuous Galerkin formulation

$$\begin{array}{cc}\hfill & {\int }_{I_n}{\widehat{x}}_i{v}^{\prime }dt-{\widehat{x}}_i\left(t_{n+1}\right)v\left(t_{n+1}\right)+{\widehat{x}}_i\left(t_n\right)v\left(t_n\right)\hfill \\ \hfill =& -{\int }_{I_n}{(t-t_i)}^{-{\alpha }_1}{t}_i^{-{\beta }_1}{\widehat{f}}_i\left(t\right)vdt-\frac{{\eta }_i}{\sqrt{\mathrm{\Delta }t}}{\int }_{I_n}{(t-t_i)}^{-{\alpha }_2}{t}_i^{-{\beta }_2}{\widehat{g}}_i\left(t\right)vdt,i=1,2,\dots ,d.\hfill \end{array}$$

We define the piecewise polynomial space $V_h^p=\left\{v:{\left.v\right|}_{I_n}\in P^p\left(I_n\right),n=0,1,\dots ,N-1\right\}$ as the space of polynomials of degree at most p in each interval $I_n$, where $P^p\left(I_n\right)$ is the set of all polynomials of degree less or equal than p on $I_n$. Because polynomials in $V_h^p$ are allowed to have discontinuities across element boundaries, we use $v\left(t_n^{\pm }\right)={lim}_{s\to 0^{\pm }}v\left(t_n+s\right)$ to denote the left limit and the right limit of v at $t_n$. Next, we approximate each ${\widehat{x}}_i\left(t\right)$ by a piecewise polynomial $x_{i,h}\left(t\right)\in V_h^p$. The discrete discontinuous Galerkin scheme consists of finding $x_{i,h}\in V_h^p$ such that: $\forall v\in V_h^p$ and $n=0,\dots ,N-1,i=1,2,\dots ,d$.

$$\begin{array}{cc}\hfill & -{\int }_{I_n}{v}^{\prime }{x}_{i,h}dt+x_{i,h}\left(t_{n+1}^{-}\right)v\left(t_{n+1}^{-}\right)-x_{i,h}\left(t_n^{-}\right)v\left(t_n^{+}\right)\hfill \\ \hfill =& {\int }_{I_n}{(t-t_i)}^{-{\alpha }_1}{t}_i^{-{\beta }_1}{\widehat{f}}_i\left(t\right)vdt+\frac{{\eta }_i}{\sqrt{\mathrm{\Delta }t}}{\int }_{I_n}{(t-t_i)}^{-{\alpha }_2}{t}_i^{-{\beta }_2}{\widehat{g}}_i\left(t\right)vdt.\hfill \end{array}$$

where we used the classical upwind numerical flux. We will refer to this discontinuous Galerkin scheme as the stochastic discontinuous Galerkin scheme.

3. Main Result

Theorem 1

([26]). Under Assumption 1, there exists a unique strong solution $x\left(t\right)$ to Equation (1). Moreover, there exists a constant $C>0$ such that

$$\underset{0\le t\le T}{sup}\mathbb{E}\left[{\left|x\left(t\right)\right|}^2\right]\le C,forgivenT>0.$$

The main result of this paper is the following theorem.

Theorem 2.

Under Assumption 1, there exists a constant $C>0$ such that

$$\mathbb{E}|x\left(t\right)-\widehat{x}{\left(t\right)|}^2\le C{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}},$$

(7)

where, for simplicity of notation, C is a positive constant which depends on ${\alpha }_1,{\alpha }_2,{\beta }_1,{\beta }_2,T,L$ and may change from line to line during the rest of the paper, while its specific form is of unimportance.

In order to prove Theorem 2, we need the following Lemma.

Lemma 3.

Let $x\left(t\right)$ be the solution of (1), then for any $0\le t_1<t_2\le T$, we have

$$\mathbb{E}|x\left(t_2\right)-x\left(t_1\right){|}^2\le C{(t_2-t_1)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}.$$

(8)

Proof. 

Observe that

$$\begin{array}{cc}\hfill x\left(t_2\right)-x\left(t_1\right)& ={\int }_0^{t_1}[{(t_2-s)}^{-{\alpha }_1}-{(t_1-s)}^{-{\alpha }_1}]s^{-{\beta }_1}f\left(x\left(s\right)\right)ds\hfill \\ \hfill & +{\int }_{t_1}^{t_2}{(t_2-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}f\left(x\left(s\right)\right)ds\hfill \\ \hfill & +{\int }_0^{t_1}[{(t_2-s)}^{-{\alpha }_2}-{(t_1-s)}^{-{\alpha }_2}]s^{-{\beta }_2}g\left(x\left(s\right)\right)dB\left(s\right)\hfill \\ \hfill & +{\int }_{t_1}^{t_2}{(t_2-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)dB\left(s\right).\hfill \end{array}$$

It follows that

$$\begin{array}{cc}\hfill \mathbb{E}|x\left(t_2\right)-x\left(t_1\right){|}^2& \le 4\mathbb{E}{\left|{\int }_0^{t_1}[{(t_2-s)}^{-{\alpha }_1}-{(t_1-s)}^{-{\alpha }_1}]s^{-{\beta }_1}f\left(x\left(s\right)\right)ds\right|}^2\hfill \\ \hfill & +4\mathbb{E}{\left|{\int }_{t_1}^{t_2}{(t_2-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}f\left(x\left(s\right)\right)ds\right|}^2\hfill \\ \hfill & +4\mathbb{E}|{\int }_0^{t_1}[{(t_2-s)}^{-{\alpha }_2}-{(t_1-s)}^{-{\alpha }_2}]s^{-{\beta }_2}g\left(x\left(s\right)\right)dB\left(s\right){|}^2\hfill \\ \hfill & +4\mathbb{E}|{\int }_{t_1}^{t_2}{(t_2-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)dB\left(s\right){|}^2\hfill \\ \hfill & :=I_1+I_2+I_3+I_4.\hfill \end{array}$$

On the one hand, Cauchy–Schwarz inequality and basic calculus imply that

$$\begin{array}{cc}\hfill I_1& \le C{(t_2-t_1)}^{1-{\alpha }_1-{\beta }_1}{\int }_0^{t_1}|{(t_2-s)}^{-{\alpha }_1}-{(t_1-s)}^{-{\alpha }_1}|s^{-{\beta }_1}\mathbb{E}{\left|f\left(x\left(s\right)\right)\right|}^{2}ds\hfill \\ \hfill & \le C{(t_2-t_1)}^{1-{\alpha }_1-{\beta }_1}{\int }_0^{t_1}|{(t_2-s)}^{-{\alpha }_1}-{(t_1-s)}^{-{\alpha }_1}|s^{-{\beta }_1}(1+\underset{0\le s\le T}{sup}\mathbb{E}|x\left(s\right){|}^2)ds\hfill \\ \hfill & \le C{(t_2-t_1)}^{2-2{\alpha }_1-2{\beta }_1}\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill I_2& \le 4\mathbb{E}{\int }_{t_1}^{t_2}{(t_2-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}ds{\int }_{t_1}^{t_2}{(t_2-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}{\left|f\left(x\left(s\right)\right)\right|}^{2}ds\hfill \\ \hfill & \le C{\left[{\int }_{t_1}^{t_2}{(t_2-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}ds\right]}^2\hfill \\ \hfill & \le C{\left[{\int }_{t_1}^{t_2}{(t_2-s)}^{-{\alpha }_1}{(s-t_1)}^{-{\beta }_1}ds\right]}^2\hfill \\ \hfill & \le C{(t_2-t_1)}^{2-2{\alpha }_1-2{\beta }_1}.\hfill \end{array}$$

On the other hand, using It$\widehat{\mathrm{o}}$ isometry, we obtain

$$\begin{array}{cc}\hfill I_3& =4\mathbb{E}{\int }_0^{t_1}{\left|[{(t_1-s)}^{-{\alpha }_2}-{(t_2-s)}^{-{\alpha }_2}]s^{-{\beta }_2}g\left(x\left(s\right)\right)\right|}^{2}ds\hfill \\ \hfill & =4\mathbb{E}{\int }_0^{t_1}[{(t_1-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}-2{(t_1-s)}^{-{\alpha }_2}{(t_2-s)}^{-{\alpha }_2}{s}^{-2{\beta }_2}+{(t_2-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}]{\left|g\left(x\left(s\right)\right)\right|}^{2}ds\hfill \\ \hfill & \le \mathbb{E}{\int }_0^{t_1}[{(t_1-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}-{(t_2-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}]{\left|g\left(x\left(s\right)\right)\right|}^{2}ds\hfill \\ \hfill & \le C{(t_2-t_1)}^{1-2{\alpha }_2-2{\beta }_2}\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill I_4& =4\mathbb{E}{\int }_{t_1}^{t_2}{(t_2-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}{\left|g\left(x\left(s\right)\right)\right|}^{2}ds\hfill \\ \hfill & \le C{\int }_{t_1}^{t_2}{(t_2-s)}^{-2{\alpha }_2}{(s-t_1)}^{-2{\beta }_2}ds\hfill \\ \hfill & \le C{(t_2-t_1)}^{1-2{\alpha }_2-2{\beta }_2}.\hfill \end{array}$$

The proof is completed. □

Proof of Theorem 2. 

Let $e\left(t\right)=x\left(t\right)-\widehat{x}\left(t\right)$, then we have

$$\begin{array}{cc}\hfill e\left(t\right)& ={\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}[f\left(x\left(s\right)\right)-f\left(\widehat{x}\left(s\right)\right)]ds\hfill \\ \hfill & +\left[{\int }_0^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)dB\left(s\right)+{\int }_0^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(\widehat{x}\left(s\right)\right)d\widehat{B}\left(s\right)\right]\hfill \\ \hfill & :=J_1+J_2.\hfill \end{array}$$

For the term $J_1$, it follows from Assumption 1 that

$$\begin{array}{cc}\hfill \mathbb{E}|J_1{|}^2& \le L{\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}ds{\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}\mathbb{E}{|x\left(s\right)-\widehat{x}\left(s\right)|}^{2}ds\hfill \\ \hfill & \le C{\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}\mathbb{E}{\left|e\left(s\right)\right|}^{2}ds.\hfill \end{array}$$

Next, let us estimate the term $J_2$. Observe that

$$\begin{array}{cc}\hfill J_2& ={\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)[dB\left(s\right)-d\widehat{B}\left(s\right)]\hfill \\ \hfill & +{\int }_0^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}[g\left(x\left(s\right)\right)-g\left(\widehat{x}\left(s\right)\right)]d\widehat{B}\left(s\right)\hfill \\ \hfill & :=J_3+J_4.\hfill \end{array}$$

In order to estimate $J_3$, we assume there exists an integer $N_0\le N$ such that $t_{N_0}\le t<t_{N_0+1}$, then

$$\begin{array}{cc}\hfill J_3& =\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)[dB\left(s\right)-d\widehat{B}\left(s\right)]\hfill \\ \hfill & +{\int }_{t_{N_0}}^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)[dB\left(s\right)-d\widehat{B}\left(s\right)]\hfill \\ \hfill & :=J_{3,1}+J_{3,2}.\hfill \end{array}$$

According to (4), we obtain

$$\begin{array}{cc}\hfill J_{3,1}& =\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)-\frac{1}{\mathrm{\Delta }t}{\int }_{t_i}^{t_{i+1}}{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}g\left(x\left(\tau \right)\right)d\tau ]dB\left(s\right)\hfill \\ \hfill & =\sum _{i=0}^{N_0-1}\frac{1}{\mathrm{\Delta }t}{\int }_{t_i}^{t_{i+1}}[{\int }_{t_i}^{t_{i+1}}[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)-{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}g\left(x\left(\tau \right)\right)]d\tau ]dB\left(s\right).\hfill \end{array}$$

Based on the fact

$$\begin{array}{cc}\hfill & {\int }_{t_i}^{t_{i+1}}[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}g\left(x\left(s\right)\right)-{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}g\left(x\left(\tau \right)\right)]d\tau \hfill \\ \hfill =& {\int }_{t_i}^{t_{i+1}}\left[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}(g\left(x\left(s\right)\right)-g\left(x\left(\tau \right)\right))\right]d\tau \hfill \\ \hfill & +{\int }_{t_i}^{t_{i+1}}g\left(x\left(\tau \right)\right)[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}-{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}]d\tau \hfill \\ \hfill \le & L{\int }_{t_i}^{t_{i+1}}{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}|x\left(s\right)-x\left(\tau \right)|d\tau \hfill \\ \hfill & +L{\int }_{t_i}^{t_{i+1}}(1+x\left(\tau \right))[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}-{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}]d\tau ,\hfill \end{array}$$

one obtains

$$\begin{array}{cc}\hfill \mathbb{E}|J_{3,1}{|}^2& \le \frac{C}{\mathrm{\Delta }{t}^2}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}\mathbb{E}{\left|{\int }_{t_i}^{t_{i+1}}{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}|x\left(s\right)-x\left(\tau \right)|d\tau \right|}^{2}ds\hfill \\ \hfill & +\frac{C}{\mathrm{\Delta }{t}^2}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}\mathbb{E}{\left|{\int }_{t_i}^{t_{i+1}}(1+x\left(\tau \right))[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}-{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}]d\tau \right|}^{2}ds\hfill \\ \hfill & \le \frac{C}{\mathrm{\Delta }t}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{\int }_{t_i}^{t_{i+1}}{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}{|s-\tau |}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}d\tau ds\hfill \\ \hfill & +\frac{C}{\mathrm{\Delta }{t}^2}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{\left|{\int }_{t_i}^{t_{i+1}}[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}-{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}]d\tau \right|}^{2}ds\hfill \\ \hfill & :=J_{3,1,1}+J_{3,1,2}.\hfill \end{array}$$

On the one hand, it is easy to see that

$$\begin{array}{cc}\hfill J_{3,1,1}& \le \frac{C}{\mathrm{\Delta }t}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}{(s-t_i)}^{1+min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}ds\hfill \\ \hfill & \le C{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}ds\hfill \\ \hfill & \le C{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}{\int }_0^{t_{N_0}}{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}ds.\hfill \end{array}$$

On the other hand,

$$\begin{array}{cc}\hfill J_{3,1,2}& \le \frac{C}{\mathrm{\Delta }{t}^2}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{\left[{\int }_{t_i}^s[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}-{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}]d\tau \right]}^{2}ds\hfill \\ \hfill & =\frac{C}{\mathrm{\Delta }t}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{\left[{\int }_{t_i}^s[{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}-{(t-t_i)}^{-{\alpha }_2}{s}^{-{\beta }_2}]d\tau \right]}^{2}ds\hfill \\ \hfill & \le \frac{C}{\mathrm{\Delta }t}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{\left[{\int }_{t_i}^s{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}d\tau \right]}^{2}ds\hfill \\ \hfill & \le \frac{C}{\mathrm{\Delta }t}\sum _{i=0}^{N_0-1}{\int }_{t_i}^{t_{i+1}}{(t-s)}^{-2{\alpha }_2}{(s-t_i)}^{2-2{\beta }_2}ds\hfill \\ \hfill & \le C{\left(\mathrm{\Delta }t\right)}^{1-2{\beta }_2}{\int }_0^{N_0}{(t-s)}^{-2{\alpha }_2}ds.\hfill \end{array}$$

Consequently,

$$\begin{array}{cc}\hfill \mathbb{E}|J_{3,1}{|}^2& \le C{\left(\mathrm{\Delta }t\right)}^{1-2{\beta }_2}{\int }_0^{t_{N_0}}{(t-s)}^{-2{\alpha }_2}ds\hfill \\ \hfill & +C{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}{\int }_0^{t_{N_0}}{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}ds.\hfill \end{array}$$

(9)

Similarly, we can obtain

$$\begin{array}{cc}\hfill \mathbb{E}|J_{3,2}{|}^2& \le C{\left(\mathrm{\Delta }t\right)}^{1-2{\beta }_2}{\int }_{t_{N_0}}^t{(t-s)}^{-2{\alpha }_2}ds\hfill \\ \hfill & +C{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}{\int }_{t_{N_0}}^t{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}ds.\hfill \end{array}$$

(10)

Combining (9) with (10), we obtain

$$\begin{array}{cc}\hfill \mathbb{E}|J_3{|}^2& \le C{\left(\mathrm{\Delta }t\right)}^{1-2{\beta }_2}{\int }_0^t{(t-s)}^{-2{\alpha }_2}ds\hfill \\ \hfill & +C{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}{\int }_0^t{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}ds\hfill \\ \hfill & \le C\left[{\left(\mathrm{\Delta }t\right)}^{1-2{\beta }_2}+{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}\right]\hfill \\ \hfill & \le C{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}.\hfill \end{array}$$

Next, we estimate the term $J_4$. For simplicity, let $t=t_{N_0+1}$, then

$$\begin{array}{cc}\hfill \mathbb{E}|J_4{|}^2& =\mathbb{E}{\left|{\int }_0^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}[g\left(x\left(s\right)\right)-g\left(\widehat{x}\left(s\right)\right)]d\widehat{B}\left(s\right)\right|}^2\hfill \\ \hfill & \le \mathbb{E}{\left|\sum _{i=0}^{N_0}\frac{1}{\mathrm{\Delta }t}{\int }_{t_i}^{t_{i+1}}\left[{\int }_{t_i}^{t_{i+1}}{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}(g\left(x\left(\tau \right)\right)-g\left(\widehat{x}\left(\tau \right)\right))d\tau \right]dB\left(s\right)\right|}^2\hfill \\ \hfill & \le \sum _{i=0}^{N_0}\frac{1}{\mathrm{\Delta }{t}^2}\mathbb{E}{\int }_{t_i}^{t_{i+1}}{\left|{\int }_{t_i}^{t_{i+1}}{(t-\tau )}^{-{\alpha }_2}{\tau }^{-{\beta }_2}(g\left(x\left(\tau \right)\right)-g\left(\widehat{x}\left(\tau \right)\right))d\tau \right|}^{2}ds\hfill \\ \hfill & \le C\sum _{i=0}^{N_0}\frac{1}{\mathrm{\Delta }t}{\int }_{t_i}^{t_{i+1}}{\int }_{t_i}^{t_{i+1}}{(t-\tau )}^{-2{\alpha }_2}{\tau }^{-2{\beta }_2}\mathbb{E}{e}^2\left(\tau \right)d\tau ds\hfill \\ \hfill & \le C{\int }_0^t{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}\mathbb{E}{e}^2\left(s\right)ds.\hfill \end{array}$$

Combining this with the above estimates, we can derive that

$$\begin{array}{cc}\hfill {\mathbb{E}\left|e\left(t\right)\right|}^2& \le C[{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}\hfill \\ \hfill & +{\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}{\mathbb{E}\left|e\left(s\right)\right|}^{2}ds+{\int }_0^t{(t-s)}^{-2{\alpha }_2}{s}^{-2{\beta }_2}\mathbb{E}{\left|e\left(s\right)\right|}^{2}ds].\hfill \end{array}$$

A simple application of the generalized Gronwall’s inequality in Lemma 2 yields

$$\begin{array}{cc}\hfill {\mathbb{E}\left|e\left(t\right)\right|}^2& \le C{\left(\mathrm{\Delta }t\right)}^{min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}}\hfill \\ \hfill & ·max\left\{E_{1-{\alpha }_1,1-{\beta }_1}{\left[\left(C\mathrm{\Gamma }(1-{\alpha }_1)\right)\right]}^{1/(1-{\alpha }_1-{\beta }_1)}T,E_{1-2{\alpha }_2,1-2{\beta }_2}{\left[\left(C\mathrm{\Gamma }(1-2{\alpha }_2)\right)\right]}^{1/(1-2{\alpha }_2-2{\beta }_2)}T\right\}.\hfill \end{array}$$

Finally, Lemma 1 helps us derive the required assertion. We complete the proof. □

4. Numerical Simulation

In this section, we will verify the numerical solution of the stochastic Volterra integral equations with doubly singular kernels and the strong convergence rate of the Galerkin method.

Example 1.

Let us consider the following stochastic Volterra integral equations with doubly singular kernels

$$x\left(t\right)=2+{\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}sin\left(2x\left(s\right)\right)ds+{\int }_0^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}sin\left(x\left(s\right)\right)dB\left(s\right),0\le t\le T.$$

(11)

When the terminal time $t_N=T=1$, the positive arguments ${\alpha }_i$ and ${\beta }_i(i=1,2)$ take the following four cases:

• Case I: ${\alpha }_1=0.4,{\beta }_1=0.2,{\alpha }_2=0.1,{\beta }_2=0.1;$
• Case II: ${\alpha }_1=0.4,{\beta }_1=0.2,{\alpha }_2=0.2,{\beta }_2=0.1;$
• Case III: ${\alpha }_1=0.4,{\beta }_1=0.4,{\alpha }_2=0.1,{\beta }_2=0.1;$
• Case IV: ${\alpha }_1=0.4,{\beta }_1=0.2,{\alpha }_2=0.1,{\beta }_2=0.2.$

Obviously, functions f and g satisfy the Assumption 1. By Theorem 2, the convergence rate of the Galerkin scheme is $min\left\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\right\}$, we obtain the convergence rate of

Table 1. Convergence rate for four cases of different values.

${\mathit{\alpha }}_1$	${\mathit{\beta }}_1$	${\mathit{\alpha }}_2$	${\mathit{\beta }}_2$	Convergence Rate
0.4	0.2	0.1	0.1	0.3
0.4	0.2	0.2	0.1	0.2
0.4	0.4	0.1	0.1	0.2
0.4	0.2	0.1	0.2	0.2

. In order to verify the numerical solution of Equation (11) and the strong convergence rate of the Galerkin method, we adopt five different step sizes $h=2^{-4},2^{-5},2^{-6},2^{-7},2^{-8}$ on the same Brownian motion to compute the four cases, respectively. The numerical solution by the Galerkin method with small step size $h^{\ast }=2^{-14}$ is used to represent the unknown true solution. Figure 1 draws the calculated error results, which are consistent with Theorems 1 and 2. Compared with the classic Euler–Maruyama method or Milstein method, the Galerkin method is more accurate in algorithm estimation by space discretization, for example, when ${\beta }_1={\beta }_2=0$, Equation (1) became SVIE with weakly singular kernels, revealing the convergence rate of the Euler–Maruyama method is $min\{1-{\alpha }_1,\frac{1}{2}-{\alpha }_2\}$ in [23].

Example 2.

We consider the following example

$$X\left(t\right)=1+{\int }_0^t{(t-s)}^{-{\alpha }_1}{s}^{-{\beta }_1}sin\left(x\left(s\right)\right)ds+\frac{1}{2}{\int }_0^t{(t-s)}^{-{\alpha }_2}{s}^{-{\beta }_2}(\sin \left(x\left(s\right)\right)+2)d{W}_s,t\in [0,1].$$

(12)

Due to the appearance of the singularity in the above stochastic integral, it is difficult for us to illustrate the convergence rate of the Galerkin scheme numerically. We set ${\alpha }_1=0.3,{\alpha }_1=0.5,{\beta }_1=0.2,{\beta }_2=0.3$. We regard the numerical solution yielded by small step size $h^{\ast }=2^{-13}$ as the ’exact’ solution. Moreover, the corresponding numerical solutions are generated by four different step sizes $h=2^{-4},2^{-5},2^{-6}$ and $2^{-7}$, respectively. The mean square errors of the Galerkin scheme are calculated at the terminal time $t_N=T=1$ by

$$e={\left(\frac{1}{M}\sum _{i=1}^M{\left|x\left(t\right)-\widehat{x}\left(t\right)\right|}^2\right)}^{\frac{1}{2}}$$

where the expectation is approximated by averaging over $M=1000$ Brownian sample paths. The mean square errors are plotted in Figure 2 in a loglog scale. In these plots, the reference lines and error lines are parallel to each other, revealing the convergence rate of the Galerkin scheme is $min\left\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\right\}$. Therefore, the convergence rate is $0.4$.

5. Conclusions

In this work, we investigated the Galerkin approximation for a class of nonlinear stochastic Volterra integral equations with doubly singular kernels. These kinds of equations are more general, for example, when ${\beta }_1={\beta }_2=0$, Equation (1) becomes the SVIE with weakly singular kernels studied in [23]. Because the It$\widehat{o}$ formula is not available, the classical proof techniques can no longer be used. With the help of new tools such as classical fractional calculus, the Mittag-Leffler-type function and generalized Gronwall inequalities with singular kernels, we obtained the strong convergence rate $min\{2-2({\alpha }_1+{\beta }_1),1-2({\alpha }_2+{\beta }_2)\}$ for the Galerkin method, which improves the corresponding result of the $\theta$-Euler–Maruyama scheme in [23]. In forthcoming works, we will consider other numerical methods for SVIEs with doubly singular kernels and verify whether the order of convergence is optimal. Moreover, we will try to analyze the stability of the numerical methods.