Open access peer-reviewed chapter

# Approximate Solutions of Some Boundary Value Problems by Using Operational Matrices of Bernstein Polynomials

Written By

Submitted: April 26th, 2019 Reviewed: October 28th, 2019 Published: January 27th, 2020

DOI: 10.5772/intechopen.90302

From the Edited Volume

## Functional Calculus

Edited by Kamal Shah and Baver Okutmuştur

Chapter metrics overview

View Full Metrics

## Abstract

In this chapter, we develop an efficient numerical scheme for the solution of boundary value problems of fractional order differential equations as well as their coupled systems by using Bernstein polynomials. On using the mentioned polynomial, we construct operational matrices for both fractional order derivatives and integrations. Also we construct a new matrix for the boundary condition. Based on the suggested method, we convert the considered problem to algebraic equation, which can be easily solved by using Matlab. In the last section, numerical examples are provided to illustrate our main results.

### Keywords

• Bernstein polynomials
• coupled systems
• fractional order differential equations
• operational matrices of integration
• approximate solutions
• 2010 MSC: 34L05
• 65L05
• 65T99
• 34G10

## 1. Introduction

Generalization of classical calculus is known as fractional calculus, which is one of the fastest growing area of research, especially the theory of fractional order differential equations because this area has wide range of applications in real-life problems. Differential equations of fractional order provide an excellent tool for the description of many physical biological phenomena. The said equations play important roles for the description of hereditary characteristics of various materials and genetical problems in biological models as compared with integer order differential equations in the form of mathematical models. Nowadays, most of its applications are found in bio-medical engineering as well as in other scientific and engineering disciplines such as mechanics, chemistry, viscoelasticity, control theory, signal and image processing phenomenon, economics, optimization theory, etc.; for details, we refer the reader to study [1, 2, 3, 4, 5, 6, 7, 8, 9] and the references there in. Due to these important applications of fractional order differential equations, mathematicians are taking interest in the study of these equations because their models are more realistic and practical. In the last decade, many researchers have studied the existence and uniqueness of solutions to boundary value problems and their coupled systems for fractional order differential equations (see [10, 11, 12, 13, 14, 15, 16, 17]). Hence the area devoted to existence theory has been very well explored. However, every fractional differential equation cannot be solved for its analytical solutions easily due to the complex nature of fractional derivative; so, in such a situation, approximate solutions to such a problem is most efficient and helpful. Recently, many methods such as finite difference method, Fourier series method, Adomian decomposition method (ADM), inverse Laplace technique (ILT), variational iteration methods (VIM), fractional transform method (FTM), differential transform method (DTM), homotopy analysis method (HAM), method of radial base function (MRBM), wavelet techniques (WT), spectral methods and many more (for more details, see [9, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38]) have been developed for obtaining numerical solutions of such differential equations. These methods have their own merits and demerits. Some of them provide a very good approximation. However, in the last few years, some operational matrices were constructed to achieve good approximation as in [39]. After this, a variety of operational matrices were developed for different wavelet methods. This method uses operational matrices, where every operation, for example differentiation and integration, involved in these equations is performed with the help of a matrix. A large variety of operational matrices are available in the literature for different orthogonal polynomials like Legendre, Laguerre, Jacobi and Bernstein polynomials [40, 41, 42, 43, 44, 45, 46, 47, 48]. Motivated by the above applications and uses of fractional differential equations, in this chapter, we developed a numerical scheme based on operational matrices via Bernstein polynomials. Our proof is more generalized and there is no need to convert the Bernstein polynomial function vector to another basis like block pulse function or Legendre polynomials. To the best of our knowledge, the method we consider provides a very good approximation to the solution. By the use of these operational matrices, we apply our scheme to a single fractional order differential equation with given boundary conditions as

D α y t + AD μ y t + By t = f t , 1 < α 2 , 0 < μ 1 , y 0 = a , y 1 = b , E1

where f t is the source term, A , B are any real numbers; then we extend our method to solve a boundary value problem of coupled system of fractional order differential equations of the form

D α x t + A 1 D μ 1 x t + B 1 D ν 1 y t + C 1 x t + D 1 y t = f t , 1 < α 2 , 0 < μ 1 , ν 1 1 , D β y t + A 2 D μ 2 x t + B 2 D ν 2 y t + C 2 x t + D 2 y t = g t , 1 < β 2 , 0 < μ 2 , ν 2 1 , y 0 = a , y 1 = b , y 0 = c , y 1 = d , E2

where f t , g t are source terms of the system, A i , B i , C i , D i i = 1 2 are any real constants. Also we compare our approximations to exact values and approximations of other methods like Jacobi polynomial approximations and Haar wavelets methods to evaluate the efficiency of the proposed method. We also provide some examples for the illustration of our main results.

This chapter is designed in five sections. In the first section of the chapter, we have cited some basic works related to the numerical and analytical solutions of differential equations of arbitrary order by various methods. The necessary definitions and results related to fractional calculus and Bernstein polynomials along with the construction of some operational matrices are given in Section 2. In Section 3, we have discussed the main theory for the numerical procedure. Section 4 contains some interesting practical examples and their images. Section 5 describes the conclusion of the chapter.

## 2. Basic definitions and results

In this section, we recall some fundamental definitions and results from the literature, which can be found in [10, 11, 12, 13, 14, 15, 16].

Definition 2.1. The fractional integral of order γ R + of a function y L 1 0 1 R is defined as

I 0 + γ y t = 1 Γ γ 0 t t τ γ 1 y τ .

Definition 2.2. The Caputo fractional order derivative of a function y on the interval 0 1 is defined by

D 0 + γ y t = 1 Γ n γ 0 t t τ n γ 1 y n τ ,

where n = γ + 1 and γ represents the integer part of γ .

Lemma 2.1. The fractional differential equation of order γ > 0

D γ y t = 0 , n 1 < γ n ,

has a unique solution of the form y t = d 0 + d 1 t + d 2 t 2 + + d n 1 t n 1 , where d k R and k = 0,1,2,3 , . , n 1 . .

Lemma 2.2. The following result holds for fractional differential equations

I γ D γ y t = y t + d 0 + c 1 t + d 2 t 2 + + d n 1 t n 1 ,

for arbitrary d k R , k = 0,1,2 , , n 1 .

Hence it follows that

D γ t k = Γ k + 1 Γ k γ + 1 t k γ , I γ t k = Γ k + 1 Γ k + γ + 1 t k + γ and D γ constant = 0 .

### 2.1 The Bernstein polynomials

The Bernstein polynomials B i , m t on 0 1 can be defined as

B i , m t = m i t i 1 t m i , for i = 0,1,2 m ,

where m i = m ! m i ! i ! , which on further simplification can be written in the most simplified form as

B i , m t = k = 0 m i Θ i k m t k + i , i = 0,1,2 m , E3

where

Θ i k m = 1 k m i m i k .

Note that the sum of the Bernstein polynomials converges to 1.

Lemma 2.3. Convergence Analysis: Assume that the function g C m + 1 0 1 that is m + 1 times continuously differentiable function and let X = B 0 , m B 1 , m B m , m . If C T Ψ x is the best approximation of g out of X , then the error bound is presented as

g C T Ψ 2 2 MS 2 m + 3 2 Γ m + 2 2 m + 3 ,

where M = max x 0 1 g m + 1 x , S = max 1 x 0 x 0 .

Proof. In view of Taylor polynomials, we have

F x = g x 0 + x x 0 g 1 x 0 + x x 0 2 Γ 3 g 2 + + x x 0 m Γ m + 1 g m ,

from which we know that

g F x = g m + 1 η x x 0 m + 1 Γ m + 2 , there exist η 0 1 .

Due the best approximation C T Ψ x of g , we have

g C T Ψ x 2 2 g F 2 2 = 0 1 g x F x 2 dx = 0 1 g m + 1 η x x 0 m + 1 Γ m + 2 2 dx M 2 Γ 2 m + 2 0 1 x x 0 2 m + 2 dx 2 M 2 S 2 m + 3 Γ 2 m + 2 2 m + 3 .

Hence we have

g C T Ψ x 2 2 MS 2 m + 3 2 Γ m + 2 2 m + 3 .

Let H = L 2 0 1 be a Hilbert space, then the inner product can be defined as

< f , g > = 0 1 f x . g x dx

and

Y = span B 0 , m B 1 , m B m , m

is a finite dimensional and closed subspace. So if f H is an arbitrary element then its best approximation is unique in Y . This terminology can be achieved by using y 0 Y and for all y Y , we have f y 0 f y . Thus any function can be approximated in terms of Bernstein polynomials as

f t = i = 0 m c i B i m , E4

where coefficient ci can easily be calculated by multiplying (4) by B j m t , j = 0,1,2 , m and integrating over 0 1 by using inner product and d i = 0 1 B i m t f t dt , θ i j = 0 1 B i m t B j m t dt , i , j = 0,1,2 . m , we have

0 1 f t B j m t dt = 0 1 i = 0 m c i B i m t . B j m t dt , j = 0,1,2 m , which implies that 0 1 f t B j m t dt = i = 0 m c i 0 1 B i m t . B j m t dt , j = 0,1,2 m which implies that d 0 d 1 . d m = c 0 c 1 . c m θ 0 0 θ 0 1 θ 0 r θ 0 m θ 1 0 θ 1 1 θ 1 r θ 1 m θ r 0 θ r 1 θ r r θ r m θ m 0 θ m 1 θ m r θ m m . E5

Let X M = d 0 d 1 . d m , C M = c 0 c 1 . c m , where M = m + 1 where M is the scale level and Φ M × M = θ 0 0 θ 0 1 θ 0 r θ 0 m θ 1 0 θ 1 1 θ 1 r θ 1 m θ r 0 θ r 1 θ r r θ r m θ m 0 θ m 1 θ m r θ m m , so

X M = C M Φ M × M C M = X M Φ M × M 1 . E6

where Φ m × M is called the dual matrix of the Bernstein polynomials. After calculating c i , (4) can be written as

f t = C M B M T t , C M is coefficient matrix

where

B M t = B 0 m B 1 m . B m m . E7

Lemma 2.4. Let B M T t be the function vector defined in (3), then the fractional order integration of B M T t is given by

I α B M T t = P M × M α B M T t , E8

where P M × M α is the fractional integration’s operational matrix defined as

P M × M α = P ̂ M × M α Φ M × M 1

and Φ M × M 1 is given in (3) and P M × M α is given by

P ̂ M × M α = Ψ 0 0 Ψ 0 1 Ψ 0 r Ψ 0 m Ψ 1 0 Ψ 1 1 Ψ 1 r Ψ 1 m Ψ r 0 Ψ r 1 Ψ r r Ψ r m Ψ m 0 Ψ m 1 Ψ m r Ψ m m , E9

where

Ψ i , j = k = 0 m i l = 0 m j θ i k m θ j l m Γ k + i + 1 i + j + k + l + α + 1 Γ k + i + α + 1 . E10

Proof. Consider

B i , m t = k = 0 m i θ i k m t k + i E11

taking fractional integration of both sides, we have

I α B i , m t = k = 0 m i θ i k m I α t k + i = k = 0 m i θ i k m Γ k + i + α Γ k + i + α + 1 t k + i + α . E12

Now to approximate right-hand sides of above

k = 0 m i θ i k m Γ k + i + α Γ k + i + α + 1 t k + i + α = C M i B M T t E13

where C M i can be approximated by using (3) as

C M i = X M i Φ M × M 1 , E14

where entries of the vector X M i can be calculated in generalized form as

X M j = 0 1 k = 0 m i θ i k m Γ k + i + α Γ k + i + α + 1 t k + i + α B j , m t dt , j = 0,1,2 . m X M j = 0 1 k = 0 m i θ i k m Γ k + i + α Γ k + i + α + 1 l = 0 m j θ j l m t k + i + α . t l + j dt , j = 0,1,2 . m = k = 0 m i θ i k m l = 0 m j θ j l m Γ k + i + α Γ k + i + α + 1 1 k + l + j + i + α + 1 , j = 0,1,2 , . m E15

evaluating this result for i = 0,1,2....m, we have

I α B 0 , m t I α B 1 , m t I α B m , m t = X M 0 Φ M × M 1 B M T t X M 1 Φ M × M 1 B M T t X M m Φ M × M 1 B M T t ) E16

further writing

Ψ i , j = k = 0 m i l = 0 m j θ i k m . θ j l m Γ k + i + α Γ k + i + α + 1 1 k + l + j + i + α + 1

we get

I α B M T t = P ̂ M × M α . Φ M × M 1 . B M T t . E17

Let us represent

P ̂ M × M α . Φ M × M 1 = P M × M α

thus

I α B M T t = P M × M α . B M T t . E18

Lemma 2.5. Let B M T t be the function vector as defined in (3), then fractional order derivative is defined as

D α B M T t = G M × M α . B M T t E19

where G M × M α is the operational matrix of fractional order derivative given by

G M × M α = G ̂ M × M α Φ M × M 1 , E20

where Φ M × M is the dual matrix given in (3) and

G ̂ M × M α = Ψ 0 0 Ψ 0 1 Ψ 0 r Ψ 0 m Ψ 1 0 Ψ 1 1 Ψ 1 r Ψ 1 m Ψ r 0 Ψ r 1 Ψ r r Ψ r m Ψ m 0 Ψ m 1 Ψ m r Ψ m m , E21

where Ψ i j is defined for two different cases as

Case I: i < α

Ψ i , j = k = α m i l = 0 m j θ i k m . θ j l m Γ k + i α Γ k + i α + 1 1 k + l + j + i α + 1 E22

Case II: i α

Ψ i , j = k = 0 m i l = 0 m j θ i k m . θ j l m Γ k + i α Γ k + i α + 1 1 k + l + j + i α + 1 . E23

Proof. Consider the general element as

D α B i , m t = D α k = 0 m i θ i k m . t k + i = k = 0 m i θ i k m D α t k + i . E24

It is to be noted in the polynomial function B i , m the power of the variable ‘ t ’ is an ascending order and the lowest power is ‘ i ’ therefore the first α 1 terms becomes zero when we take derivative of order α .

Case I: i < α By the use of definition of fractional derivative

D α B i , m t = k = α m i θ i k m Γ k + i + 1 Γ k + i α + 1 t k + i α . E25

Now approximating RHS of (25) as

k = α m i θ i k m Γ k + i + 1 Γ k + i α + 1 t k + i α = C M i B M T t E26

further implies that

X M j = 0 1 k = α m i θ i k m Γ k + i + 1 Γ k + i α + 1 t k + i α B j , m t dt , j = 0,1,2 . m X M j = k = α m i θ i k m l = 0 m j θ j l m Γ k + i + 1 Γ k + i α + 1 k + i + l + j α + 1 , j = 0,1,2 . m E27

Case II: i α if i α then

X M j = k = 0 m i θ i k m l = 0 m j θ j l m Γ k + i + 1 Γ k + i α + 1 k + i + l + j α + 1 , j = 0,1,2 . m . E28

After careful simplification, we get

D α B 0 , m t D α B 1 , m t D α B m , m t = X M 0 Φ M × M 1 B M T t X M 1 Φ M × M 1 B M T t X M m Φ M × M 1 B M T t ) . E29

On further simplification, we have

Ψ i , j = k = α m i l = 0 m j θ i k m . θ j l m Γ k + i + 1 Γ k + i α + 1 1 k + l + j + i α + 1 i < α Ψ i , j = k = 0 m i l = 0 m j θ i k m . θ j l m Γ k + i + 1 Γ k + i α + 1 1 k + l + j + i α + 1 we get D α B M T t = G ̂ M × M α Φ M × M 1 . B M T t . E30

Let

G ̂ M × M α Φ M × M 1 = G M × M α

so

D α B M T t = G M × M α B M T t

which is the desired result. □

Lemma 2.6. An operational matrix corresponding to the boundary condition by taking B M T t is function vector and K is coefficient vector by taking the approximation

u t = K B ̂ t

is given by

Q M × M α , ϕ = Ω 0 0 Ω 0 1 Ω 0 r Ω 0 m Ω 1 0 Ω 1 1 Ω 1 r Ω 1 m Ω r 0 Ω r 1 Ω r r Ω r m Ω m 0 Ω m 1 Ω m r Ω m m , E31

where

Ω i , j = 0 1 Δ i , m ϕ t B j t dt , i , j = 0,1,2 . m .

Proof. Let us take u t = K B ̂ t , then

0 I 1 α K B ̂ t = K 0 I 1 α B ̂ t = K 0 I 1 α B 0 t 0 I 1 α B 1 t 0 I 1 α B m t .

Let us evaluate the general terms

0 I 1 α B i t dt = 1 Γ α 0 1 1 s α 1 B i , m s ds = 1 Γ α k = 0 m i Θ i , k , m 0 1 1 s α 1 s k + i ds . E32

Since by

L 0 1 1 s α 1 s k + i ds = Γ α Γ k + i + 1 τ k + α + i

taking inverse Laplace of both sides, we get

0 1 1 s α 1 . s k + i ds = L 1 Γ α . Γ k + i + 1 τ k + α + i = Γ α . Γ k + i + 1 Γ k + i + α + 1

now Eq. (32) implies that

0 I 1 α B i t dt = k = 0 m i Θ i , k , m Γ k + i + 1 Γ k + i + α + 1 = Δ i , m . E33

Now using the approximation Δ i , m ϕ t = i = 0 m c ̂ i B i t = C M i B M T , and using Eq. (3) we get C M i = K M i Φ M × M 1 B M T and using c j = 0 1 ϕ t B j t dt ,

ϕ t KI α B ̂ t = K Δ 0 , m ϕ t Δ 1 , m ϕ t Δ m , m ϕ t = K C M 0 Φ M × M 1 B M T t C M 1 Φ M × M 1 B M T t C M m Φ M × M 1 B M T t = K c 0 0 c 1 0 c r 0 c m 0 c 0 1 c 1 1 c r 1 c m 1 c 0 r c 1 r c r r c m r c 0 m c 1 m c r m c m m Φ M × M 1 B M T t Φ M × M 1 B M T t Φ M × M 1 B M T t . E34

On further simplification, we get

ϕ t KI α B ̂ t = K Ω 0 0 Ω 0 1 Ω 0 r Ω 0 m Ω 1 0 Ω 1 1 Ω 1 r Ω 1 m Ω r 0 Ω r 1 Ω r r Ω r m Ω m 0 Ω m 1 Ω m r Ω m m B 0 t B 1 t B m t . E35

So

ϕ t 0 I 1 α u t = KQ M × M α , ϕ B M T t ,

and

Ω i , j = 0 1 Δ i , m ϕ t B j t dt , i , j = 0,1,2 . m . E36

which is the required result. □

## 3. Applications of operational matrices

In this section, we are going to approximate a boundary value problem of fractional order differential equation as well as a coupled system of fractional order boundary value problem. The application of obtained operational matrices is shown in the following procedure.

### 3.1 Fractional differential equations

Consider the following problem in generalized form of fractional order differential equation

D α y t + AD μ y t + By t = f t , 1 < α 2 , 0 < μ 1 , subject to the boundary conditions y 0 = a , y 1 = b , E37

where f t is a source term; A , B are any real constants and y t is an unknown solution which we want to determine. To obtain a numerical solution of the above problem in terms of Bernstein polynomials, we proceed as

Let D α y t = K M B M T t . E38

Applying fractional integral of order α we have

y t = K M P M × M α B M T t c 0 + c 1 t

using boundary conditions, we have

c 0 = a c 1 = b a K M P M × M α B M T t t = 1 .

Using the approximation and Lemma 2.2

a + t b a = F M 1 B M T t tP M × M α B M T t t = 1 = Q M × M α , ϕ B M T t .

Hence

y t = K M P M × M α B M T t + a + t b a tK M P M × M α B M T t t = 1 , which gives y t = K M P M × M α B M T t + F M 1 B M T t Q M × M α , ϕ B M T t = K M P M × M α Q M × M α , ϕ B M T t + F M 1 B M T t . E39

Now

D μ y t = D μ K M P M × M α Q M × M α , ϕ B M T t + F M 1 B M T t = K M P M × M α Q M × M α , ϕ G M × M μ B M T t + F M 1 G M × M μ B M T t E40

and

f t = F M 2 B M T t . E41

Putting Eqs. (38)(41) in Eq. (37), we get

K M B M T ( t ) + A K M ( P M × M α Q M × M α , ϕ ) G M × M μ B M T ( t ) + A F M (1) G M × M μ B M T ( t ) + B K M ( P M × M α Q M × M α , ϕ ) B M T ( t ) + B F M (1) B M T ( t ) = F M (2) B M T ( t ). E42

In simple form, we can write (42) as

K M B M T t + AK M P M × M α Q M × M α , ϕ G M × M μ B M T t + AF M 1 G M × M μ B M T t +   BK M P M × M α Q M × M α , ϕ B M T t + BF M 1 B M T t F M 2 B M T t = 0 K M + K M A P ̂ M × M α G M × M μ + B P ̂ M × M α + AF M 1 G M × M μ + BF M 1 F M 2 , E43

where

P ̂ M × M α = P M × M α Q M × M α , ϕ .

Eq. (43) is an algebraic equation of Lyapunov type, which can be easily solved for the unknown coefficient vector KM . When we find KM , then putting this in Eq. (39), we get the required approximate solution of the problem.

### 3.2 Coupled system of boundary value problem of fractional order differential equations

Consider a coupled system of a fractional order boundary value problem

D α x t + A 1 D μ 1 x t + B 1 D ν 1 y t + C 1 x t + D 1 y t = f t , 1 < α 2 , 0 < μ 1 , ν 1 1 , D β y t + A 2 D μ 2 x t + B 2 D ν 2 y t + C 2 x t + D 2 y t = g t , 1 < β 2 , 0 < μ 2 , ν 2 1 , E44

subject to the boundary conditions

x 0 = a , x 1 = b y 0 = c , y 1 = d , E45

where A i , B i , C i , D i i = 1 2 are any real constants, f t , g t are given source terms. We approximate the solution of the above system in terms of Bernstein polynomials such as

D α x t = K M B M T t , D β y t = L M B M T t x t = K M P M × M α B M T t + c 0 + c 1 t , y t = L M ( P M × M β B M T t + d 0 + d 1 t

applying boundary conditions, we have

x t = K M ( P M × M α B M T t + a + t b a tK M P M × M α B M T t t = 1 , y t = K M ( P M × M β B M T t + c + t d c tK M P M × M β B M T t t = 1 .

let us approximate

a + t b a = F M 1 B M T t , c + t d c = F M 2 B M T t tP M × M α B M T t t = 1 = Q M × M α , ϕ B M T t tP M × M β B M T t t = 1 = Q M × M β , ϕ B M T t

then

x t = K M P M × M α B M T t + F M 1 B M T t K M Q M × M α , ϕ B M T t y t = L M P M × M β B M T t + F M 2 B M T t L M Q M × M β , ϕ B M T t D μ 1 x t = K M P M × M α B M T t + F M 1 B M T t K M Q M × M α , ϕ B M T t = K M P M × M α Q M × M α , ϕ G M × M μ 1 + F M 1 G M × M μ 1 B M T t D ν 1 y t = D ν 1 L M P M × M β B M T t + F M 2 B M T t L M Q M × M β , ϕ B M T t = L M P M × M β Q M × M β , ϕ G M × M ν 1 + F M 2 G M × M ν 1 B M T t D μ 2 x t = D μ 2 K M P M × M α B M T t + F M 1 B M T t K M Q M × M α , ϕ B M T t = K M P M × M α Q M × M α , ϕ G M × M μ 2 + F M 1 G M × M μ 2 B M T t

and

D ν 2 y t = D ν 2 K M P M × M β B M T t + F M 2 B M T t K M Q M × M β , ϕ B M T t = L M P M × M β Q M × M β , ϕ G M × M ν 2 + F M 2 G M × M ν 2 B M T t f t = F 3 B M T t , g t = F 4 B M T t .

Thus system (44) implies that

K M B M T t + A 1 K M P M × M α Q M × M α , ϕ G M × M μ 1 + A 1 F M 1 G M × M μ 1 B M T t +   B 1 L M P M × M β Q M × M β , ϕ G M × M ν 1 + B 1 F M 2 G M × M ν 1 B M T t + C 1 K M P M × M α B M T t +   C 1 F M 1 B M T t C 1 K M Q M × M α , ϕ B M T t + D 1 L M P M × M β B M T t + D 1 F M 2 B M T t   D 1 L M Q M × M β , ϕ B M T t = F 3 B M T t L M B M T t + A 2 K M P M × M α Q M × M α , ϕ G M × M μ 2 + A 2 F M 1 G M × M μ 2 B M T t +   B 2 L M P M × M β Q M × M β , ϕ G M × M ν 2 + B 2 F M 2 G M × M ν 2 B M T t + C 2 K M P M × M α B M T t +   C 2 F M 1 B M T t C 2 K M Q M × M α , ϕ B M T t + D 2 L M P M × M β B M T t + D 2 F M 2 B M T t   D 2 L M Q M × M β , ϕ B M T t = F 4 B M T t . E46

Rearranging the terms in the above system and using the following notation for simplicity in Eq. (46)

Q ̂ M × M α = A 1 P M × M α Q M × M α , ϕ G M × M μ 1 + C 1 P M × M α Q M × M α , ϕ Q ̂ M × M β = B 1 P M × M β Q M × M β , ϕ G M × M ν 1 + D 1 P M × M β Q M × M β , ϕ R ̂ M × M α = A 2 P M × M α Q M × M α , ϕ G M × M μ 2 + C 2 P M × M α Q M × M α , ϕ R ̂ M × M β = B 2 P M × M β Q M × M β , ϕ G M × M ν 2 + D 2 P M × M β Q M × M β , ϕ F M = A 1 F M 1 G M × M μ 1 + B 1 F M 2 G M × M ν 1 + C 1 F M 1 + F M 2 D 1 F M 3 G M = A 2 F M 1 G M × M μ 2 + B 2 F M 2 G M × M ν 2 + C 2 F M 1 + D 2 F M 2 F M 4 ,

the above system (46) becomes

K M B M T ( t ) + K M Q ^ M × M α B M T ( t ) + L M Q ^ M × M β B M T ( t ) + F M B M T ( t ) =0 L M B M T ( t ) + K M R ^ M × M α B M T ( t ) + L M R ^ M × M β B M T ( t ) + G M B M T ( t ) =0 [ K M L M ] [ B M T ( t ) 0 0 B M T ( t ) ] + [ K M L M ] [ Q ^ M × M α 0 0 R ^ M × M β ] [ B M T ( t ) 0 0 B M T ( t ) ] + [ K M L M ] [ 0 R ^ M × M α Q ^ M × M β 0 ] [ B M T ( t ) 0 0 B M T ( t ) ] + [ F M G M ] [ B M T ( t ) 0 0 B M T ( t ) ] =0 [ K M L M ] + [ K M L M ] [ Q ^ M × M α R ^ M × M α Q ^ M × M β R ^ M × M β ] + [ F M G M ] =0, E47

which is an algebraic equation that can be easily solved by using Matlab functional solver or Mathematica for unknown matrix K M L M . Calculating the coefficient matrix K M , L M and putting it in equations

x t = K M P M × M α B M T t + F M 1 B M T t K M Q M × M α , ϕ B M T t y t = L M P M × M β B M T t + F M 2 B M T t L M Q M × M β , ϕ B M T t ,

we get the required approximate solution.

## 4. Applications of our method to some examples

Example 4.1. Consider

D α y t + c 1 D ν y t + c 2 y t = f t , 1 < α < 2 E48

subject to the boundary conditions

y 0 = 0 , y 1 = 0 .

Solution: We solve this problem under the following parameters sets defined as

S 1 = α = 2 ν = 1 c 1 = 1 c 2 = 1 , S 2 = α = 1.8 ν = 0.8 c 1 = 10 c 2 = 100 , S 3 = α = 1.5 ν = 0.5 c 1 = 1 / 10 c 2 = 1 / 100 , and select source term for S 1 as

f 1 t = t 6 t 1 3 + t 6 72 t 168 + 126 30 t 4 + 3 t 5 3 t 2 t 1 2 E49
f 2 ( t ) = 11147682583723703125 x 21 5 ( 1750 x 3 4200 x 2 + 3255 x 806 ) 406548945561989414912 + 278692064593092578125 x 26 5 ( 5250 x 3 14350 x 2 + 12915 x 3813 ) 25002760152062349017088 + 100 x 6 ( x 1 ) 3 , E50
f 3 ( t ) = 5081767996463981 x 9 2 ( 1344 x 3 3360 x 2 + 2730 x 715 ) 264146673456906240 + 5081767996463981 x 11 2 ( 1344 x 3 3808 x 2 + 3570 x 1105 ) 22452467243837030400 + x 6 ( x 1 ) 3 100 . E51

The exact solution of the above problem is

y t = t 6 t 1 3 .

We solve this problem with the proposed method under different sets of parameters as defined in S 1 , S 2 , S 3 . The observation and simulation demonstrate that the solution obtained with the proposed method is highly accurate. The comparison of exact solution with approximate solution obtained using the parameters set S 1 is displayed in Figure 1 subplot (a), while in Figure 1 subplot (b) we plot the absolute difference between the exact and approximate solutions using different scale levels. One can easily observe that the absolute error is much less than 10 12 . The order of derivatives in this set is an integer.

By solving the problem under parameters set S 2 and S 3, we observe the same phenomena. The approximate solution matches very well with the exact solution. See Figures 2 and 3 respectively.

Example 4.2. Consider

D α y t 2 D 0.9 y t 3 y t = 4 cos 2 t 7 sin 2 t E52

subject to the boundary conditions

y 0 = 0 , y 1 = sin 2 .

Solution: The exact solution of the above problem is y t = sin 2 t , when α = 2 . However the exact solution at fractional order is not known. We use the well-known property of FDEs that when α 2 , the approximate solution approaches the exact solution for the evaluations of approximate solutions and check the accuracy by using different scale levels. By increasing the scale level M, the accuracy is also increased. By the proposed method, the graph of exact and approximate solutions for different values of M and at α = 1.7 is shown in Figure 4 . From the plot, we observe that the approximate solution becomes equal to the exact solution at α = 2 . We approximate the error of the method at different scale levels and record that when scale level increases the absolute error decreases as shown in Figure 4 subplot (b) and accuracy approaches 10 −9, which is a highly acceptable figure. For convergence of our proposed method, we examined the quantity 0 1 y exact y approx dt for different values of M and observed that the norm of error decreases with a high speed with the increase of scale level M as shown in Figure 4b .

Example 4.3. Consider the following coupled system of fractional differential equations

D 1.8 x t + Dx t + 9 D 0.8 y t + 2 x t y t = f t D 1.8 y t 6 D 0.8 x t + Dy t x t = g t E53

subject to the boundary conditions

x 0 = 1 , x 1 = 2 and y 0 = 2 , y 1 = 2 .

Solution: The exact solution is

x t = t 5 1 t , y t = t 4 1 t .

We approximate the solution of this problem with this new method. The source terms are given by

f ( t ) =2 x 5 ( x 1 ) x 4 ( x 1 ) + x 4 ( 6 x 5 ) 2229536516744740625 x 16 5 ( 10 x 7 ) 1008806316530991104 + 1337721910046844375 x 16 5 ( 25 x 21 ) 1008806316530991104 E54
g ( t ) = x 3 ( 5 x 4 ) x 5 ( x 1 ) 11147682583723703125 x 21 5 ( 15 x 13 ) 6557241057451442176 89181460669789625 x 11 5 ( 25 x 16 ) 144115188075855872 . E55

In the given Figure 5 , we have shown the comparison of exact x t , y t and approximate x t , y t in subplot (a) and (b) respectively.

As expected, the method provides a very good approximation to the solution of the problem. At first, we approximate the solutions of the problem at α = 2 because the exact solution at α = 2 is known. We observe that at very small scale levels, the method provides a very good approximation to the solution. We approximate the absolute error by the formula

X error = x exact x approx .

and

Y error = y exact y approx .

We approximate the absolute error at different scale level of M, and observe that the absolute error is much less than 10 10 at scale level M = 7 , see Figure 6 . We also approximate the solution at some fractional value of α and observe that as α 2 the approximate solution approaches the exact solution, which guarantees the accuracy of the solution at fractional value of α . Figure 6 shows this phenomenon. In Figure 6 , the subplot (a) represents the absolute error of x t and subplot (b) represents the absolute error of y t .

Example 4.4. Consider the following coupled system

D 1.8 x t x t + 3 y t = f t D 1.8 y t + 4 x t 2 y t = g t , E56

subject to the boundary conditions

x 0 = 1 , x 1 = 1 and y 0 = 1 , y 1 = 1 .

Solution: The exact solution for α = β = 2 is

x t = t 5 t 4 1 , and y t = t 4 t 3 2 . E57

The source terms are given by

f t = 445907303348948125 t 3.2 25 t 21 3026418949592973312 + t 3 4 t 4 + 3 t 5 2 , g t = 89181460669789625 t 2.5 5 t 4 14411518807585872 4 t 3 + 6 t 4 2 t 5 2 .

Approximating the solution with the proposed method, we observe that our scheme gives high accuracy of approximate solution. In Figure 7 , we plot the exact solutions together with the approximate solutions in Figure 7(a) and (b) for x t and y t , respectively. We see from the subplots (a) and (b) that our approximations have close agreement to that of exact solutions. This accuracy may be made better by increasing scale level. Further, one can observe that absolute error is below 10 10 in Figure 8 , which indicates better accuracy of our proposed method for such types of practical problems of applied sciences.

In Figure 8 , the subplot (a) represents absolute error for x t while subplot (b) represents the same quantity for y t . From the subplots, we see that maximum absolute error for our proposed method for the given problem (4.4) is below 10 10 . This is very small and justifies the efficiency of our constructed method.

Example 4.5. Consider the boundary value problem

D α x t + ωπ 2 D ν x t + x t = ωπ sin ωπt + ωπ x 0 = 0 , x 1 = 2 . E58

Taking α = 2 , ν = 1 and ω = 1,3,5 , , the exact solution is given by

x t = cos ωπt 1 .

We plot the comparison between exact and approximate solutions to the given example at M = 10 and corresponding to ω = 3.5 , α = 2 , β = 1 . Further, we approximate the solution through Legendre wavelet method (LWM) [47], Jacobi polynomial method (JM) and Bernstein polynomials method (BM), as shown in Figure 9 .

From Table 1 , we see that Bernstein polynomials also provide excellent solutions to fractional differential equations [48].

ω M α ν x app x ex at BM x app x ex at WM x app x ex at JM
0.5 10 2 1 7.000 3 2.966 1 1.500 2
1.5 15 1.6 0.9 6.091 3 4.918 2 1.623 1
2.0 20 1.8 0.8 1.237 3 2.108 2 2.723 2
3.5 25 1.9 0.7 1.008 3 5.795 2 1.813 3

### Table 1.

Comparison of solution between Legendre wavelet method (LWM) [47], Jacobi polynomial method (JM) and Bernstein polynomials method (BM) for Example 4.5.

## 5. Conclusion and future work

The above analysis and discussion take us to the conclusion that the new method is very efficient for the solution of boundary value problems as well as initial value problems including coupled systems of fractional differential equations. One can easily extend the method for obtaining the solution of such types of problems with other kinds of boundary and initial conditions. Bernstein polynomials also give best approximate solutions to fractional order differential equations like Legendre wavelet method (LWM), approximation by Jacobi polynomial method (JPM), etc. The new operational matrices obtained in this method can easily be extended to two-dimensional and higher dimensional cases, which will help in the solution of fractional order partial differential equations. Also, we compare our result to that of approximate methods for different scale levels. We observed that the proposed method is also an accurate technique to handle numerical solutions.

## Acknowledgments

This research work has been supported by Higher Education Department (HED) of Khyber Pakhtankhwa Government under grant No: HEREF-46 and Higher Education Commission of Pakistan under grant No: 10039.

## References

1. 1. Lakshmikantham V, Leela S, Vasundhara J. Theory of Fractional Dynamic Systems. Cambridge, UK: Cambridge Academic Publishers; 2009
2. 2. EL-Sayed AMA, Bin-Taher EO. Positive solutions for a nonlocal multi-point boundary-value problem of fractional and second order. Electronic Journal of Differential Equations. 2013;64:1-8
3. 3. Magin RL. Fractional calculus in bioengineering—Part 2. Journal of Critical Reviews in Biomedical Engineering. 2004;32(2):105-193
4. 4. Magin RL. Fractional calculus in bioengineering. Journal of Critical Reviews in Biomedical Engineering. 2004;32(1):1-104
5. 5. Rossikhin YA, Shitikova MV. Applications of fractional calculus to dynamic problems of linear and nonlinear hereditary mechanics of solids. Journal of Applied Mechanics Reviews. 1997;50:15-67
6. 6. Ichise M, Nagayanagi Y, Kojima T. An analog simulation of non-integer order transfer functions for analysis of electrode processes. Journal of Electroanalytical Chemistry and Interfacial Electrochemistry. 1971;33:253-265
7. 7. Inca M, Kiliça B. Classification of traveling wave solutions for time-fractional fifth-order KdV-like equation. Waves in Random and Complex Media. 2014:393-404
8. 8. He JH. Some applications of nonlinear fractional differential equations and their approximations. Bulletin of Science and Technology. 1999;15(2):86-90
9. 9. Rehman M, Khan RA. The Legender wavelet method for solving fractional differential equation. Communications in Nonlinear Science and Numerical Simulation. 2011;10:4163-4173
10. 10. Kilbas AA, Marichev OI, Samko SG. Fractional Integrals and Derivatives (Theory and Applications). Switzerland: Gordon and Breach; 1993
11. 11. Miller KS, Ross B. An Introduction to the Fractional Calculus and Fractional Differential Equations. New York: Wiley; 1993
12. 12. Podlubny I. Fractional Differential Equations, Mathematics in Science and Engineering. New York: Academic Press; 1999
13. 13. Hilfer R. Applications of Fractional Calculus in Physics. Singapore: World Scientific; 2000
14. 14. Kilbas AA, Srivastava HM, Trujillo JJ. Theory and Applications of Fractional Differential Equations, North-Holland Mathematics Studies. Vol. 204. Amsterdam: Elsevier; 2006
15. 15. Benchohra M, R. Graef J, Hamani S. Existence results for boundary value problems with nonlinear fractional differential equations. Journal of Applied Analysis. 2008;87:851-863
16. 16. Rehman M, Khan R. A note on boundary value problems for a coupled system of fractional differential equations. Computers & Mathematcs with Applications. 2011;61:2630-2637
17. 17. Su X. Boundary value problem for a coupled system of nonlinear fractional differential equations. Applied Mathematics Letters. 2009;22:64-69
18. 18. Idrees Bhattia M, Bracken P. Solutions of differential equations in a Bernstein polynomial basis. Journal of Computational and Applied Mathematics. 2007;205:272-280
19. 19. Henryk G, Jose’ LP. On the approximation properties of Bernstein polynomials via probabilistic tools. Boletin de la Asociacion Matematica Venezolana. 2003;X(1):1
20. 20. Ordokhani Y, Davaei S. Approximate solutions of differential equations by using the Bernstein polynomials. Applied Mathematics. 2011. DOI: 10.5402/2011/787694. Article ID 787694
21. 21. Mandal BN, Bhattacharya S. Numerical solution of some classes of integral equations using Bernstein polynomials. Applied Mathematics and Computation. 2007;190(2):1707-1716
22. 22. Yousefi SA, Behroozifar M. Operational matrices of Bernstein polynomials and their applications. International Journal of Systems Science. 2010;41(6):709-716
23. 23. Blank L. Numerical treatment of differential equations of fractional order. In: Numerical Analysis Report 287. Manchester Centre for Computational Mathematics; 1996
24. 24. Sun HH, Onaral B, Tsao Y. Application of positive reality principle to metal electrode linear polarization phenomena. IEEE Transactions on Biomedical Engineering. 1984;31(10):664-674
25. 25. Baillie RT. Long memory processes and fractional integration in econometrics. Journal of Econometrics. 1996;73:5-59
26. 26. Momani S, Odibat Z. Numerical approach to differential equations of fractional order. Journal of Computational and Applied Mathematics. 2007;207:96-110
27. 27. El-Wakil SA, Elhanbaly A, Abdou MA. Adomian decomposition method for solving fractional nonlinear differential equations. Applied Mathematics and Computation. 2006;182:313-324
28. 28. Odibat Z, Momani S. Numerical methods for nonlinear partial differential equations of fractional order. Applied Mathematical Modelling. 2008;32:28-39
29. 29. Erturk VS, Momani S. Solving systems of fractional differential equations using differential transform method. Journal of Computational and Applied Mathematics. 2008;215:142-151
30. 30. Khalil H, Khan RA. A new method based on legender polynomials for solution of system of fractional order partial differential equation. International Journal of Computer Mathematics. 2014;91(12):2554-2567
31. 31. Sweilam NH, Khader MM, Al-Bar RF. Numerical studies for a multi-order fractional differential equation. Physics Letters A. 2007;371:26-33
32. 32. Das S. Analytical solution of a fractional diffusion equation by variational iteration method. Computers & Mathematcs with Applications. 2009;57:483-487
33. 33. Momani S, Odibat Z. Homotopy perturbation method for nonlinear partial differential equations of fractional order. Physics Letters A. 2007;365:345-350
34. 34. Hashim I, Abdulaziz O, Momani S. Homotopy analysis method for fractional IVPs. Communications in Nonlinear Science and Numerical Simulation. 2009;14:674-684
35. 35. Arikoglu A, Ozkol I. Solution of fractional differential equations by using differential transform method. Chaos, Solitons & Fractals. 2007;34:1473-1481
36. 36. Arikoglu A, Ozkol I. Solution of fractional integro-differential equations by using fractional differential transform method. Chaos, Solitons & Fractals. 2009;40:521-529
37. 37. Erturk VS, Momani S, Odibat Z. Application of generalized differential transform method to multi-order fractional differential equations. Communications in Nonlinear Science and Numerical Simulation. 2008;13:1642-1654
38. 38. Odibat Z, Shawagfeh N. Generalized Taylor’s formula. Applied Mathematics and Computation. 2007;186:286-293
39. 39. Saadatmandi A, Deghan M. A new operational matrix for solving fractional-order differential equation. Computers & Mathematics with Applications. 2010;59:1326-1336
40. 40. Kumar P, Agrawal OP. An approximate method for numerical solution of fractional differential equations. Signal Processing. 2006;86:2602-2610
41. 41. Liua F, Anh V, Turner I. Numerical solution of the space fractional Fokker-Planck equation. Journal of Computational and Applied Mathematics. 2004;146:209-219
42. 42. Yang W. Positive solutions for a coupled system of nonlinear fractional differential equations with integral boundary conditions. Computers & Mathematcs with Applications. 2012;63:288-297
43. 43. Diethelm K. An algorithm for the numerical solution of differential equations of fractional order. Electronic Transactions on Numerical Analysis. 1997;5:1-10
44. 44. Diethelm K, Walz G. Numerical solution of fractional order differential equations by extropolation. Numerical Algorithms. 1997;16:231-253
45. 45. Li Y, Shah K. Numerical solutions of coupled systems of fractional order partial differential equations. Advances in Mathematical Physics. 2017;2017:14
46. 46. Shah K, Wang J. A numerical scheme based on nondiscretization of data for boundary value problems of fractional order differential equations. RACSAM. 2018;2018:1-18
47. 47. Rehman M. Boundary value problems for fractional differential equations: Existence theory and numerical solutions [PhD dissertation]. Islamabad, Pakistan: NUST University; 2016
48. 48. Darania P, Ebadian A. A method for the numerical solution of the integro-differential equations. Applied Mathematics and Computation. 2007;188:657-668

Written By