Monte-Carlo Galerkin Approximation of Fractional Stochastic Integro-Differential Equation

Abstract

A stochastic differential equation, SDE, describes the dynamics of a stochastic process defined on a space-time continuum. This paper reformulates the fractional stochastic integro-differential equation as a SDE. Existence and uniqueness of the solution to this equation is discussed. A numerical method for solving SDEs based on the Monte-Carlo Galerkin method is presented.

1. Introduction

Recently, many applications in numerous fields, such as viscoelastic materials, signal processing, controlling, quantum mechanics, meteorology, finance, and life science have been remodeled in terms of fractional calculus where derivatives and integrals of fractional order is introduced and so differential equation of fractional order are involved in these models, see [1–5]. Generally, a fractional integro-differential equation with Caputo's definition of fractional derivative takes the form 𝑢^𝛼(𝑡)=𝑓(𝑡,𝑢)+^𝑡_𝑡0ℎ(𝑡,𝑠,𝑢)𝑑𝑠,𝑢𝑡₀=𝑢₀,𝑡₀≤𝑡≤𝑇,0≤𝛼≤1. (1.1) Also, in recent years, development of adequate statistical techniques for stochastic systems has constantly been in the focus of scientific attention because of its outstanding importance for a number of physical applications such as turbulence, heterogeneous flows and materials, random noises, and so forth, see [6–8]. A stochastic integro-differential equation (SIDE) takes the form 𝑢^𝛼(𝑡;𝑤)=𝑓(𝑡,𝑢;𝑤)+^𝑡_𝑡0ℎ(𝑡,𝑠,𝑢;𝑤)𝑑𝑊_𝑠,𝑢𝑡₀=𝑢₀,𝑡₀≤𝑡≤𝑇,0≤𝛼≤1, (1.2) where 𝑊 is a Brownian motion.

In this paper we study the FSIDE: 𝑢^𝛼(𝑡;𝑤)+𝐴(𝑡)𝑢(𝑡;𝑤)=^𝑡_𝑡0𝐵(𝑡,𝑠;𝑤)𝑢(𝑠;𝑤)𝑑𝑊_𝑠+𝑓(𝑡;𝑤),𝑡>𝑡₀ (1.3) with initial condition 𝑢(𝑡₀;𝑤)=𝑢₀ and 0≤𝛼≤1.

The FSDE is a generalization of the fractional Fokker-Planck equation which describes the random walk of a particle, see [9]. The aim of this paper is threefold. First we rewrite the above equation as stochastic differential equation 𝑑𝑢(𝑡)=𝑎(𝑡,𝑢;𝑤)𝑑𝑡+𝑏(𝑡,𝑢;𝑤)𝑑𝑊_𝑠. (1.4) Secondly we prove the existence and the uniqueness of the solution of this equation. Thirdly we present a numerical method using finite element method and the Monte-Carlo method.

2. Preliminaries

In this section we collect a few known results to which we refer frequently in the sequel.

Let Ω={𝜔₁,𝜔₂,…} be the collection of all outcomes, 𝜔_𝑖,𝑖=1,2,… of an experiment, and let 𝐴⊂Ω. Define the measure 𝑃∶𝐴→[0,1] with 𝑃(Ω)=1. The triplet (Ω,𝐴,𝑃) is called a probability space.

Definition 2.1. A collection {𝑋(𝑡,⋅)∶𝑡≥0} of random variables 𝑋∶Ω×𝑇→ℛ, 𝑇=[0,∞] is called a stochastic process. For each point 𝜔∈Ω, the mapping 𝑡→𝑋(𝑡,𝜔) is a realization, sample path, or trajectory of the stochastic process.

Definition 2.2. A real-valued stochastic process that depends continuously on 𝑡∈[0,𝑇]𝑊={𝑊(𝑡),𝑡≥0},𝑊∶𝑇×Ω→ℛ is called a Brownian motion (BM) or Wiener process if (1)𝑊(0)=0 a.s., (2)𝑊(𝑡)−𝑊(𝑠) is 𝑁(0,𝑡−𝑠), the Gaussian distribution with mean 0 and variance 𝑡−𝑠, for all 𝑡>𝑠>0, (3)for all times 0<𝑡₁<𝑡₂<⋯<𝑡_𝑛, the random variables 𝑊(𝑡₁),𝑊(𝑡₂)−𝑊(𝑡₁),…,𝑊(𝑡_𝑛)−𝑊(𝑡_𝑛−1) are independents. Also, the expected values of 𝑊(𝑡) and 𝑊²(𝑡) are given by 𝐸(𝑊(𝑡))=0 and 𝐸(𝑊²(𝑡))=𝑡 for each time 𝑡>0.

Definition 2.3. (i) Suppose that {𝐴_𝑡≥0} is an increasing family of 𝜎-algebras of 𝐴 such that 𝑊_𝑡 is 𝐴_𝑡-measurable 𝐸(𝑊(𝑡)∣𝐴₀)=0 and 𝐸(𝑊(𝑡)−𝑊(𝑠)∣𝐴_𝑠)=0 with probability one, w.p.1, for all 0≤𝑠≤𝑡≤𝑇.
(ii) Let ℒ be the 𝜎-algebra of Lebesgue subsets of ℛ. Let ℒ²_𝑇 be a class of functions 𝑓∶[0,𝑇]×Ω→ℛ satisfying the following: (a)𝑓 is jointly ℒ×𝐴-measurable; (b)∫^𝑇₀𝐸((𝑓(𝑡,⋅))²)𝑑𝑡<∞; (c)𝐸((𝑓(𝑡,⋅))²) for each 0≤𝑠≤𝑡≤𝑇; (d)𝑓(𝑡,⋅) is 𝐴_𝑡 measurable for each 0≤𝑡≤𝑇.

Definition 2.4 (Convergence). If 𝑢(ℎ) is the exact value of a random variable and 𝑢_𝑛(ℎ) is its approximate value, then we have (1)strong convergence when 𝐸𝑢(ℎ)−𝑢_𝑛(ℎ)⟶0,𝑛⟶∞, (2.1)(2)weak convergence when 𝐸(𝑢(ℎ))−𝐸𝑢_𝑛(ℎ)⟶0,𝑛⟶∞, (2.2)(3)mean square convergence when 𝐸𝑢(ℎ)−𝑢_𝑛(ℎ)²⟶0,𝑛⟶∞. (2.3)

Definition 2.5. Let 𝑦(𝑥)∈𝐿₁([𝑎,𝑏]), we define the fractional integral of order 𝛼, 1≥𝛼>0, of the function 𝑦(𝑥) over the interval [𝑎,𝑏] as 𝐽^𝛼_𝑎+𝑦(𝑥)=1Γ(𝛼)^𝑥_𝑎𝑦(𝑡)(𝑥−𝑡)^1−𝛼𝐽𝑑𝑡,𝑥>𝑎,_𝛼𝑏−𝑦(𝑥)=1Γ(𝛼)^𝑏_𝑥𝑦(𝑡)(𝑡−𝑥)^1−𝛼𝑑𝑡,𝑥<𝑏. (2.4) For completion, we define 𝐽⁰=𝐼 (identity operator); that is, we mean 𝐽⁰𝑦(𝑥)=𝑦(𝑥). Furthermore, by 𝐽^𝛼_𝑎+𝑦(𝑥) we mean the limit (if it exists) of 𝐽^𝛼_𝑐 as 𝑐→𝑎⁺. This definition is according to Riemann-Liouville definition of fractional integral of arbitrary order 𝛼. For simplicity, we define the fractional integral by the equation 𝐽^𝛼𝑦(𝑥)=1Γ(𝛼)^𝑥_𝑎𝑦(𝑡)(𝑥−𝑡)^1−𝛼𝑑𝑡,𝑥>𝑎,0<𝛼≤1, (2.5) where we dropped 𝑎⁺. The folllowing equation 𝐷^𝛼𝑦(𝑥)=𝑑𝐽𝑑𝑥^1−𝛼𝑦(𝑥) (2.6) denotes the fractional derivative of 𝑦(𝑥) of order 𝛼. For numerical computations, we usually use another definition of the fractional derivative, due to Caputo [2], given by ^𝐶𝐷^𝛼𝑦(𝑥)=𝐷^𝛼𝑦(𝑡)−𝑦(𝑎)(𝑥). (2.7)

3. Formulation of the Problem

In this paper, we present a numerical method to solve the following FSIDE: ^𝐶𝐷^𝛼_𝑡𝑢(𝑡;𝑤)−𝐴(𝑡)𝑢(𝑡;𝑤)=^𝑡_𝑡0𝐵(𝑡,𝑠;𝑤)𝑢(𝑠;𝑤)𝑑𝑊_𝑠+𝑓(𝑡;𝑤),𝑡>𝑡₀ (3.1) with initial condition 𝑢(𝑡₀;𝑤)=𝑢_𝑡0.

For simplicity of notations, we drop the variable 𝑤, and so this equation takes the form ^𝐶𝐷^𝛼𝑢(𝑡)+𝐴(𝑡)𝑢(𝑡)=^𝑡_𝑡0𝐵(𝑡,𝑠)𝑢(𝑠)𝑑𝑊_𝑠+𝑓(𝑡),𝑡>𝑡₀,𝑢𝑡₀=𝑢_𝑡0, (3.2) where ^𝐶𝐷^𝛼 is the Caputo-fractional derivative operator of order 𝛼 with 0<𝛼≤1,𝑡∈[𝑡₀,𝑇]. The integral term of the right hand side is an It𝑜̂ integral, 𝑊 is a BM, and 𝑑𝑊/𝑑𝑡 is a white noise, the derivative of a BM (the derivative in the sense of distribution). The functions 𝐴(𝑡),𝐵(𝑡,𝑠;𝑤),𝑓(𝑡;𝑤),𝑢_𝑡0(𝑤), and 𝑢(𝑡;𝑤) satisfy the following conditions: C1 (measurability): 𝐴=𝐴(𝑡),𝐵=𝐵(𝑠,𝑡) and 𝑢(𝑡) are ℒ²-measurable in [0,1]×ℛ; C2: there exists constants 𝐾₁,𝐾₂,𝐾₃>0 such that 𝐴(𝑡)(𝑡−𝑠)^1−𝛼<𝐾₁,𝑓(𝑡)(𝑡−𝑠)^1−𝛼<𝐾₂,_𝑡𝑠𝐵(𝑧,𝑠)𝑑𝑧(𝑡−𝑧)^1−𝛼<𝐾₃ (3.3) for all 𝑡₀≤𝑠≤𝑡≤𝑇 and 0<𝛼≤1; C3 (initial value): 𝑢_𝑡0 is 𝐴_𝑡0-measurable with 𝐸(𝑢_𝑡0)<∞.

Let 𝑘=max{𝑘₁,𝑘₂,𝑘₃}/Γ(𝛼).

Equation (3.2) may be written as 𝑢(𝑡)=𝑢_𝑡0+^𝐶𝐷^−𝛼𝑓(𝑡)+^𝐶𝐷^−𝛼𝐴(𝑡)𝑢(𝑡)+^𝐶𝐷^−𝛼𝑡_𝑡0𝐵(𝑡,𝑠)𝑢(𝑠)𝑑𝑊(𝑠). (3.4) Using the definition of the operator ^𝐶𝐷^−𝛼, (2.7), we obtain that 𝑢(𝑡)=𝑢_𝑡0+1Γ(𝛼)^𝑡_𝑡0𝑓(𝑧)𝑑𝑧(𝑡−𝑧)^1−𝛼+1Γ(𝛼)^𝑡_𝑡0𝐴(𝑧)𝑢(𝑧)𝑑𝑧(𝑡−𝑧)^1−𝛼+1Γ(𝛼)^𝑡_𝑡0^𝑧₀𝐵(𝑧,𝑠)𝑢(𝑠)(𝑡−𝑧)^1−𝛼𝑑𝑊_𝑠𝑑𝑧. (3.5) This is a Fredholm integral equation that we are going to solve instead of solving (3.2).

Without loss of generality, we will let 𝑡₀=0 and 𝑇=1 throughout; that is, we assume that the time parameter set of the processes considered is [0,1].

Lemma 3.1. The FSIDE (3.5) can be written as a stochastic integral equation 𝑢(𝑡)=𝑢_𝑡0+^𝑡_𝑡0𝑎𝑡,𝑠,𝑢_𝑠𝑑𝑠+^𝑡_𝑡0𝑏𝑡,𝑠,𝑢_𝑠𝑑𝑊_𝑠, (3.6) where the drift is given by 𝑎𝑡,𝑠,𝑢_𝑠=1Γ(𝛼)𝐴(𝑠)𝑢_𝑠+𝑓(𝑠)(𝑡−𝑠)^1−𝛼0≤𝑠<𝑡≤1,0else, (3.7) and the diffusion is 𝑏𝑡,𝑠,𝑢_𝑠=𝑢_𝑠Γ(𝛼)^𝑡_𝑧=𝑠𝐵(𝑧,𝑠)𝑑𝑧(𝑡−𝑧)^1−𝛼0≤𝑠<𝑡≤1,0else. (3.8)

The proof of this lemma and the next lemma are easy to see.

Lemma 3.2. The conditions C1–C3 imply the conditions A1–A4 where A1 (measurability): the drift and the diffusion are jointly ℒ²-measurable in [0,1]×ℛ; A2 (lipschitz condition): there is a constant 𝐾>0 such that |𝑎(𝑡,𝑠,𝑢)−𝑎(𝑡,𝑠,𝑣)|≤|𝑢−𝑣|,|𝑏(𝑡,𝑠,𝑢)−𝑏(𝑡,𝑠,𝑣)|≤|𝑢−𝑣|, (3.9) for 0≤𝑠≤𝑡≤1,𝑢,𝑣∈𝑅; A3 (linear growth bound): there is a constant 𝐾>0 such that |𝑎(𝑡,𝑠,𝑢)|²≤𝐾²1+|𝑢|²,|𝑏(𝑡,𝑠,𝑢)|²≤𝐾²|𝑢|²; (3.10)A4 (Initial value): 𝑢₀ is 𝐴₀-measurable with 𝐸(𝑢₀)<∞.
With the help of these two lemmas, we can establish the existence and the uniqueness theorem to the (3.5); for proof see Theorem 4.5.3, Kloeden and Platen [10].

Theorem 3.3. Under the assumptions C1–C3 (or A1–A4) the stochastic integral equation (3.6), and so (3.2), has a pathways unique strong solution 𝑢_𝑡 on [0,1].
The existence of a unique solution of the SDE (3.2) ensures the existence of the integrals on the right hand side of (3.5) at each point in the domain of the definition.

4. The Monte Carlo Galerkin Finite Element Method

It is known that, introducing a finite element method (FEM) that approximates a solution of differential equation (DE), we first need to obtain a weak formulation in the standard sense of DE and FEM which is not possible with the presence of the white noise. In our method, we do not require to approximate the white noise using FEM, instead we follow the approach in Allen et al. [11] who have suggested a smoother approximation for the white noise process when computing the approximate solutions of stochastic differential equations.

They have suggested the following approximation for the one-dimensional white noise process 𝑊̇(𝑡),𝑡∈[0,1].

Consider the uniform time discretization of the interval [0,1]ℑ_𝜌=𝐼_𝑛∶𝐼_𝑛=𝑡_𝑛−1,𝑡_𝑛,𝑛=[1∶𝑁],𝑡_𝑛−𝑡_𝑛−1=𝜌,𝑡₀=0,𝑡_𝑁=1. (4.1) Then the following approximation is defined for the white noise process 𝑊̇(𝑡) on this discretization, 𝑑𝑊_𝑡≈𝑑𝑊𝑑𝑡^𝑡=𝑑𝑡^𝑁_𝑛=1𝜂_𝑛𝜒_𝑛(𝑡), (4.2) where the coefficients are random variables defined by 𝜂_𝑛=1𝜌_ℑ𝜌𝜒_𝑛(𝑡)(𝑡)1𝑑𝑊∈𝑁0,𝜌,𝜒_𝑛(𝑡)=1,𝑡_𝑛−1≤𝑡≤𝑡_𝑛,0,else. (4.3) As a direct result of this approximation of the white noise, it is easily to see that 𝐸(∫_ℑ𝜌𝑔[𝑑𝑊_𝑡−𝑑𝑊_𝑡])=0=𝐸((∫_ℑ𝜌𝑔[𝑑𝑊_𝑡−𝑑𝑊_𝑡])²) for any bounded function 𝑔.

Now 𝑑𝑊_𝑡 can be substituted for 𝑑𝑊_𝑡 to obtain the following smoothed version of the SDE (3.5): 𝑢(𝑡)=𝑢₀+1Γ(𝛼)^𝑡₀𝑓(𝑧)𝑑𝑧(𝑡−𝑧)^1−𝛼+1Γ(𝛼)^𝑡₀𝐴(𝑧)𝑢(𝑧)𝑑𝑧(𝑡−𝑧)^1−𝛼+1Γ(𝛼)^𝑡₀^𝑧₀𝐵(𝑧,𝑠)𝑢(𝑠)(𝑡−𝑧)^1−𝛼𝑑𝑊_𝑠𝑑𝑧. (4.4) The solution of this equation, 𝑢(𝑡), is smoother than 𝑢(𝑡) and therefore standard numerical procedures can be applied to compute its approximate solution. Once the approximate solution, 𝑢, is obtained for 𝑀 realizations of such approximate solution, the Monte Carlo method then uses these approximations to compute corresponding sample averages of these 𝑀 realizations.

Now we show that 𝑢(𝑡) is a good approximation of 𝑢(𝑡). To show this the following lemma is required.

Lemma 4.1. Given a nonrandom function 𝐻(𝑡) such that |𝐻(𝑡)−𝐻(𝜏)|≤𝜆|𝑡−𝜏|^𝛿,0≤𝑡,𝜏≤1,0≤𝛿≤1, (4.5) where 𝜆>0 is a constant, then 𝐸^𝑧₀𝐻(𝑡)𝑑𝑊_𝑡−^𝑧₀𝐻(𝜏)𝑑𝑊_𝜏²≤𝜆²𝜌^2𝛿 (4.6) with 𝜌≥|𝑡−𝜏| for 0≤𝑡,𝜏,𝑧≤1.

Proof. Let ℑ^∗_𝜌=𝐼^∗_𝑛∶𝐼^∗_𝑛=𝑡^∗_𝑛−1,𝑡^∗_𝑛,𝑛=1∶𝑁^∗,𝑡^∗_𝑛−𝑡^∗_𝑛−1=𝜌^∗,𝑡₀=0,𝑡^∗_𝑁=𝑧,𝜌^∗≤𝜌, (4.7) be a partition of the interval [0,𝑧]. We have 𝐸^𝑧₀𝐻(𝑡)𝑑𝑊_𝑡−^𝑧₀𝐻(𝜏)𝑑𝑊_𝜏²=𝐸^𝑁∗_𝑛=1^𝑡∗𝑛_{𝑡∗𝑛−1}𝐻(𝑡)−1𝜌^𝑡∗𝑛_𝑡^∗_𝑛−1𝐻(𝜏)𝑑𝜏𝑑𝑊_𝑡²=^𝑁∗_𝑛=1^𝑡∗𝑛_{𝑡∗𝑛−1}𝐻(𝑡)−1𝜌^𝑡∗𝑛_𝑡^∗_𝑛−1𝐻(𝜏)𝑑𝜏²=𝑑𝑡_{𝑁∗𝑛=1𝑡∗𝑛𝑡∗𝑛−1}1𝜌^𝑡∗𝑛_𝑡^∗_𝑛−1{𝐻(𝑡)−𝐻(𝜏)}𝑑𝜏²≤𝜆𝑑𝑡₂𝜌^2𝑁∗_𝑛=1^𝑡∗𝑛_𝑡^∗_𝑛−11𝜌^𝑡∗𝑛_𝑡^∗_𝑛−1|𝑡−𝜏|^𝛿𝑑𝜏²𝑑𝑡≤𝜆²𝜌^∗2𝛿≤𝜆²𝜌^2𝛿. (4.8) With the help of this lemma, we can show that 𝑢(𝑡)→𝑢(𝑡) as 𝜌→0.

Theorem 4.2. 𝐸¹₀(𝑢(𝑡)−𝑢(𝑡))²𝜌𝑑𝑡≤𝐶²1−3𝜆², (4.9) where 𝜆²=1Γ²(𝛼)¹₀^𝑡₀1(𝑡−𝑧)^2−2𝛼𝐴²(𝑧)+^𝑧₀𝐵²(𝑧,𝑠)𝑑𝑠𝑑𝑧𝑑𝑡. (4.10)

Proof. Let 𝜖(𝑡)=𝑢(𝑡)−𝑢(𝑡), then (3.5) and (4.4) lead to 𝜖(𝑡)=1Γ(𝛼)^𝑡₀𝐴(𝑧)𝜖(𝑧)𝑑𝑧(𝑡−𝑧)^1−𝛼+1Γ(𝛼)^𝑡₀^𝑧₀𝐵(𝑧,𝑠)𝜖(𝑠)(𝑡−𝑧)^1−𝛼𝑑𝑊_𝑠+1𝑑𝑧Γ(𝛼)^𝑡₀^𝑧₀𝐵(𝑧,𝑠)𝑢(𝑠)(𝑡−𝑧)^1−𝛼𝑑𝑊−𝑑𝑊_𝑠𝑑𝑧. (4.11) Applying the inequality (𝑎+𝑏+𝑐)²≤3(𝑎²+𝑏²+𝑐²) and the Hölder's inequality to this equation, we obtain ¹₀𝜖²(𝑡)𝑑𝑡≤3𝜆²¹₀𝜖²(𝑡)𝑑𝑡+𝜙, (4.12) where 𝜆² is given by (4.10) and 1𝜙=Γ²(𝛼)¹₀^𝑡₀^𝑧₀𝐵(𝑧,𝑠)𝑢(𝑠)(𝑡−𝑧)^1−𝛼𝑑𝑊−𝑑𝑊_𝑠²𝑑𝑧𝑑𝑡. (4.13) Taking expectations on both side, letting 𝑒=𝐸(∫¹₀𝜖²(𝑡)𝑑𝑡) and using the Burkholder-Gundy-type inequality (𝐸𝑋)²≤𝐸(𝑋²), we get 𝜙𝑒≤3𝐸1−3𝜆². (4.14) Applying Lemma 4.1, there is constant 𝐶 such that 𝐸¹₀(𝑢(𝑡)−𝑢(𝑡))²𝑑𝑡≤𝐶𝜌²1−3𝜆² (4.15) as we claimed. In the rest of this section, we construct our numerical method which enables us to solve (4.4) numerically.
Let {𝜓_𝑗(𝑡),𝑗=0,1,2,…,𝐽} be a set of deterministic orthogonal functions with weight function 𝜈(𝑡) and [𝑎,𝑏] its interval of orthogonality. Also, assume that ^𝑏_𝑎𝜈(𝑡)𝜓_𝑘(𝑡)𝜓_𝑗(𝑡)ℎ𝑑𝑡=^𝑘,𝑘=𝑗,0,𝑘≠𝑗. (4.16) Now, for 𝑡𝜖[0,1], we assume that 𝑢(𝑡)≈𝑢_𝐽(𝑡)=𝑢₀+𝜈(𝑡)^𝐽_𝑗=0𝑐_𝑗𝜓_𝑗(𝑡), (4.17) then (4.4) reduces to 𝜈(𝑡)^𝐽_𝑗=0𝑐_𝑗𝜓_𝑗(𝑡)=^𝑡₀1Γ(𝛼)(𝑡−𝑧)^1−𝛼𝑓(𝑧)+𝐴(𝑧)𝑢₀+^𝑧₀𝐵(𝑧,𝑠)𝑢₀𝑑𝑊_𝑠+1Γ(𝛼)^𝐽_𝑗=0𝑐_𝑗^𝑡₀1(𝑡−𝑧)^1−𝛼𝐴(𝑧)𝑢₀+𝜈(𝑧)𝜓_𝑗(𝑧)+^𝑧₀𝐵(𝑧,𝑠)𝑢₀+𝜈(𝑠)𝜓_𝑗(𝑠)𝑑𝑊_𝑠𝑑𝑧. (4.18) Multiply this equation by 𝜓_𝑘(𝑡) and integrating the resultant over the interval [𝑎,𝑏], we obtain that ℎ_𝑘−𝛾_𝑘𝑘𝑐_𝑘−^𝐽_{𝑗=0,𝑗≠𝑘}𝑐_𝑗𝛾_𝑘𝑗=𝛽_𝑘,𝑘=0,1,2,…,𝐽, (4.19) where 𝛽_𝑘=^𝑏_𝑎𝜓_𝑘(𝑡)𝐽^𝛼𝑉₀(𝑡),𝛾_𝑘𝑗=^𝑏_𝑎𝜓_𝑘(𝑡)𝐽^𝛼𝑉_𝑗(𝑡)𝑉𝑑𝑡,₀(𝑧)=𝑓(𝑧)+𝐴(𝑧)𝑢₀+^𝑧₀𝐵(𝑧,𝑠)𝑢₀𝑑𝑊_𝑠,𝑉_𝑗(𝑧)=𝐴(𝑧)𝜈(𝑧)𝜓_𝑗(𝑧)+^𝑧₀𝐵(𝑧,𝑠)𝜈(𝑠)𝜓_𝑗(𝑠)𝑑𝑊_𝑠 (4.20) for 𝑗,𝑘=0,1,2,…,𝐽. In case of white noise, we may evaluate the function 𝑉_𝑗(𝑧) if we regularize the stochastic term by replacing the white noise with a smoother stochastic term. In other words, the last two equations may be replaced with 𝑉₀(𝑧)≈𝑓(𝑧)+𝐴(𝑧)𝑢₀+^𝑁1_𝑛=1𝜂_𝑛^𝑧₀𝐵(𝑧,𝑠)𝑢₀𝜒_𝑛(𝑠)𝑉𝑑𝑠,_𝑗(𝑧)≈𝐴(𝑧)𝜈(𝑧)𝜓_𝑗(𝑧)+^𝑁2_𝑛=1𝜂^∗_𝑛^𝑧₀𝐵(𝑧,𝑡)𝜈(𝑡)𝜓_𝑗(𝑡)𝜒_𝑛(𝑠)𝑑𝑠. (4.21)

4.1. The Solution of the Stochastic Linear System

The linear system (4.19) may be rewritten as 𝑐_𝑘−^𝐽_𝑗=0𝛾^∗_𝑘𝑗𝑐_𝑗=𝛽^∗_𝑘,𝑘=0,1,2,…,𝐽, (4.22) where 𝛽^∗_𝑘=𝛽_𝑘ℎ_𝑘−𝛾_𝑘𝑘,𝛾^∗_𝑘𝑗=𝛾_𝑘𝑗ℎ_𝑘−𝛾_𝑘𝑘. (4.23)

Theorem 4.3. This system has a unique solution if 1−𝛾^∗_𝑖𝑖>^𝐽_𝑗≠𝑖𝛾^∗_𝑖𝑗,𝑖=0,1,2,…,𝐽. (4.24)

Proof. The above system can be written, in the matrix form as (𝐼−𝑃)𝐶=𝑄, (4.25) where 𝛾𝑃=^∗00𝛾^∗₀₁𝛾^∗₀₂⋯𝛾^∗_0𝐽𝛾^∗₁₀𝛾^∗₁₁𝛾^∗₁₂⋯𝛾^∗_1𝐽𝛾^∗₂₀𝛾^∗₂₁𝛾^∗₂₂⋯𝛾^∗_2𝐽𝛾⋯⋯⋯⋯⋯_∗𝐽0𝛾^∗_𝐽1𝛾^∗_𝐽2⋯𝛾^∗_𝐽𝐽,𝑄=𝑢₀𝛽^∗₀𝛽^∗₁𝛽^∗₂⋯𝛽^∗_𝐽, (4.26) and the unknowns matrix 𝐶=(𝑐₀𝑐₁𝑐₂⋯𝑐_𝐽)^𝑇. Cramer's rule gives the solution of this linear system; namely, 𝑐_𝑗=1det(𝐼−𝑃)^𝐽_𝑖=0det𝑃^∗_𝑖𝑗𝑄_𝑗,𝑗=0,1,2,…,𝐽, (4.27) where det(𝐼−𝑃) is the determinant of 𝐼−𝑃, 𝑄_𝑗 is the 𝑗th component of the matrix 𝑄, and det(𝑃^∗_𝑖𝑗) agrees with det(𝐼−𝑃) except for the 𝑖th column where all but the 𝑗th term is 0, and the 𝑗th term is 1, providing that det(𝐼−𝑃)≠0; see [12]. It is known that the inverse of the matrix 𝐼−𝜆𝑃 will always exist for all values of 𝜆 with the exception of at most 𝐽+1 values. These are the roots of the characteristic determinantal equation det(𝐼−𝜆𝑃)=0. Therefore, the condition det(𝐼−𝑃)≠0 holds if 𝜆=1 is not an eigenvalue of the matrix 𝑃. According to Gershgorin theorem, all eigenvalues of the matrix 𝑃 lie in the circles 𝑆_𝑖={𝑧∶|𝑧−𝛾^∗_𝑖𝑖|≤∑^𝐽_𝑗≠𝑖|𝛾^∗_𝑖𝑗|},𝑖=0,1,2,…𝐽. Hence the system (4.22) has a unique solution given by (4.27) if (4.24) is satisfied. As a conclusion of this theorem, our numerical method depends on the choice of the orthonormal set {𝜓_𝑗(𝑡),𝑗=0,1,2,…,𝐽}. This set should be chosen in such a way that the requirement of the above theorem, (4.24), should be satisfied.

5. Numerical Examples

In this section, we give two examples. The first example is fractional differential equation in which there is no stochastic term. The second example is a stochastic differential equation in which the differentiation is ordinary not of fractional order.

For each of these examples, we use the Jacobi polynomials of degree 𝑘, see [13] 𝐺_𝑘=Γ𝑝,𝑞,𝑥𝑞+𝑘Γ𝑝+2𝑘^𝑘_𝑚=0(−1)^𝑘𝑘𝑚Γ𝑝+2𝑘−𝑚Γ𝑥𝑞+𝑘−𝑚^𝑘−𝑚,=𝑘!Γ𝑝+𝑘Γ𝑃𝑝+2𝑘^{𝑘(𝑝−𝑞,𝑞−1)}(2𝑥−1),𝑝−𝑞>−1,𝑞>0 (5.1) with weight 𝜈(𝑡)=𝑥^𝑞−1(1−𝑥)^𝑝−𝑞 and satisfy the orthonormality condition ¹₀𝜈(𝑥)𝐺_𝑘𝑝,𝑞,𝑥𝐺_𝑗ℎ𝑝,𝑞,𝑥𝑑𝑥=^𝑘,𝑘=𝑗,0,𝑘≠𝑗, (5.2) where ℎ_𝑘=𝑘!Γ𝑝+𝑘Γ𝑞+𝑘Γ𝑝−𝑞+𝑘+12𝑘+𝑝Γ²𝑝+2𝑘. (5.3)

Example 5.1. Consider the fractional differential equation ^𝐶𝐷^𝛼𝑢(𝑡)(𝑡)Γ(2)+𝑢=𝑡+Γ()𝑡2−𝛼^1−𝛼,𝑢(0)=0,0≤𝑡≤1. (5.4) With exact solution 𝑢(𝑡)=𝑡, Table 1 gives the error (error = absolute value of the difference between the exact solution and the approximate solution) for different values of 𝐽 and Figure 1 compares the graph of the exact solution with the graphs of the approximate solutions for different values of 𝐽= 2, 4, and 8. As an application of this example, the logistic population growth 𝐷^𝛽𝑁(𝑡)=𝑟𝛼𝑁(𝑡𝑁(𝑡))1−𝐾^𝛼,0≤𝛽≤1. (5.5)

Table 1 

The error term, Example 5.1, for different values of 𝐽.

Figure 1 

Graph of the error term, Example 5.1, for different values of 𝐽.

Example 5.2. Lemma 3.1 shows that the FSIDE (3.2) can be rewritten as (3.6). Therefore, we consider the stochastic differential equation 𝑑𝑢(𝑡)=𝑢(𝑡)𝑑𝑊(𝑡), (5.6) whose exact solution is 𝑢(𝑡)=𝑢₀𝑒^{(−𝑡/2)+𝑊(𝑡)}. (5.7) Denoting by 𝑢_𝐽 the approximating solution as in (4.17), by 𝑢(𝑡) the exact solution, by 𝑒^𝐽_𝑁=max{|𝑢(𝑡)−𝑢_𝐽|,0≤𝑡≤0.5} the errors, and by 𝛼^𝐽=log₂𝑒^𝐽_𝑁/𝑒^𝐽_2𝑁 an estimate of a convergence order, the results are contained in Table 2.

Table 2 

The convergence order, Example 5.2.

6. Conclusion

Our presented numerical method is applied for many different FSIDE of the form (3.2). From our numerical computations, we see the following. (1)There is no restriction on the choice of the orthogonal polynomials. Moreover, the values of the parameters, 𝑝 and 𝑞, of the chosen orthogonal polynomials, see (5.1), do not play any role in the computations yet more information about the expected value of the exact solution will be helpful in determining appropriate values of these parameters. In all cases the assumed approximation should agree with the initial condition with and the exact solution. (2)If the FSIDE (3.2) is free of the stochastic term while the differentiation is of fractional order as in the case of the first example, the method works quite well. (3)In case of the existence of the stochastic term, as in the case of the second example, although the method is of higher-order accuracy yet in practice the obtained results are not quite well as with the previous case. (4)We note that for every path, even with a series of Monte Carlo simulations, the stochastic linear system, (4.22), yields a unique deterministic solution.

Acknowledgment

The authors would like to thank Professor Bill Mclean for his valuable discussions regarding this work.