Finite Element Method for Linear Multiterm Fractional Differential Equations

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.

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<𝛼₁<𝛼₂<⋯<𝛼_𝑚<𝛼<𝑛,𝑖=1,2,…,𝑛=[𝛼]+1, (1.1) with initial conditions 𝐷^𝛼−𝑘𝑢(𝑥)_𝑥=𝑎=𝑏_𝑘,𝑏_𝑘∈𝑘=1,…,𝑛, (1.2) where 𝑏_𝑘 are given constants and (𝑥) means the integer part of 𝑥.

The operator 𝐷^𝛼 denotes the 𝛼-derivative of the function 𝑓(𝑡). There are various ways of defining the derivative of a given function 𝑓(𝑥) of order 𝛼. We mention only the following definition due to Caputo’s definition.

Definition 1.1. ^𝑐𝐷^𝛼_𝑎+𝑓(𝑥)=1Γ(𝑛−𝛼)^𝑥_𝑎𝑓^(𝑛)(𝜏)(𝑥−𝜏)^{𝛼−𝑛+1}𝑑𝜏,𝑛−1<𝛼<𝑛,𝑥>𝑎. (1.3) Similarly for ^𝑐𝐷^𝛼_𝑏−𝑓(𝑥).
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)if 𝑥>𝑎,^𝑐𝐷^𝛼1=0,(ii)if 𝑥>𝑏>𝑎,^𝑐𝐷^𝛼(𝑥−𝑏)=(1/Γ(2−𝛼))(𝑥−𝑏)^1−𝛼.In this paper, we write 𝐷^𝛼𝑓(𝑥) instead of ^𝐶𝐷^𝛼_𝑎+𝑓(𝑥).

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: 𝐷^1+𝛼𝑢(𝑥)+^𝑀_𝑚=1𝐴_𝑚(𝑥)𝐷^𝛼𝑚𝑢(𝑥)=𝑓(𝑥),1>𝑥>0, (2.1)𝑛≥𝛼,𝛼_𝑚>𝑛−1 and incorporated given initial conditions data 𝑢^𝑘(0)=𝑏_𝑘,𝑘=0,1,2,…,𝑛−1. (2.2) 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 𝑢 be the solution of (2.1) with initial conditions (2.2) and let 𝑣 be the solution of 𝐷^1+𝛼𝑣(𝑥)+^𝑀_𝑚=1𝐴_𝑚(𝑥)𝐷^𝛼𝑚𝑣(𝑥)=𝑓(𝑥),1>𝑥>0, (2.3)𝑛≥𝛼,𝛼_𝑚>𝑛−1, and incorporated given initial conditions data 𝑣^𝑘(0)=𝑏_𝑘,𝑘=0,1, (2.4) where |𝛼−𝛼|<𝜖,|𝛼_𝑚−𝛼_𝑚|<𝜖. For 𝑇<∞, we have ‖𝑢−𝑣‖_{𝐿∞[0,𝑇]}=𝑂(𝜖),𝜖→0. (2.5)

As a consequence of this theorem, we can assume that the fractional orders 𝛼,𝛼_𝑚, are irrational numbers.

Theorem 2.2 (Diethelm 2001). Let the function 𝑓 in (1.1) satisfy Lipschitz condition with Lipschitz constant 𝐿 in all its arguments except for the first. Assume that the orders 𝛼,𝛼_𝑚∈𝑄. Then (1.1) subject to (1.2) has a unique solution on the interval [0,𝑇] 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<𝛼,𝛼_𝑚≤1. To perform such approach, we segment the domain [0,1] into 𝑁 linear elements, say 𝑒_𝑖,𝑖=1,2,…,𝑁,𝑒_𝑖=(𝑥_𝑖,𝑥_𝑖+1) with 𝑥₁=0 and 𝑥_𝑁+1=1. These points are called the nodal points. Let the length of each element 𝑒_𝑖 be equal to ℓ=𝑥_𝑖+1−𝑥_𝑖=1/𝑁. At each nodal point 𝑥_𝑖,𝑖=1,2,…,𝑁+1, we define the roof function 𝑅_𝑖(𝑥) as follows: 𝑅₁=𝑥−𝑥₂−ℓ,𝑥₁=0≤𝑥<𝑥₂,𝑅0,else,_𝑁+1=𝑥−𝑥_𝑁ℓ,𝑥_𝑁<𝑥≤1−𝑥_𝑁+1,0,else, (3.1) and for 𝑖≠1,𝑁+1 we have 𝑅_𝑖=𝑥−𝑥_𝑖+1−ℓ,𝑥∈𝑒_𝑖,𝑥−𝑥_𝑖−1ℓ,𝑥∈𝑒_𝑖−1. (3.2) Note that 𝑅_𝑖(𝑥_𝑖)=1 and 𝑅_𝑖(𝑥_𝑗≠𝑖)=0.

We assume the approximate solution of (2.1) is a linear combination of these roof functions 𝑅_𝑗(𝑥). In other words, let 𝑢(𝑥)≈𝑢(𝑥)=^𝑁+1_𝑗=1𝑐_𝑗𝑅_𝑗(𝑥), (3.3) where 𝑐_𝑗 are constants to be determined. We choose the constant 𝑐_𝑗 such that 𝐷^1+𝛼𝑢(𝑥)+^𝑀_𝑚=1𝐴_𝑚(𝑥)𝐷^𝛼𝑚𝑢(𝑥)−𝑓(𝑥) (3.4) is minimized. This quantity is minimized if ¹₀𝑅_𝑖(𝑥)𝐷^1+𝛼𝑢(𝑥)+^𝑀_𝑚=1𝐴_𝑚(𝑥)𝐷^𝛼𝑚𝑢(𝑥)−𝑓(𝑥)𝑑𝑥=0,𝑖=1,2,…,𝑁+1. (3.5) Integration by parts the first term of the above equation, we obtain −¹₀𝑑𝑅_𝑖(𝑥)𝐷𝑑𝑥^𝛼𝑢(𝑥)𝑑𝑥+^𝑀_𝑚=1¹₀𝐴_𝑚(𝑥)𝑅_𝑖(𝑥)𝐷^𝛼𝑚𝑢(𝑥)𝑑𝑥=¹₀𝑓(𝑥)𝑅_𝑖(𝑥)𝑑𝑥−𝑅_𝑖(𝑥)𝐷^𝛼𝑢(𝑥)¹₀,𝑖=1,2,…,𝑁+1. (3.6) Using (3.3) in the last equation, we obtain ^𝑁+1_𝑗=1𝑐_𝑗¹₀𝑑𝑅_𝑖(𝑥)𝐷𝑑𝑥^𝛼𝑅_𝑗(𝑥)𝑑𝑥−^𝑀_𝑚=1¹₀𝐴_𝑚(𝑥)𝑅_𝑖(𝑥)𝐷^𝛼𝑚𝑅_𝑗(𝑥)𝑑𝑥=−¹₀𝑓(𝑥)𝑅_𝑖(𝑥)𝑐𝑑𝑥−^𝑁+1()𝛿ℓΓ2−𝛼^𝑖,𝑁+1,𝑖=1,2,…,𝑁+1, (3.7) where 𝛿_𝑖,𝑁+1=1,𝑖=𝑁+1,0,else. (3.8) Note that the first and the last equations of the above equations are invalid.

Now, define 𝐾¹_𝑖𝑗=¹₀𝑑𝑅_𝑖(𝑥)𝐷𝑑𝑥^𝛼𝑅_𝑗(𝑥)𝐾𝑑𝑥,_2𝑖𝑗=^𝑀_𝑚=1¹₀𝐴_𝑚(𝑥)𝑅_𝑖(𝑥)𝐷^𝛼𝑚𝑅_𝑗(𝑥)𝐹𝑑𝑥,_𝑖=−¹₀𝑓(𝑥)𝑅_𝑖(𝑥)𝐵𝑑𝑥,_𝑖=𝑐_𝑁+1()𝛿ℓΓ2−𝛼^𝑖,𝑁+1. (3.9) Then (3.7) has the matrix form ^𝑁+1_𝑗=1𝐾_𝑖𝑗𝑐_𝑗=𝐹_𝑖,𝑖=2,…,𝑁, (3.10) where 𝐾_𝑖𝑗=𝐾¹_𝑖𝑗+𝐾²_𝑖𝑗+𝐵_𝑖. (3.11) Or, simply, [𝐾]𝐶=𝐹. (3.12)

The matrix 𝐾=(𝐾_𝑖𝑗+𝐵_𝑖) is called the global stiffness matrix and the vector 𝐹=(𝐹_𝑖) is called the global force vector. In calculating the elements of the two matrices 𝐾 and 𝐹, we have to integrate over each element 𝑒_𝜈,𝜈=1,2,…,𝑁. Therefore, we express 𝐾 as a sum of element stiffness matrixes 𝐾^[𝑒𝜈] and 𝐹 as a sum of element force vectors 𝐹^[𝑒𝜈]. Namely, 𝐾=^𝑁+1_𝜈=1𝐾^[𝑒𝜈],𝐹=^𝑁+1_𝜈=1𝐹^[𝑒𝜈], (3.13) where 𝐾^[𝑒𝜈]_𝑖𝑗=𝐾^1,[𝑒𝜈]_𝑖𝑗+𝐾^2,[𝑒𝜈]_𝑖𝑗,𝐾^1,[𝑒𝜈]_𝑖𝑗=_𝑒𝜈𝑑𝑅_𝑖(𝑥)𝐷𝑑𝑥^𝛼𝑅_𝑗(𝑥)𝐾𝑑𝑥,_{2,[𝑒𝜈]𝑖𝑗}=^𝑀_{𝑚=1𝑒𝜈}𝐴_𝑚(𝑥)𝑅_𝑖(𝑥)𝐷^𝛼𝑚𝑅_𝑗(𝑥)𝐹𝑑𝑥,_{[𝑒𝜈]𝑖}=−_𝑒𝜈𝑓(𝑥)𝑅_𝑖(𝑥)𝑑𝑥. (3.14) Using the properties of the roof functions, we obtain 𝐾¹_𝑖𝑗=∫0,𝑖−𝑗≥2,_𝑒𝑖𝑑𝑅_𝑖(𝑥)𝐷𝑑𝑥^𝛼𝑅_𝑗(𝑥)∫𝑑𝑥,𝑗=𝑖+1,_𝑒𝑗𝑑𝑅_𝑖(𝑥)𝐷𝑑𝑥^𝛼𝑅_𝑗(𝑥)∫𝑑𝑥,𝑖=𝑗+1,_𝑒𝑖−1+∫_𝑒𝑖𝑑𝑅_𝑖(𝑥)𝐷𝑑𝑥^𝛼𝑅_𝑖(𝑥)𝑑𝑥,𝑗=𝑖. (3.15) Similarly, 𝐾²_𝑖𝑗=0,𝑖−𝑗≥2,_𝑀𝑚∑₌₁∫_𝑒𝑖𝐴_𝑚(𝑥)𝑅_𝑖(𝑥)𝐷^𝛼𝑚𝑅_𝑗(𝑥)𝑑𝑥,𝑗=𝑖+1,_𝑀𝑚∑₌₁∫_𝑒𝑗𝐴_𝑚(𝑥)𝑅_𝑖(𝑥)𝐷^𝛼𝑚𝑅_𝑗(𝑥)∫𝑑𝑥,𝑖=𝑗+1,_𝑒𝑖−1+∫_𝑒𝑖𝐴_𝑚(𝑥)𝑅_𝑖(𝑥)𝐷^𝛼𝑚𝑅_𝑗(𝑥)𝐹𝑑𝑥,𝑗=𝑖,_𝑖=−_𝑒𝑖−1+_𝑒𝑖𝑓(𝑥)𝑅_𝑖(𝑥)𝑑𝑥. (3.16) Note that ∫_𝑒0𝐺(𝑥)𝑑𝑥=0=∫_𝑒𝑛+1𝐺(𝑥)𝑑𝑥, for any function 𝐺(𝑥).

Also, the initial conditions given by (1.2) gives 𝑐₁=𝑏₀,𝑐₂=ℓ𝑏₁+𝑏₀. (3.17) Therefore, the final linear system which gives the unknown 𝑐_𝑖,𝑖=2,3,4,…,𝑁 takes the form ^𝑁+1_𝑗=3𝐾_𝑖𝑗𝑐_𝑗=𝐹^∗_𝑖,𝑖=2,3,4,…,𝑁, (3.18) where 𝐾_𝑖𝑗=𝐾¹_𝑖𝑗+𝐾²_𝑖𝑗,𝐹^∗_𝑖=𝐹_𝑖−𝐾_𝑖1𝑐₁−𝐾_𝑖2𝑐_𝑖2. (3.19)

4. Numerical Experiments

We consider the following Cauchy problem: 𝐷^1.5𝑦(𝑡)+2𝐷¹𝑦(𝑡)+3√𝑡𝐷^0.5𝑦(𝑡)+(1−𝑡)𝑦(𝑡)=𝑓(𝑡),𝑦(0)=𝑦^′(0)=0, (4.1) where 𝑓(𝑡)=2Γ()√1.54𝑡+4𝑡+Γ()𝑡1.5²+(1−𝑡)𝑡². (4.2) It is easy to check that all assumptions of Theorem 2.2 are fulfilled. The exact solution of the fractional differential equation is 𝑦(𝑡)=𝑡².

This problem was solved numerically by the modified Galerkin method on the interval [0,1] using different values of 𝑁, the number of nodal points. In Table 1, some results for different values of the parameters 𝑁 are presented.

Table 1 

Denoting by 𝑒_𝑁=max{|𝑢(𝑥)−𝑢_𝑁(𝑥)|,0<𝑥<1} the errors and by 𝛼_𝑁=log(𝑒_𝑁/𝑒_2𝑁) 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, 𝐾, is a tridiagonal matrix. Moreover, the coefficients of the matrices 𝐾 and 𝐹, 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 𝑚 (with 𝑚 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.