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

Authors: Abdallah Ali Badr; Hanan Salem El-Hoety
Journal: Mathematical Problems in Engineering
 (2012)

---

## 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.


---

## Body

## 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
(1.1)uα(t)=f(t,u)+∫t0th(t,s,u)ds,u(t0)=u0,t0≤t≤T,0≤α≤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
(1.2)uα(t;w)=f(t,u;w)+∫t0th(t,s,u;w)dWs,u(t0)=u0,t0≤t≤T,0≤α≤1,
where W is a Brownian motion.In this paper we study the FSIDE:(1.3)uα(t;w)+A(t)u(t;w)=∫t0tB(t,s;w)u(s;w)dWs+f(t;w),t>t0
with initial condition u(t0;w)=u0 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)du(t)=a(t,u;w)dt+b(t,u;w)dWs.
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Ω={ω1,ω2,…} be the collection of all outcomes, ωi,i=1,2,… of an experiment, and let A⊂Ω. Define the measure P:A→[0,1] with P(Ω)=1. The triplet (Ω,A,P) is called a probability space.Definition 2.1.
A collection{X(t,·):t≥0} of random variables X:Ω×T→ℛ, T=[0,∞] is called a stochastic process. For each point ω∈Ω, the mapping t→X(t,ω) is a realization, sample path, or trajectory of the stochastic process.Definition 2.2.
A real-valued stochastic process that depends continuously ont∈[0,T]W={W(t),t≥0},W:T×Ω→ℛ is called a Brownian motion (BM) or Wiener process if  (1)
W(0)=0 a.s.,(2)
W(t)-W(s) is N(0,t-s), the Gaussian distribution with mean 0 and variance t-s, for all t>s>0,(3)
for all times0<t1<t2<⋯<tn, the random variables W(t1),W(t2)-W(t1),…,W(tn)-W(tn-1) are independents.Also, the expected values of W(t) and W2(t) are given by E(W(t))=0 and E(W2(t))=t for each time t>0.Definition 2.3.
(i) Suppose that{At≥0} is an increasing family of σ-algebras of A such that Wt is At-measurable E(W(t)∣A0)=0 and E(W(t)-W(s)∣As)=0 with probability one, w.p.1, for all 0≤s≤t≤T.
(ii) Letℒ be the σ-algebra of Lebesgue subsets of ℛ. Let ℒT2 be a class of functions f:[0,T]×Ω→ℛ satisfying the following:  (a)
f is jointly ℒ×A-measurable;(b)
∫0TE((f(t,·))2)dt<∞;(c)
E((f(t,·))2) for each 0≤s≤t≤T;(d)
f(t,·) is At measurable for each 0≤t≤T.Definition 2.4 (Convergence).
Ifu(h) is the exact value of a random variable and ũn(h) is its approximate value, then we have  (1)
strong convergence when(2.1)E(|u(h)-ũn(h)|)⟶0,n⟶∞,(2)
weak convergence when(2.2)|E(u(h))-E(ũn(h))|⟶0,n⟶∞,(3)
mean square convergence when(2.3)E(u(h)-ũn(h))2⟶0,n⟶∞.Definition 2.5.
Lety(x)∈L1([a,b]), we define the fractional integral of order α, 1≥α>0, of the function y(x) over the interval [a,b] as
(2.4)Ja+αy(x)=1Γ(α)∫axy(t)(x-t)1-αdt,x>a,Jb-αy(x)=1Γ(α)∫xby(t)(t-x)1-αdt,x<b.
For completion, we define J0=I (identity operator); that is, we mean J0y(x)=y(x). Furthermore, by Ja+αy(x) we mean the limit (if it exists) of Jcα as c→a+. This definition is according to Riemann-Liouville definition of fractional integral of arbitrary order α. For simplicity, we define the fractional integral by the equation
(2.5)Jαy(x)=1Γ(α)∫axy(t)(x-t)1-αdt,x>a,0<α≤1,
where we dropped a+. The folllowing equation
(2.6)Dαy(x)=ddxJ1-αy(x)
denotes the fractional derivative of y(x) of order α. For numerical computations, we usually use another definition of the fractional derivative, due to Caputo [2], given by
(2.7)(CDαy)(x)=(Dα(y(t)-y(a)))(x).


## 3. Formulation of the Problem

In this paper, we present a numerical method to solve the following FSIDE:(3.1)CDtαu(t;w)-A(t)u(t;w)=∫t0tB(t,s;w)u(s;w)dWs+f(t;w),t>t0
with initial condition u(t0;w)=ut0.For simplicity of notations, we drop the variablew, and so this equation takes the form (3.2)CDαu(t)+A(t)u(t)=∫t0tB(t,s)u(s)dWs+f(t),t>t0,u(t0)=ut0,
where CDα is the Caputo-fractional derivative operator of order α with 0<α≤1,t∈[t0,T]. The integral term of the right hand side is an Itô integral, W is a BM, and dW/dt is a white noise, the derivative of a BM (the derivative in the sense of distribution). The functions A(t),B(t,s;w),f(t;w),ut0(w), and u(t;w) satisfy the following conditions:  C1 (measurability):A=A(t),B=B(s,t) and u(t) are ℒ2-measurable in [0,1]×ℛ;C2: there exists constantsK1,K2,K3>0 such that
(3.3)|A(t)(t-s)1-α|<K1,|f(t)(t-s)1-α|<K2,|∫stB(z,s)dz(t-z)1-α|<K3
for all t0≤s≤t≤T and 0<α≤1;C3 (initial value):ut0 is At0-measurable with E(ut0)<∞.Letk=max{k1,k2,k3}/Γ(α).Equation (3.2) may be written as
(3.4)u(t)=ut0+CD-αf(t)+CD-αA(t)u(t)+CD-α∫t0tB(t,s)u(s)dW(s).
Using the definition of the operator CD-α, (2.7), we obtain that (3.5)u(t)=ut0+1Γ(α)∫t0tf(z)dz(t-z)1-α+1Γ(α)∫t0tA(z)u(z)dz(t-z)1-α+1Γ(α)∫t0t∫0zB(z,s)u(s)(t-z)1-αdWsdz.
This is a Fredholm integral equation that we are going to solve instead of solving (3.2).Without loss of generality, we will lett0=0 and T=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
(3.6)u(t)=ut0+∫t0ta(t,s,us)ds+∫t0tb(t,s,us)dWs,
where the drift is given by
(3.7)a(t,s,us)={1Γ(α){A(s)us+f(s)(t-s)1-α}0≤s<t≤1,0else,
and the diffusion is
(3.8)b(t,s,us)={usΓ(α)∫z=st{B(z,s)dz(t-z)1-α}0≤s<t≤1,0else.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 whereA1 (measurability): the drift and the diffusion are jointlyℒ2-measurable in [0,1]×ℛ;A2 (lipschitz condition): there is a constantK>0 such that
(3.9)|a(t,s,u)-a(t,s,v)|≤|u-v|,|b(t,s,u)-b(t,s,v)|≤|u-v|,
for 0≤s≤t≤1,u,v∈R;A3 (linear growth bound): there is a constantK>0 such that
(3.10)|a(t,s,u)|2≤K2(1+|u|2),|b(t,s,u)|2≤K2|u|2;A4 (Initial value):u0 is A0-measurable with E(u0)<∞.
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 ut 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 processẆ(t),t∈[0,1].Consider the uniform time discretization of the interval[0,1](4.1)Iρ={In:In=[tn-1,tn),n=[1:N],tn-tn-1=ρ,t0=0,tN=1}.
Then the following approximation is defined for the white noise process Ẇ(t) on this discretization, (4.2)dWtdt≈dŴtdt=∑n=1Nηnχn(t),
where the coefficients are random variables defined by
(4.3)ηn=1ρ∫Iρχn(t)dW(t)∈N(0,1ρ),χn(t)={1,tn-1≤t≤tn,0,else.
As a direct result of this approximation of the white noise, it is easily to see that E(∫ℑρg[dWt-dŴt])=0=E((∫ℑρg[dWt-dŴt])2) for any bounded function g.NowdŴt can be substituted for dWt to obtain the following smoothed version of the SDE (3.5):
(4.4)û(t)=u0+1Γ(α)∫0tf(z)dz(t-z)1-α+1Γ(α)∫0tA(z)û(z)dz(t-z)1-α+1Γ(α)∫0t∫0zB(z,s)û(s)(t-z)1-αdŴsdz.
The solution of this equation, û(t), is smoother than u(t) and therefore standard numerical procedures can be applied to compute its approximate solution. Once the approximate solution, û, is obtained for M realizations of such approximate solution, the Monte Carlo method then uses these approximations to compute corresponding sample averages of these M realizations.Now we show thatû(t) is a good approximation of u(t). To show this the following lemma is required.Lemma 4.1.
Given a nonrandom functionH(t) such that
(4.5)|H(t)-H(τ)|≤λ|t-τ|δ,0≤t,τ≤1,0≤δ≤1,
where λ>0 is a constant, then
(4.6)E[∫0zH(t)dWt-∫0zH(τ)dWτ]2≤λ2(ρ)2δ
with ρ≥|t-τ| for 0≤t,τ,z≤1.Proof.
Let(4.7)Iρ*={In*:In*=[tn-1*,tn*),n=[1:N*],tn*-tn-1*=ρ*,t0=0,tN*=z},ρ*≤ρ,
be a partition of the interval [0,z]. We have
(4.8)E[∫0zH(t)dWt-∫0zH(τ)dWτ]2=E[∑n=1N*∫tn-1*tn*(H(t)-1ρ∫tn-1*tn*H(τ)dτ)dWt]2=∑n=1N*∫tn-1*tn*[H(t)-1ρ∫tn-1*tn*H(τ)dτ]2dt=∑n=1N*∫tn-1*tn*[1ρ∫tn-1*tn*{H(t)-H(τ)}dτ]2dt≤λ2(ρ)2∑n=1N*∫tn-1*tn*[1ρ∫tn-1*tn*{|t-τ|δ}dτ]2dt≤λ2ρ*2δ≤λ2ρ2δ.
With the help of this lemma, we can show that u(t)→û(t) as ρ→0.Theorem 4.2.
(4.9)E[∫01(u(t)-û(t))2dt]≤Cρ21-3λ2,
where
(4.10)λ2=1Γ2(α)∫01∫0t1(t-z)2-2α{A2(z)+∫0zB2(z,s)ds}dzdt.Proof.
Letϵ(t)=u(t)-û(t), then (3.5) and (4.4) lead to
(4.11)ϵ(t)=1Γ(α)∫0tA(z)ϵ(z)dz(t-z)1-α+1Γ(α)∫0t∫0zB(z,s)ϵ(s)(t-z)1-αdWsdz+1Γ(α)∫0t∫0zB(z,s)û(s)(t-z)1-α[dW-dŴs]dz.
Applying the inequality (a+b+c)2≤3(a2+b2+c2) and the Hölder's inequality to this equation, we obtain
(4.12)∫01ϵ2(t)dt≤3(λ2∫01ϵ2(t)dt+ϕ),
where λ2 is given by (4.10) and
(4.13)ϕ=1Γ2(α)∫01∫0t∫0z|B(z,s)û(s)(t-z)1-α[dW-dŴs]|2dzdt.
Taking expectations on both side, letting e=E(∫01ϵ2(t)dt) and using the Burkholder-Gundy-type inequality (EX)2≤E(X2), we get
(4.14)e≤3E(ϕ)1-3λ2.
Applying Lemma 4.1, there is constant C such that
(4.15)E[∫01(u(t)-û(t))2dt]≤Cρ21-3λ2
as we claimed. In the rest of this section, we construct our numerical method which enables us to solve (4.4) numerically.
Let{ψj(t),j=0,1,2,…,J} be a set of deterministic orthogonal functions with weight function ν(t) and [a,b] its interval of orthogonality. Also, assume that
(4.16)∫abν(t)ψk(t)ψj(t)dt={hk,k=j,0,k≠j.
Now, for tϵ[0,1], we assume that
(4.17)u(t)≈ûJ(t)=u0+ν(t)∑j=0Jcjψj(t),
then (4.4) reduces to
(4.18)ν(t)∑j=0Jcjψj(t)=∫0t1Γ(α)(t-z)1-α{f(z)+A(z)u0+∫0zB(z,s)u0dŴs}+1Γ(α)∑j=0Jcj∫0t1(t-z)1-α[A(z)u0+ν(z)ψj(z)+∫0zB(z,s)u0+ν(s)ψj(s)dŴs]dz.
Multiply this equation by ψk(t) and integrating the resultant over the interval [a,b], we obtain that
(4.19)(hk-γkk)ck-∑j=0,j≠kJcjγkj=βk,k=0,1,2,…,J,
where
(4.20)βk=∫abψk(t)JαV0(t),γkj=∫abψk(t)JαVj(t)dt,V0(z)=f(z)+A(z)u0+∫0zB(z,s)u0dŴs,Vj(z)=A(z)ν(z)ψj(z)+∫0zB(z,s)ν(s)ψj(s)dŴs
for j,k=0,1,2,…,J. In case of white noise, we may evaluate the function Vj(z) 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
(4.21)V0(z)≈f(z)+A(z)u0+∑n=1N1ηn∫0zB(z,s)u0χn(s)ds,Vj(z)≈A(z)ν(z)ψj(z)+∑n=1N2ηn*∫0zB(z,t)ν(t)ψj(t)χn(s)ds.


### 4.1. The Solution of the Stochastic Linear System

The linear system (4.19) may be rewritten as (4.22)ck-∑j=0Jγkj*cj=βk*,k=0,1,2,…,J,
where
(4.23)βk*=βk(hk-γkk),γkj*=γkj(hk-γkk).Theorem 4.3.
This system has a unique solution if(4.24)|1-γii*|>∑j≠iJ|γij*|,i=0,1,2,…,J.Proof.
The above system can be written, in the matrix form as(4.25)(I-P)C=Q,
where
(4.26)P=(γ00*γ01*γ02*⋯γ0J*γ10*γ11*γ12*⋯γ1J*γ20*γ21*γ22*⋯γ2J*⋯⋯⋯⋯⋯γJ0*γJ1*γJ2*⋯γJJ*),Q=u0(β0*β1*β2*⋯βJ*),
and the unknowns matrix C=(c0c1c2⋯cJ)T. Cramer's rule gives the solution of this linear system; namely,
(4.27)cj=1det⁡⁡(I-P)∑i=0Jdet⁡⁡(Pij*)Qj,j=0,1,2,…,J,
where det(I-P) is the determinant of I-P, Qj is the jth component of the matrix Q, and det(Pij*) agrees with det(I-P) except for the ith column where all but the jth term is 0, and the jth term is 1, providing that det(I-P)≠0; see [12]. It is known that the inverse of the matrix I-λP will always exist for all values of λ with the exception of at most J+1 values. These are the roots of the characteristic determinantal equation det(I-λP)=0. Therefore, the condition det(I-P)≠0 holds if λ=1 is not an eigenvalue of the matrix P. According to Gershgorin theorem, all eigenvalues of the matrix P lie in the circles Si={z:|z-γii*|≤∑j≠iJ|γij*|},i=0,1,2,…J. 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 {ψj(t),j=0,1,2,…,J}. This set should be chosen in such a way that the requirement of the above theorem, (4.24), should be satisfied.


## 4.1. The Solution of the Stochastic Linear System

The linear system (4.19) may be rewritten as (4.22)ck-∑j=0Jγkj*cj=βk*,k=0,1,2,…,J,
where
(4.23)βk*=βk(hk-γkk),γkj*=γkj(hk-γkk).Theorem 4.3.
This system has a unique solution if(4.24)|1-γii*|>∑j≠iJ|γij*|,i=0,1,2,…,J.Proof.
The above system can be written, in the matrix form as(4.25)(I-P)C=Q,
where
(4.26)P=(γ00*γ01*γ02*⋯γ0J*γ10*γ11*γ12*⋯γ1J*γ20*γ21*γ22*⋯γ2J*⋯⋯⋯⋯⋯γJ0*γJ1*γJ2*⋯γJJ*),Q=u0(β0*β1*β2*⋯βJ*),
and the unknowns matrix C=(c0c1c2⋯cJ)T. Cramer's rule gives the solution of this linear system; namely,
(4.27)cj=1det⁡⁡(I-P)∑i=0Jdet⁡⁡(Pij*)Qj,j=0,1,2,…,J,
where det(I-P) is the determinant of I-P, Qj is the jth component of the matrix Q, and det(Pij*) agrees with det(I-P) except for the ith column where all but the jth term is 0, and the jth term is 1, providing that det(I-P)≠0; see [12]. It is known that the inverse of the matrix I-λP will always exist for all values of λ with the exception of at most J+1 values. These are the roots of the characteristic determinantal equation det(I-λP)=0. Therefore, the condition det(I-P)≠0 holds if λ=1 is not an eigenvalue of the matrix P. According to Gershgorin theorem, all eigenvalues of the matrix P lie in the circles Si={z:|z-γii*|≤∑j≠iJ|γij*|},i=0,1,2,…J. 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 {ψj(t),j=0,1,2,…,J}. 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 degreek, see [13]
(5.1)Gk(p,q,x)=Γ(q+k)Γ(p+2k)∑m=0k(-1)k(km)Γ(p+2k-m)Γ(q+k-m)xk-m,=k!Γ(p+k)Γ(p+2k)Pk(p-q,q-1)(2x-1),p-q>-1,q>0
with weight ν(t)=xq-1(1-x)p-q and satisfy the orthonormality condition (5.2)∫01ν(x)Gk(p,q,x)Gj(p,q,x)dx={hk,k=j,0,k≠j,
where
(5.3)hk=k!Γ(p+k)Γ(q+k)Γ(p-q+k+1)(2k+p)Γ2(p+2k).Example 5.1.
Consider the fractional differential equation(5.4)CDαu(t)+u(t)=t+Γ(2)Γ(2-α)t1-α,u(0)=0,0≤t≤1.
With exact solution u(t)=t, Table 1 gives the error (error = absolute value of the difference between the exact solution and the approximate solution) for different values of J and Figure 1 compares the graph of the exact solution with the graphs of the approximate solutions for different values of J= 2, 4, and 8. As an application of this example, the logistic population growth
(5.5)DβN(t)=rαN(t)[1-(N(t)K)α],0≤β≤1.Table 1
The error term, Example5.1, for different values of J.
tJ=2J=4J=80.010.023860826300.007392765000.001119306570.080.02963852220.000362107730.000423089070.120.02398314610.00309400540.00045735510.180.01416987450.00421700170.00055115020.230.00628046080.00344455840.00001469720.350.00880944210.00046986180.00038131400.460.01660342670.00288838600.00038675470.520.01827180280.00304017520.00039283410.660.01510024590.00050021610.00039807820.730.00987789570.00141507460.00022833390.870.00759823760.00260859030.00029758510.960.02362977480.00176605640.00046625271.000.0319257630.0058847590.001152260Figure 1
Graph of the error term, Example5.1, for different values of J.Example 5.2.
Lemma3.1 shows that the FSIDE (3.2) can be rewritten as (3.6). Therefore, we consider the stochastic differential equation
(5.6)du(t)=u(t)dW(t),
whose exact solution is
(5.7)u(t)=u0e(-t/2)+W(t).
Denoting by ûJ the approximating solution as in (4.17), by u(t) the exact solution, by eNJ=max{|u(t)-ûJ|,0≤t≤0.5} the errors, and by αJ=log2eNJ/e2NJ an estimate of a convergence order, the results are contained in Table 2.Table 2
The convergence order, Example5.2.
NeN2α2eN4α4eN8α81280.62401.2285.135×10-21.56038.531×10-32.9132560.51031.29193.291×10-21.7322.928×10-33.3215120.3951.50191.910×10-21.8737.533×10-43.45110240.2631.014×10-22.183×10-4


## 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,p and q, 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.



---

*Source: 709106-2012-04-12.xml*