# Finite Element Method for Linear Multiterm Fractional Differential Equations

Authors: Abdallah A. Badr
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

We consider the linear multiterm fractional differential equation (fDE). Existence and uniqueness of the solution of such equation are discussed. We apply the finite element method (FEM) to obtain the numerical solution of this equation using Galerkin approach. A comparison, through examples, between our techniques and other previous numerical methods is established.


---

## Body

## 1. Introduction

Recently, many applications in numerous fields of science, engineering, viscoelastic materials, signal processing, controlling, quantum mechanics, meteorology, finance, life science, applied mathematics, and economics have been remodeled in terms of fractional calculus where derivatives and integrals of fractional order are introduced and so differential equation of fractional order are involved in these models, see [1–4]. Fractional-order derivatives provide an excellent instrument for the description of memory and hereditary properties of various materials and processes. They have been successfully used to model many problems. As an example which will give us a physical understanding of the fractional derivatives: in dynamical systems with fractional-order derivatives, fractional-order derivatives have been successfully used to model damping forces with memory effect or to describe state feedback controllers. In particular, the BagleyTorvik equation with 1/2-order derivative or 3/2-order derivative describes motion of real physical systems, an immersed plate in a Newtonian fluid, and a gas in a fluid, respectively [5]. Recently, it is found in [6] that in fractional-order vibration systems of single degree of freedom, the term of fractional-order derivative whose order is between 0 and 2 acts always as damping force. In addition, almost all systems containing internal damping are not suitable to be described properly by the classical methods, but the fractional calculus represents one of the promising tools which describe such systems. Therefore, mainly, a considerable importance is given to the field of fractional calculus. For analytical solution of fDEs, we refer to a domain decomposition method [7], the homotopy-perturbation method [8], variational iteration method [9–11], the fractional complex transform [12, 13], and the exp-function method [14]. For the numerical solution of fDEs, many approaches has been considered, for example, FDM [15], wavelet operational method [16], and recently series solution [17]. Also many authors used the fact that the solution of a fDE is the same as the solution of a singular integral equation, and so they solved this integral equation instead, see [18].A multiterm fDE may take the form(1.1)Dαu(x)=F(x,u(x),Dα1u(x),…,Dαℓu(x)),a≤x≤b,n-1<α1<α2<⋯<αm<α<n,i=1,2,…,n=[α]+1,
with initial conditions
(1.2)Dα-ku(x)|x=a=bk,bk∈k=1,…,n,
where bk are given constants and (x) means the integer part of x.The operatorDα denotes the α-derivative of the function f(t). There are various ways of defining the derivative of a given function f(x) of order α. We mention only the following definition due to Caputo’s definition.Definition 1.1.
(1.3)
D
c
a
+
α
f
(
x
)
=
1
Γ
(
n
-
α
)
∫
a
x
f
(
n
)
(
τ
)
(
x
-
τ
)
α
-
n
+
1
d
τ
,
n
-
1
<
α
<
n
,
x
>
a
.
Similarly for Dcb-αf(x).
The advantages of Caputo’s approach is that the initial conditions for the fDE with Caputo’s definition take the same form as the initial conditions of differential equation of integer order.Examples 1.
From the previous definition, we deduce(i)
ifx>a,Dcα1=0,(ii)
ifx>b>a,Dcα(x-b)=(1/Γ(2-α))(x-b)1-α.In this paper, we write Dαf(x) instead of DCa+αf(x).


## 2. Results for the Linear Multiterm fDE

In this paper, we write a linear multiterm fDE with Caputo’s derivatives in the following form; because of its importance in fluid mechanics:(2.1)D1+αu(x)+∑m=1MAm(x)Dαmu(x)=f(x),1>x>0,n≥α,αm>n-1 and incorporated given initial conditions data
(2.2)uk(0)=bk,k=0,1,2,…,n-1.
Equations of the form of (1.1) and (2.1) have been studied extensively by many authors, see [18]. For our concerns, we state the following two theorems, see [19].Theorem 2.1 (Diethelm 2001).
Letu be the solution of (2.1) with initial conditions (2.2) and let v be the solution of
(2.3)D1+α-v(x)+∑m=1MAm(x)Dα-mv(x)=f(x),1>x>0,n≥α-,α-m>n-1, and incorporated given initial conditions data
(2.4)vk(0)=bk,k=0,1,
where |α--α|<ϵ,|α-m-αm|<ϵ. For T<∞, we have
(2.5)∥u-v∥L∞[0,T]=O(ϵ),ϵ→0.As a consequence of this theorem, we can assume that the fractional ordersα,αm, are irrational numbers.Theorem 2.2 (Diethelm 2001).
Let the functionf in (1.1) satisfy Lipschitz condition with Lipschitz constant L in all its arguments except for the first. Assume that the orders α,αm∈Q. Then (1.1) subject to (1.2) has a unique solution on the interval [0,T] of the real line.


## 3. Modified Galerkin Method

In this section, we present our approach by using FEM to get the numerical solution of the general linear multiterm fDE (2.1) with initial conditions (2.2) and we restrict our self to the case 0<α,αm≤1. To perform such approach, we segment the domain [0,1] into N linear elements, say ei,i=1,2,…,N,ei=(xi,xi+1) with x1=0 and xN+1=1. These points are called the nodal points. Let the length of each element ei be equal to ℓ=xi+1-xi=1/N. At each nodal point xi,i=1,2,…,N+1, we define the roof function Ri(x) as follows:
(3.1)R1={x-x2-ℓ,x1=0≤x<x2,0,else,RN+1={x-xNℓ,xN<x≤1-xN+1,0,else,
and for i≠1,N+1 we have
(3.2)Ri={x-xi+1-ℓ,x∈ei,x-xi-1ℓ,x∈ei-1.
Note that Ri(xi)=1 and Ri(xj≠i)=0.We assume the approximate solution of (2.1) is a linear combination of these roof functions Rj(x). In other words, let
(3.3)u(x)≈u~(x)=∑j=1N+1cjRj(x),
where cj are constants to be determined. We choose the constant cj such that
(3.4)D1+αu~(x)+∑m=1MAm(x)Dαmu~(x)-f(x)
is minimized. This quantity is minimized if
(3.5)∫01Ri(x){D1+αu~(x)+∑m=1MAm(x)Dαmu~(x)-f(x)}dx=0,i=1,2,…,N+1.
Integration by parts the first term of the above equation, we obtain
(3.6)-∫01dRi(x)dxDαu~(x)dx+∑m=1M∫01Am(x)Ri(x)Dαmu~(x)dx=∫01f(x)Ri(x)dx-[Ri(x)Dαu~(x)]01,i=1,2,…,N+1.
Using (3.3) in the last equation, we obtain
(3.7)∑j=1N+1cj{∫01dRi(x)dxDαRj(x)dx-∑m=1M∫01Am(x)Ri(x)DαmRj(x)dx}=-∫01f(x)Ri(x)dx-cN+1ℓΓ(2-α)δi,N+1,i=1,2,…,N+1,
where
(3.8)δi,N+1={1,i=N+1,0,else.
Note that the first and the last equations of the above equations are invalid.Now, define(3.9)Kij1=∫01dRi(x)dxDαRj(x)dx,Kij2=∑m=1M∫01Am(x)Ri(x)DαmRj(x)dx,Fi=-∫01f(x)Ri(x)dx,Bi=cN+1ℓΓ(2-α)δi,N+1.
Then (3.7) has the matrix form
(3.10)∑j=1N+1Kijcj=Fi,i=2,…,N,
where
(3.11)Kij=Kij1+Kij2+Bi.
Or, simply,
(3.12)[K]C=F.The matrixK=(Kij+Bi) is called the global stiffness matrix and the vector F=(Fi) is called the global force vector. In calculating the elements of the two matrices K and F, we have to integrate over each element eν,ν=1,2,…,N. Therefore, we express K as a sum of element stiffness matrixes K[eν] and F as a sum of element force vectors F[eν]. Namely,
(3.13)K=∑ν=1N+1K[eν],F=∑ν=1N+1F[eν],
where
(3.14)Kij[eν]=Kij1,[eν]+Kij2,[eν],Kij1,[eν]=∫eνdRi(x)dxDαRj(x)dx,Kij2,[eν]=∑m=1M∫eνAm(x)Ri(x)DαmRj(x)dx,Fi[eν]=-∫eνf(x)Ri(x)dx.
Using the properties of the roof functions, we obtain
(3.15)Kij1={0,|i-j|≥2,∫eidRi(x)dxDαRj(x)dx,j=i+1,∫ejdRi(x)dxDαRj(x)dx,i=j+1,(∫ei-1+∫ei‍)dRi(x)dxDαRi(x)dx,j=i.
Similarly,
(3.16)Kij2={0,|i-j|≥2,∑m=1M∫eiAm(x)Ri(x)DαmRj(x)dx,j=i+1,∑m=1M∫ejAm(x)Ri(x)DαmRj(x)dx,i=j+1,(∫ei-1+∫ei‍)Am(x)Ri(x)DαmRj(x)dx,j=i,Fi=-(∫ei-1+∫ei‍)f(x)Ri(x)dx.
Note that ∫e0G(x)dx=0=∫en+1G(x)dx, for any function G(x).Also, the initial conditions given by (1.2) gives
(3.17)c1=b0,c2=ℓb1+b0.
Therefore, the final linear system which gives the unknown ci,i=2,3,4,…,N takes the form
(3.18)∑j=3N+1Kijcj=Fi*,i=2,3,4,…,N,
where
(3.19)Kij=Kij1+Kij2,Fi*=Fi-Ki1c1-Ki2ci2.


## 4. Numerical Experiments

We consider the following Cauchy problem:(4.1)D1.5y(t)+2D1y(t)+3tD0.5y(t)+(1-t)y(t)=f(t),y(0)=y′(0)=0,
where
(4.2)f(t)=2Γ(1.5)t+4t+4Γ(1.5)t2+(1-t)t2.
It is easy to check that all assumptions of Theorem 2.2 are fulfilled. The exact solution of the fractional differential equation is y(t)=t2.This problem was solved numerically by the modified Galerkin method on the interval[0,1] using different values of N, the number of nodal points. In Table 1, some results for different values of the parameters N are presented.Table 1
N
e
N (error)
α
N (convergence order)
4
8.5231043E−2
2.979325336
8
4.3320479E−3
3.409608789
16
1.4319207E−4
3.416117538
32
4.7023827E−6
3.850441833
64
1.00021262E−7
3.896749819
128
2.03121262E−9Denoting byeN=max{|u(x)-uN(x)|,0<x<1} the errors and by αN=log(eN/e2N) an estimate of a convergence order, the results is contained in Table 1.


## 5. Discussion

In general, finding the exact solutions of fractional differential equations is difficult and needs more computational work or mostly impossible. In this study, the finite element method is generalized and applied to fDE with multilinear terms. The method described in this paper considers only fDE of the form of (2.1) with conditions of the form of (2.2), but the basic ground work has been laid for extension to any fDE (linear or nonlinear) with any initial (or boundary) conditions. The roof functions defined by (3.1), (3.2) are chosen to be linear, yet we could choose them to be of higher order; quadratic, cubic,…, and so forth. Singularities of the fDE is the key behind the difficulties of the numerical solution of such equation. In our approach for solving fDE, such difficulties has been eliminated. The obtained linear system is easy to solve since the coefficients matrix, K, is a tridiagonal matrix. Moreover, the coefficients of the matrices K and F, given by (3.15), (3.16), are easily computable explicitly either by hand or by using software such as Maple which can calculate them symbolically. In solving differential equations of integer order using the modified Galerkin techniques described above in conjunction with piecewise linear shape functions, the terms derivatives of order m (with m greater than 2) in the given differential equation would make no contribution to the approximation leading to a poor result. In contrast to this situation, derivatives of fractional order in the fDE will have contributions even with linear shape functions.Also, the described method gave us a good agreement with other numerical methods with a relatively simple procedure and little computational efforts. It is also noted that this procedure transforms linear differential equations into an algebraic system, which depends on the roof functions. Therefore, it can provide with some advantages in writing computer codes of the desired system.



---

*Source: 482890-2012-10-24.xml*