1 Introduction
As we all know, fractional differential equations (FDEs) can be viewed as a generalization of ordinary differential equations, which can be used to model complex physical phenomena and processes with nonlocal properties. In recent twenty years, with the continuous development of fractional differential equations, they have very important applications in the fields of mechanics, electrical engineering, electromagnetic wave, population system
[31, 34, 14, 30, 37] and so on. Among them, multiterm FDEs are known as an important tool in describing viscoelastic damping materials, modelling nonlinear wave phenomenon in plasma and simulating anomalous diffusive process, such as anomalous relaxation in magnetic resonance imaging signal magnitude, mechanical models of oxygen delivery through capillaries [21, 29, 32, 1], etc. Thus, there are lots of excellent researches focusing on the analytical methods and numerical methods for multiterm FDEs, see [22, 8, 10, 5, 13] and the references therein.Meanwhile, with scientific research constantly deepening, researchers find that there are always some noise disturbances that could not be ignored in real life, and in order to better describe the phenomena and process influenced by these noisy factors, researchers pay more attention to stochastic differential equations (SDEs). As we all know, many classical SDEs play an extremely important role in describing physical phenomena. For instance, Langevin equations are used to describe Brownian motion and explain Einstein relations [19] and the stochastic NavierStokes equation are often used to simulate various problems in fluid motion[35]. Nowadays, in addition to the field of physics, SDEs are widely used in option pricing, population growth and many other fields [11, 20, 28]. Especially, stochastic fractional differential equations (SFDEs) appeal more scholars’ attention and many studies have been carried out, such as the random motion of harmonically trapped charged particles in a constant external magnetic field, the relationship between fluctuationdissipation theorem and physical behavior, the noise driving in some financial models, epidemiological research [25, 23, 27, 36] and so on. Obviously, it is quite difficult or even rarely impossible to obtain the exact solutions of SFDEs, thus there has been a growing interest to construct numerical methods for these equations. Until now, many numerical methods have been developed for SFDEs. For example, Kamrani in [18] investigated a numerical solution of SFDEs driven by additive noise and proved the convergence of the proposed method. Jin et al. in [17] studied the stochastic timefractional diffusion problem driven by fractionally integrated Gaussian noise and developed a numerical scheme by employing the Galerkin finite element method in space, GrnwaldLetnikov formula in time and projection for the noise. Doan et al. in [12]
constructed an EulerMaruyama (EM) method for a kind of SFDEs driven by a multiplicative white noise with the Caputo fractional order
. Zhou et al. in [40] used the finite difference method to solve the stochastic fractional nonlinear wave equation and presented the performance of numerical solution and property of energy under effect of two different types of noise, i.e., additive noise and multiplicative noise. Additionally, there also exist many other good works of SFDEs, see [4, 39, 2, 3], for examples.Motivated by the above numerical methods for SFDEs, we in this paper will construct and analyze a modified EM method for the following multiterm RiemannLiouville SFDEs,
(1.1) 
where

are the RiemannLiouville fractional derivatives with ;
In addition to the purpose of numerically solving the multiterm RiemannLiouville SFDEs (1.1), we hope to construct more efficient numerical methods to avoid the huge computational cost caused by the nonlocal property of RiemannLiouville fractional integral operators. Recently, lots of efficient methods based on the sumofexponentials (SOEs) technique proposed in [16] by Jiang et al. are presented to overcome this difficulty for FDEs, see [41, 6, 15, 33], for examples. There are also some articles about applying the SOEs technique to efficiently solve the SFDEs. For example, Dai et al in [9] proposed an EM method and its fast implementation based on the SOEs approximation for Lévydriven stochastic Volterra integral equations with doubly singular kernels. Ma et al. [24] used the SOEs approximation to simulate the rough volatility and developed a fast twostep iteration algorithm. We in [38] also presented a fast EM method for a class of nonlinear SFDEs by the SOEs technique. Herein, based on the SOEs approximation and the proposed EM method, a fast EM method for (1.1) is also introduced and analyzed in this paper.
The rest of this paper is organized as follows. In Section 2, some notations and preliminaries are given. In Section 3, the modified EM method is derived and the strong convergence of this method is proved. Section 4 aims to construct a fast EM method based on SOEs technique and present the corresponding numerical theoretical results. Two numerical examples are shown to examine the theoretical results and illustrate the effectiveness of the proposed two methods. Finally, a brief conclusion is given.
2 Notations and preliminaries
In this section, we first introduce some useful notations that will be used throughout this paper. Let denote the expectation corresponding to . Define as the space of all measurable, mean square integral functions with standard Euclidean norm . For two real numbers and , we denote and . Moreover, the capital letter will be used to represent a positive constant whose value may change when it appears in different places. The following two definitions are taken from [26].
Definition 2.1
The RiemannLiouville fractional integral of order of function is defined by
with , where is the Euler gamma function.
Definition 2.2
The RiemannLiouville derivative of order of function can be written as
In the rest of this section, to guarantee the existence and uniqueness of the solution of equation (1.1), we impose three important hypotheses.
Assumption 2.1
(Hlder continuity) There exists a positive constant such that for all and , and satisfy the condition:
Assumption 2.2
(Lipschitz continuity) There exists a positive constant such that for all and , and satisfy the inequality:
Assumption 2.3
(Linear growth condition) There exists a positive constant such that for all and all , and satisfy the condition:
3 The modified EM method for multiterm RiemannLiouville SFDEs
3.1 The modified EM method
In order to get the numerical method of the RiemannLiouville SFDEs (1.1), we first transform equation (1.1) into the following stochastic Volterra integral equation (SVIE) by taking the RiemannLiouville fractional integral operator on both sides of equation (1.1), that is
(3.1) 
where the property of RiemannLiouville fractional integral and RiemannLiouville fractional derivative is used, i.e.
In another way, for each , when the above equality holds for , a adapted process is called a solution of equation (1.1) on the interval with the initial condition .
For every integer , the modified EM method for equation (3.1) can be presented as
(3.2)  
where denote the grid points with step size , denote the increments of Brownian motion. To facilitate its convergence analysis, we introduce a continuoustime version
(3.3)  
where for .
3.2 Strong convergence of the modified EM method
To prove the strong convergence of the modified EM method (3.3), we first list some necessary lemmas.
Lemma 3.1
If , then for any , the following inequality holds
Proof. Define with . Then we have . According to the property of power function, it holds that . Thus it implies that together with .
Lemma 3.2
Suppose , then there exists a positive constant independent of such that for any ,
where have been defined in the above subsection.
Proof. For , it is easily derived that
The following inequality can easily be obtained by the property of norm.
Lemma 3.3
For all , , the following inequality holds
To carry out the proof of the strong convergence of the modified EM method (3.3
), the bounded estimate for the numerical solutions is given.
Theorem 3.1
Proof. From equation (3.3) and Lemma 3.3, it is deduced that
Then using the Hlder inequality and It isometry as well as Lemma 3.3 implies that
Applying the linear growth condition of Assumption 2.3 and the property of power function with the fact , it can be deduced that
Arranging the above inequality, it becomes
There exists a positive constant such that
Taking the supremum on both sides of the above inequality and according to the Gronwall’s inequality, we can arrive at
Similarly, we can prove that . This proof is completed.
Lemma 3.4
Suppose Assumption 2.3 holds, then for all , there exists a positive constant independent of such that for all ,
Proof. For arbitrary , it follows from equation (3.3) that
Applying Lemma 3.3, we obtain
By Hlder inequality and It isometry, the above inequality turns into
This together with Lemma 3.1, Assumption 2.3 and Theorem 3.1 as well as the fact implies that
Noticing , thus we can get
This proof is completed.
Theorem 3.2
Proof. From SVIE (3.1) and the modified EM method (3.3), by adding some intermediate items and using Lemma 3.3, Hlder inequality and It isometry, we derive that
Based on Assumptions 2.1 and 2.2, Lemmas 3.2 and 3.4 and Theorem 3.1 as well as the fact , it is easy to obtain that
Arranging the above inequality yields
where . Then using the Gronwall’s inequality and arbitrariness of , we arrive at
The proof is completed.
4 The fast EM method for multiterm RiemannLiouville SFDE
4.1 The fast EM method
To construct the fast EM method, it is necessary to introduce the sumofexponentials (SOEs) approximation.
Lemma 4.1
[16] For a given , let denote tolerance error, denote cutoff time restriction and denote final time, there are a positive integer and positive constants and , such that for any
where .