Finite Element Method for Linear Multiterm Fractional Differential Equations

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.


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 fractionalorder 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 with initial conditions where b k are given constants and x means the integer part of x.
The operator D α 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. 1.3 Similarly for c D α b− 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
In this paper, we write

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: 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 .Let u be the solution of 2.1 with initial conditions 2.2 and let v be the solution of and incorporated given initial conditions data As a consequence of this theorem, we can assume that the fractional orders α, α m , are irrational numbers.

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 e i , i 1, 2, . . ., N, e i x i , x i 1 with x 1 0 and x N 1 1.These points are called the nodal points.Let the length of each element e i be equal to x i 1 − x i 1/N.At each nodal point x i , i 1, 2, . . ., N 1, we define the roof function R i x as follows: and for i / 1, N 1 we have

3.2
Note that R i x i 1 and R i x j / i 0. We assume the approximate solution of 2.1 is a linear combination of these roof functions R j x .In other words, let where c j are constants to be determined.We choose the constant c j such that A m x D α m u x − f x dx 0, i 1, 2, . . ., N 1.

3.5
Integration by parts the first term of the above equation, we obtain

3.6
Using 3.3 in the last equation, we obtain where Note that the first and the last equations of the above equations are invalid.Now, define 3.9 Then 3.7 has the matrix form where Or, simply, K C F.

3.12
The matrix K K ij B i is called the global stiffness matrix and the vector F F i 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, where

3.14
Using the properties of the roof functions, we obtain

3.16
Note that e 0 G x dx 0 e n 1 G x dx, for any function G x .Also, the initial conditions given by 1.2 gives Therefore, the final linear system which gives the unknown c i , i 2, 3, 4, . . ., N takes the form where

Numerical Experiments
We consider the following Cauchy problem: where 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 t 2 .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.
Denoting by e N max{|u x − u N x |, 0 < x < 1} the errors and by α N log e N /e 2N an estimate of a convergence order, the results is contained in Table 1.

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.

Theorem 2 . 2
Diethelm 2001  .Let the function f 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.