Open AccessArticle

Extension of Operational Matrix Technique for the Solution of Nonlinear System of Caputo Fractional Differential Equations Subjected to Integral Type Boundary Constrains

Department of Mathematics, Attock Campus, University of Education, Lahore 43600, Pakistan

Department of Basic Sciences, University of Engineering and Technology Peshawar, Peshawar 25000, Pakistan

School of Mathematical Sciences, Universiti Kebangsaan Malaysia, Bangi Selangor 43600, Malaysia

⁴

Department of Mathematics, Anand International College of Engineering, Jaipur 303012, India

Author to whom correspondence should be addressed.

Entropy 2021, 23(9), 1154; https://doi.org/10.3390/e23091154

Submission received: 31 May 2021 / Revised: 12 August 2021 / Accepted: 18 August 2021 / Published: 2 September 2021

(This article belongs to the Special Issue Dynamical Systems, Differential Equations and Applications)

Download

Browse Figures

Figure 1
(a) Comparison of exact and approximate solution at different scale levels of Test Problem 1. (b) Absolute difference in the exact and approximate solutions at different scale levels of Test Problem 1. "> Figure 2
Errornorm for different values of <math display="inline"><semantics> <msup> <mi>σ</mi> <mn>1</mn> </msup> </semantics></math> and <math display="inline"><semantics> <msup> <mi>σ</mi> <mn>2</mn> </msup> </semantics></math> in Test Problem 2. "> Figure 3
Errornorm for different values of <math display="inline"><semantics> <msup> <mi>σ</mi> <mn>1</mn> </msup> </semantics></math> and <math display="inline"><semantics> <msup> <mi>σ</mi> <mn>2</mn> </msup> </semantics></math> in Test Problem 3. "> Figure 4
Error norm on boundary for different values of <math display="inline"><semantics> <msup> <mi>σ</mi> <mn>1</mn> </msup> </semantics></math> and <math display="inline"><semantics> <msup> <mi>σ</mi> <mn>2</mn> </msup> </semantics></math> in Test Problem 4. "> Figure 5
Error norm of Test Problem 4. "> Figure 6
Error norm of Test Problem 5. "> Figure 7
Error norm of Test Problem 6. ">

Review Reports Versions Notes

Abstract

We extend the operational matrices technique to design a spectral solution of nonlinear fractional differential equations (FDEs). The derivative is considered in the Caputo sense. The coupled system of two FDEs is considered, subjected to more generalized integral type conditions. The basis of our approach is the most simple orthogonal polynomials. Several new matrices are derived that have strong applications in the development of computational scheme. The scheme presented in this article is able to convert nonlinear coupled system of FDEs to an equivalent S-lvester type algebraic equation. The solution of the algebraic structure is constructed by converting the system into a complex Schur form. After conversion, the solution of the resultant triangular system is obtained and transformed back to construct the solution of algebraic structure. The solution of the matrix equation is used to construct the solution of the related nonlinear system of FDEs. The convergence of the proposed method is investigated analytically and verified experimentally through a wide variety of test problems.

Keywords:

approximation; numerical simulation; iterative methods

1. Introduction

Fractional calculus has a long and rich history. Its wide range of applications made this field an active area of mathematical research. Frequently investigated properties of FDEs include existence and uniqueness problems, multiplicity of solutions, stability of solutions and analytical study of the analytical properties of solution. Parallel to these area of research, one of the most explored and interested areas of research in this field is the design of new numerical methods for finding the approximate solutions to problems of this category. Many scientists and mathematicians are trying to design efficient and reliable techniques to find possible estimates to solutions of FDEs or their coupled systems.

There are many analytical, semi-analytical, numerical, and spectral approximations of solution to FDEs and their coupled systems. Among others, one of the easiest method is the differential transformation method (DTM). In [1], DTM is successfully applied to solve simple nonlinear FDEs with simple initial conditions. In [2], DTM is designed to solve the fractional-order counterpart of Korteweg De Vries (KDV) equation. The method is further improved for the solution of fractional-order boundary value problems in [3]. Solutions to fractional-order boundary value problems are also attempted with analytical methods such as the homotopy perturbation method; see for example [4,5,6]. The Adomian decomposition method is also a powerful analytical method [7,8,9]. Spectral methods have gained the attention of scholars in recent decades. Compared to other methods, spectral methods are easy to design and implement [10,11,12,13,14,15,16,17,18,19,20]. The availability of a wide range of orthogonal polynomials makes this method more interesting. They have the ability to solve fractional order problems, whose solutions are difficult or sometimes impossible to obtain with other traditional methods. For new readers, we strongly recommend studying the results obtained in [21,22,23,24,25,26] for a clear understanding and developing a good base. However, to the best of our knowledge, the spectral method becomes difficult and sometimes fails to handle the situation when boundary conditions are given in more complicated forms such as local conditions, nonlocal m-point terminal conditions, integral type terminal conditions, and radiation boundary conditions. Such conditions have solid application in various problems of traveling waves, heat conduction and electromagnetism. One can find a good example of application of integral type boundary condition to heat conduction phenomena in a rod of fixed length in a recent article [27].

Nonlocal FDEs arise in mathematical modeling of various problems in physics, engineering, ecology, and biological sciences [28,29,30]. Some of the numerical investigations regarding FDEs with nonlocal constrains are discussed in [31,32,33,34,35]. Numerical approaches such as finite difference and radial base function also remain a focus of interest. Application of these methods to one-dimensional heat-like equations has been studied in [32,36,37,38]. Two-dimensional diffusion problems [33,39,40] and Laplace equations with integral constraints are explored in [31].

Keeping in view the increasing interest of mathematicians in fractional calculus, we are strongly motivated to design a new spectral approximation technique for complicated problems such as

\begin{matrix} ^{c} D^{σ_{1}} u (t) = f (u, v, u^{(γ_{1})}, v^{(γ_{2})}), \\ ^{c} D^{σ_{2}} v (t) = g (u, v, u^{(γ_{1})}, v^{(γ_{2})}), \\ u (0) = u_{0}, u (τ) = m_{1} \int_{0}^{ω_{1}} s (t) u (t) d t, 0 < ω_{1} \leq τ, \\ v (0) = v_{0}, v (τ) = m_{2} \int_{0}^{ω_{2}} r (t) v (t) d t, 0 < ω_{2} \leq τ . \end{matrix}

(1)

where

0 < γ_{1}, γ_{2} \leq 1,

1 < σ_{1}, σ_{2} \leq 2,

t \in [0, τ] .

The scalar functions

u (t)

and

v (t)

are the solution to be determined. f and g are nonlinear functions of

u (t),

v (t)

and its fractional order derivatives and is assumed to be in such a form that the solution of the problem exists.

We start our discussion by introducing some definitions and preliminary results.

2. Preliminaries

The following definitions and notations are important for our further analysis. More details and theoretical understanding of these results, see [41,42,43,44,45].

Definition 1.

The Riemann–Liouville

α -

order integral of

ϕ \in (L^{1} [a, b], R)

is defined by the following relation.

I^{α} ϕ (t) = \frac{1}{Γ (α)} \int_{a}^{t} {(t - s)}^{α - 1} ϕ (s) d s,

Definition 2.

For

ϕ (t) \in C^{n} [a, b]

, the α order Caputo derivative is defined as

\begin{matrix} ^{c} D^{α} ϕ (t) = & \frac{1}{Γ (n - α)} \int_{a}^{t} \frac{ϕ^{(n)} (s)}{{(t - s)}^{α + 1 - n}} d s, \\ n - 1 \leq α < n, n \in N, \end{matrix}

where

n = [α] + 1

From the above definition, it is easy to extract

I^{α} t^{b} = \frac{Γ (1 + b)}{Γ (1 + b + α)} t^{b + α} for α > 0, k \geq 0

^{c} D^{α} C = 0,

and

^{c} D^{α} t^{b} = \frac{Γ (1 + b)}{Γ (1 + b - α)} t^{b - α}, for b \geq [α] .

(2)

The Shifted Legendre Polynomials (LP)

These polynomials play a central role in approximation theory. Generally, these polynomials are defined on the domain

[0, τ]

, which is given by

ρ_{l}^{τ} (t) = \sum_{b = 0}^{l} ℷ_{(l, b)} t^{b},

(3)

where

ℷ_{(l, b)} = \frac{{(- 1)}^{l + b} (l + b)!}{(l - b)! τ^{b} {(b!)}^{2}} .

(4)

These polynomials enjoys a very important property of orthogonality on the domain

[0, τ],

which is expressed mathematically as

\int_{0}^{τ} ρ_{i}^{τ} (t) ρ_{j}^{τ} (t) d t = \{\begin{matrix} \frac{τ}{2 j + 1}, & if i = j, \\ 0, & if i \neq j . \end{matrix}

(5)

By using Equation (5), any

s (t)

can be expressed in terms of LP as

s (t) \approx \sum_{l = 0}^{m} c_{l} ρ_{l}^{τ} (t), where c_{l} = \frac{(2 l + 1)}{τ} \int_{0}^{τ} s (t) ρ_{l}^{τ} (t) d t .

(6)

The above equation has an equivalent vector representation given as

s (t) \approx C_{M} Λ_{M}^{τ} (t),

(7)

where

Λ_{M}^{τ} (t) = {[\begin{matrix} ρ_{0}^{τ} (t) & ρ_{1}^{τ} (t) & \dots & ρ_{i}^{τ} (t) & \dots & ρ_{m}^{τ} (t) \end{matrix}]}^{T}

(8)

and

C_{M} = [\begin{matrix} c_{0} & c_{1} & \dots & c_{i} & \dots & c_{m} \end{matrix}] .

(9)

The following useful constant is important in the derivation of the operational matrices. We only recall the definition of the constant. The detailed derivation of which can be found in [25].

Lemma 1.

The integral of product of any three LP is a constant, represented by,

℘^{(i, j, k)},

defined as

\int_{0}^{τ} ρ_{i}^{τ} (t) ρ_{j}^{τ} (t) ρ_{k}^{τ} (t) d t = ℘^{(i, j, k)},

where

℘^{(i, j, k)} = \sum_{a = 0}^{i} \sum_{p = 0}^{j} \sum_{r = 0}^{k} ℷ_{(i, a)} ℷ_{(i, p)} ℷ_{(i, r)} Y_{(a, p, r)} .

ℷ_{(., .)}

are as defined in (4) and

Y_{(a, p, r)} = \frac{τ^{(a + p + r + 1)}}{(a + p + r + 1)} .

The constant defined above is important in the solution of FDEs. We recall one more important result of the Legendre polynomials, which is their application in the study of convergence and developing of error bounds.

Theorem 1

([21]). Let

\prod_{M}

be the space of M terms Legendre polynomials and let

u (t) \in C^{m} [0, 1],

then

u_{m} (t)

is in space

\prod_{M}

; then, for the m term approximation,

u (t) = \sum_{i = 0}^{m} c_{i} ρ_{i}^{τ} (t),

there exists a constant C such that

c_{k} ≃ \frac{C}{λ_{k}} | u^{(m)} | .

and

| u (t) - \sum_{i = 0}^{m} c_{i} ρ_{i}^{τ} {(t) |}^{2} \leq \sum_{k = m + 1}^{\infty} λ_{k} c_{k}^{2},

where

c_{k} = \frac{λ_{k} + 1}{τ} \int_{0}^{τ} u (t) ρ_{k}^{τ} (t) d t, λ_{k} = k (k + 1) .

C is constant and can be chosen in such a way that

u^{(2 m)}

belongs to

\prod_{M}

, where

u^{(m)}

is defined as

u^{(m)} = Z (u^{(m - 1)}) = Z^{m} (u^{(0)})

where Z is storm livoliel operator of legendre polynomials with

u^{(0)} = u (t) .

In the next section, we first recall some of our previously designed operational matrices and then develop new operational matrices.

3. Operational Matrices (OP)

OP acts as a basic block in developing approximation techniques. The purpose of operational matrices is to replace a given derivative term with its matrix notation. The following matrices are important in our further investigation.

Lemma 2

([24]). Let

Λ_{M}^{τ} (t)

be the function vector; the α order integration is generalized as

I^{α} (Λ_{M}^{τ} (t)) ≃ H_{M \times M}^{τ, α} Λ_{M}^{τ} (t),

where

H_{M \times M}^{τ, α}

is the OP of integration, defined as

H_{M \times M}^{τ, α} = [\begin{matrix} Θ_{i, j, τ} \end{matrix}],

where

Θ_{i, j, τ} = \sum_{a = 0}^{i} s_{(a, j)} ℷ_{(i, a)} \frac{Γ (a + 1)}{Γ (a + α + 1)} .

where

s_{(a, j)} = \frac{(2 j + 1)}{τ} \sum_{l = 0}^{j} \frac{{(- 1)}^{j + l} (j + l)! {(τ)}^{a + l + α + 1}}{(τ^{l}) (j - l) {(l!)}^{2} (a + l + α + 1)} .

Corollary 1.

The error

| E_{M} | = | I^{α} u (t) - C_{M} H_{M \times M}^{τ, α} Λ_{M}^{τ} (t) |

is bounded by the following relation

| E_{M} | \leq | \sum_{a = m + 1}^{\infty} c_{k} {\sum_{i = 0}^{m} Θ_{a, j, τ}} | .

Proof.

Consider

u (t) = \sum_{a = 0}^{\infty} c_{a} ρ_{a}^{τ} (t) .

Then, using the previous result, we get

I^{α} u (t) = \sum_{a = 0}^{\infty} c_{a} \sum_{j = 0}^{m} Θ_{a, j, τ} ρ_{j}^{τ} (t) .

Simplifying the above estimate

I^{α} u (t) - \sum_{a = 0}^{m} c_{a} \sum_{j = 0}^{m} Θ_{a, j, τ} ρ_{j}^{τ} (t) = \sum_{a = m + 1}^{\infty} c_{a} \sum_{j = 0}^{m} Θ_{a, j, τ} ρ_{j}^{τ} (t) .

In matrix notation

I^{α} u (t) - C_{M} H_{M \times M}^{τ, α} Λ_{M}^{τ} (t) = \sum_{a = m + 1}^{\infty} c_{a} \sum_{j = 0}^{m} Θ_{a, j, τ} ρ_{j}^{τ} (t) .

Using the fact

| ρ_{j}^{τ} (t) | \leq 1

for

t \in [0, τ]

, therefore, we can write

| I^{α} u (t) - C_{M} H_{M \times M}^{τ, α} Λ_{M}^{τ} (t) | \leq | \sum_{a = m + 1}^{\infty} c_{a} \sum_{j = 0}^{m} Θ_{a, j, τ} | .

□

Lemma 3

([24]). Let

Λ_{M}^{τ} (t)

be the function vector as defined in (8); then the derivative of order σ of

Λ_{M}^{τ} (t)

is generalized as

^{c} D^{σ} (Λ_{M}^{τ} (t)) ≃ G_{M \times M}^{τ, σ} Λ_{M}^{τ} (t),

where

G_{M \times M}^{τ, σ}

is the operational matrix of derivative of order σ, and

G_{M \times M}^{τ, σ}

is defined as

G_{M \times M}^{τ, σ} = [\begin{matrix} Φ_{i, j, τ} \end{matrix}]

where

Φ_{i, j, τ} = \sum_{k = ⌈ σ ⌉}^{i} s_{(k, j)} ℷ_{(i, k)} \frac{Γ (k + 1)}{Γ (k - σ + 1)}

with

Φ_{i, j, τ} = 0

i < ⌈ σ ⌉

Furthermore,

ℷ_{(i, k)}

is similar as defined in (4) and

s_{(k, j)} = \frac{(2 j + 1)}{τ} \sum_{l = 0}^{j} \frac{{(- 1)}^{j + l} (j + l)! {(τ)}^{k + l - σ + 1}}{(τ^{l}) (j - l) {(l!)}^{2} (k + l - σ + 1)} .

Corollary 2.

The error

| E_{M} {| = |}^{c} D^{σ} u (t) - C_{M} G_{M \times M}^{τ, σ} Λ_{M}^{τ} (t) |

in approximating

D^{σ} u (t)

with operational matrix of derivative is bounded by the following.

| E_{M} | \leq | \sum_{k = m + 1}^{\infty} u_{k} {\sum_{i = ⌈ σ ⌉}^{m} Φ_{i, j, τ}} |

Proof.

The proof of this corollary is similar as Corollary 1. □

Lemma 4

([24]). Let

u (t)

and

ϕ_{n} (t)

be smooth functions that are well-defined on

[0, τ]

. Then

ϕ_{n} (t)^{c} D^{σ} u (t) = W_{M} B_{ϕ_{n}}^{σ} Λ_{M}^{τ} (t) .

where

W_{M}

is the Legendre coefficients vector of

u (t)

as defined in (7) and

B_{ϕ_{n}}^{σ} = G_{M \times M}^{τ, σ} J_{M \times M}^{τ, ϕ_{n}} .

B_{ϕ_{n}}^{σ} = [\begin{matrix} \overset{⏞}{Θ_{(r, j)}} \end{matrix}]

where

\overset{⏞}{Θ_{(i, j)}} = \sum_{l = 0}^{m} Φ_{(i, l)} Ω_{(l, j)}

The matrix

G_{M \times M}^{τ, σ}

is the operational matrix of derivative as defined in Lemma 3; the entries of matrix

J_{M \times M}^{τ, ϕ_{n}}

are defined by the following relation

Ω_{(l, j)} = \frac{2 j + 1}{τ} \sum_{i = 0}^{m} c_{i} ℘_{}^{(i, l, j)} .

and

c_{i} = \int_{0}^{τ} ϕ_{n} (t) ρ_{i} (t) d t

Corollary 3.

The error

| E_{M} | = | ϕ_{n} {(t)}^{c} D^{σ} u (t) - C_{M} B_{ϕ_{n}}^{σ} Λ_{M}^{τ} (t) |

in approximating

ϕ_{n} {(t)}^{c} D^{σ} u (t)

with operational matrix of fractional derivative with variable coefficient is bounded by the following.

| E_{M} | \leq | \sum_{k = m + 1}^{\infty} c_{k} \sum_{j = 0}^{m} \overset{⏞}{Θ_{(k, j)}} | .

Proof.

Consider

u (t) = \sum_{k = 0}^{\infty} c_{k} ρ_{k}^{τ} (t) .

Then, using the relation above, we get

ϕ_{n} {(t)}^{c} D^{σ} u (t) = \sum_{k = 0}^{\infty} c_{k} \sum_{j = 0}^{m} \overset{⏞}{Θ_{(k, j)}} ρ_{j}^{τ} (t)

Truncating the sum and writing in modified form we get

ϕ_{n} {(t)}^{c} D^{σ} u (t) - \sum_{k = 0}^{m} c_{k} \sum_{j = 0}^{m} \overset{⏞}{Θ_{(k, j)}} ρ_{j}^{τ} (t) = \sum_{k = m + 1}^{\infty} c_{k} \sum_{j = 0}^{m} \overset{⏞}{Θ_{(k, j)}} ρ_{j}^{τ} (t) .

We can also write it in matrix form as

ϕ_{n} {(t)}^{c} D^{σ} u (t) - C_{M} B_{ϕ_{n}}^{σ} Λ_{M}^{τ} (t) = \sum_{k = m + 1}^{\infty} c_{k} \sum_{j = 0}^{m} \overset{⏞}{Θ_{(k, j)}} ρ_{j}^{τ} (t)

Using the fact

ρ_{j}^{τ} (t) \leq 1

for

t \in [0, τ]

, therefore, we can write

| ϕ_{n} {(t)}^{c} D^{σ} u (t) - C_{M} B_{ϕ_{n}}^{σ} Λ_{M}^{τ} (t) | \leq | \sum_{k = m + 1}^{\infty} c_{k} \sum_{j = 0}^{m} \overset{⏞}{Θ_{(k, j)}} | .

Furthermore, hence the proof is complete. □

The above matrices have been successfully applied to solve fractional-order differential equations (FDEs) under the effect of initial conditions. However these matrices do not have the ability to solve FDEs with integral types of boundary conditions. Therefore, the invention of new matrices that can easily handle integral types of boundary conditions are of basic importance. In the forthcoming discussion, we will design two new operational matrices having the ability to deal with integral type boundary conditions.

Lemma 5.

Let

s (t)

be a known function well defined on

[0, τ]

and

ϕ_{c}^{n} (t) = c t^{n}

be a polynomial then, for a function vector

Λ_{M}^{τ} (t)

as defined in (8), the following result holds

ϕ_{c}^{n} (t) \int_{0}^{τ} s (t) Λ_{M}^{τ} (t) d t = Q_{M \times M}^{c, n, s (t)} Λ_{M}^{τ} (t),

where the matrix

Q_{M \times M}^{c, n, s (t)}

is given as

Q_{M \times M}^{c, n, s (t)} = (\begin{matrix} Ω_{(0, 0)} & Ω_{(0, 1)} & \dots & Ω_{(0, m)} \\ Ω_{(1, 0)} & Ω_{(1, 1)} & \dots & Ω_{(1, m)} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ Ω_{(m, 0)} & Ω_{(m, 1)} & \dots & Ω_{(m, m)} \end{matrix}) .

where

Ω_{(i, j)} = \sum_{r = 0}^{j} \frac{τ c d_{i} τ^{r + n + 1} ℷ_{(j, r)}}{(2 i + 1) (r + n + 1)},

d_{i}

is the Legendre coefficients of the function

s (t)

and

ℷ_{(j, r)}

is as defined in (4).

Proof.

Let

s (t)

be approximated with Legendre polynomials, as

s (t) = \sum_{l = 0}^{m} d_{l} ϕ_{c}^{n} (t) .

We can write the

i^{t h}

term of

ϕ_{c}^{n} (t) \int_{0}^{τ} s (t) Λ_{M}^{τ} (t) d t

\begin{matrix} c t^{n} \int_{0}^{τ} \sum_{l = 0}^{m} d_{l} ρ_{l}^{τ} (t) ρ_{i}^{τ} (t) d t & = c t^{n} \sum_{l = 0}^{m} d_{l} \int_{0}^{τ} ρ_{l}^{τ} (t) ρ_{i}^{τ} (t) d t, \\ = \frac{τ c d_{i} t^{n}}{2 i + 1} . \end{matrix}

The polynomial

\frac{τ c d_{i} t^{n}}{2 i + 1}

can be expressed as a series of Legendre polynomials as

\frac{τ c d_{i} t^{n}}{2 i + 1} = \sum_{j = 0}^{m} Ω_{(i, j)} ρ_{j} (t) .

where the constant

Ω_{(i, j)}

is given by

Ω_{(i, j)} = \frac{τ c d_{i}}{2 i + 1} \int_{0}^{τ} ρ_{j}^{τ} (t) t^{n} d t .

Using the definition of

ρ_{j}^{τ} (t)

and after simplification, we can write

\begin{matrix} Ω_{(i, j)} & = \frac{τ c d_{i}}{2 i + 1} \sum_{r = 0}^{j} ℷ_{(j, r)} \int_{0}^{τ} t^{r + n} d t, \\ = \sum_{r = 0}^{j} \frac{τ c d_{i} τ^{r + n + 1} ℷ_{(j, r)}}{(2 i + 1) (r + n + 1)} . \end{matrix}

Simulating the result for

i = 0 : M

and

j = 0 : M

completes the proof of the Lemma. □

Lemma 6.

Let

ϕ_{c}^{n} (n) = c t^{n}

be a polynomial then for a function vector

Λ_{M}^{τ} (t)

as defined in (8); the following holds

ϕ_{c}^{n} (t) Λ_{M}^{τ} (τ) = R_{M \times M}^{c, n, τ} Λ_{M}^{τ} (t)

The matrix

R_{M \times M}^{c, n, τ}

is defined as

R_{M \times M}^{c, n, τ} = (\begin{matrix} ℑ_{(0, 0)} & ℑ_{(0, 1)} & \dots & ℑ_{(0, m)} \\ ℑ_{(1, 0)} & ℑ_{(1, 1)} & \dots & ℑ_{(1, m)} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ ℑ_{(m, 0)} & ℑ_{(m, 1)} & \dots & ℑ_{(m, m)} \end{matrix}) .

Furthermore, the entries are defined as

ℑ_{(i, j)} = \sum_{k = 0}^{i} \sum_{l = 0}^{j} \frac{c τ^{k + n + l + 1} ℷ_{(i, k)} ℷ_{(j, l)}}{n + l + 1} .

Proof.

The general term of the equality can be written as

ϕ_{c}^{n} (t) ρ_{i}^{τ} (τ) = \sum_{k = 0}^{i} ℷ_{(i, k)} τ^{k} c t^{n} .

It can be approximated with Legendre polynomials

\sum_{k = 0}^{i} ℷ_{(i, k)} τ^{k} c t^{n} = \sum_{j = 0}^{m} ℑ_{(i, j)} ρ_{j}^{τ} (t),

where

ℑ_{(i, j)} = \int_{0}^{τ} \sum_{k = 0}^{i} ℷ_{(i, k)} τ^{k} c t^{n} ρ_{j}^{τ} (t) d t .

Using the definition of Legendre polynomials we can write

ℑ_{(i, j)} = \sum_{k = 0}^{i} \sum_{l = 0}^{j} \frac{c τ^{k + n + l + 1} ℷ_{(i, k)} ℷ_{(j, l)}}{n + l + 1} .

which completes the proof of the lemma. □

4. Application of Operational Matrices

The operational matrix method for the solution of fractional differential equations is, in fact, a spectral method. The main aim of the spectral method is to convert a typical differential equation to system of easily solvable algebraic equations, which can be solved to obtain the solution in the series form of some orthogonal polynomials. The application of these methods to nonlinear differential equations results in a nonlinear system of algebraic equations, which can be solved using some iterative algorithms (the Newton–Raphson method is a frequently used method), see for example [46,47,48,49,50,51,52,53].

In this paper, we implement a different approach. We first use the Taylor series method to linearize the nonlinear part f and g to convert the nonlinear fractional differential equation into a recurrence relation of linear fractional differential equations with variable coefficients.

4.1. Linear FDEs with Variable Coefficients

We first consider the following coupled system of linear fractional differential equations with variable coefficients

\begin{matrix} ^{c} D^{σ_{1}} u (t) = c_{0} (t) u + c_{1} (t) u^{(γ_{1})} + c_{2} (t) v + c_{3} (t) v^{(γ_{1})} + h (t), \\ ^{c} D^{σ_{2}} v (t) = d_{0} (t) u + d_{1} (t) u^{(γ_{1})} + d_{2} (t) v + d_{3} (t) v^{(γ_{1})} + k (t), \\ u (0) = u_{0}, u (τ) = m_{1} \int_{0}^{τ} s (t) u (t) d t, \\ v (0) = v_{0}, v (τ) = m_{2} \int_{0}^{τ} r (t) v (t) d t, \end{matrix}

(10)

where

m_{1}

and

m_{2}

are some real constants.

0 < γ_{1}, γ_{2} \leq 1,

1 < σ_{1}, σ_{2} \leq 2,

t \in [0, τ]

and

c_{i} (t), d_{i} (t), s (t), r (t), h (t)

and

k (t)

are bounded, continuous, and well-defined functions on the domain

[0, τ] .

Assume the solution of (10) in terms of shifted Legendre polynomials, such that the following hold

^{c} D^{σ_{1}} u (t) = K_{M} Λ_{M}^{τ} (t),^{c} D^{σ_{2}} v (t) = L_{M} Λ_{M}^{τ} (t) .

(11)

Applying fractional integral of order

σ_{1}

on the first equation of (11) and making use of Lemma 2, we can write

u (t) = K_{M} H_{M \times M}^{τ, σ_{1}} Λ_{M}^{τ} (t) + e_{0} + e_{1} t .

(12)

By the application of initial conditions, we get

e_{0} = u_{0},

and to get the value of

e_{1}

, we use the integral-type boundary conditions, and after simplification it follows that

\begin{matrix} K_{M} & H_{M \times M}^{τ, σ_{1}} Λ_{M}^{τ} (τ) + u_{0} + e_{1} τ = \\ m_{1} K_{M} H_{M \times M}^{τ, σ_{1}} \int_{0}^{τ} s (t) Λ_{M}^{τ} (t) d t \\ + m_{1} \int_{0}^{τ} s (t) u_{0} d t + m_{1} e_{1} \int_{0}^{τ} s (t) t d t, \end{matrix}

\begin{matrix} e_{1} = \frac{1}{s_{1}} (K_{M} H_{M \times M}^{τ, σ_{1}} Λ_{M}^{τ} (τ) \\ - m_{1} K_{M} H_{M \times M}^{τ, σ_{1}} \int_{0}^{τ} s (t) Λ_{M}^{τ} (t) d t - s_{0}), \end{matrix}

where

s_{1} = (m_{1} \int_{0}^{τ} s (t) t d t - τ)

and

s_{0} = m_{1} \int_{0}^{τ} s (t) u_{0} d t - u_{0} .

Now using the values of

e_{0}

and

e_{1}

in (12), we can write

u (t)

\begin{matrix} u (t) = K_{M} H_{M \times M}^{τ, σ_{1}} Λ_{M}^{τ} (t) + u_{0} + \\ \frac{t}{s_{1}} (K_{M} H_{M \times M}^{τ, σ_{1}} Λ_{M}^{τ} (τ) - m_{1} K_{M} H_{M \times M}^{τ, σ_{1}} \int_{0}^{τ} s (t) Λ_{M}^{τ} (t) d t - s_{0}), \end{matrix}

(13)

In view of Lemma 5 and Lemma 6, we can write Equation (13) as

\begin{matrix} u (t) = & K_{M} H_{M \times M}^{τ, σ_{1}} Λ_{M}^{τ} (t) + K_{M} H_{M \times M}^{τ, σ_{1}} R_{M \times M}^{\frac{1}{s_{1}}, 1, τ} Λ_{M}^{τ} (t) \\ - K_{M} H_{M \times M}^{τ, σ_{1}} Q_{M \times M}^{\frac{m_{1}}{s_{1}}, 1, s (t)} Λ_{M}^{τ} (t) + {F_{1}}_{M} Λ_{M}^{τ} (t) . \end{matrix}

where

u_{0} - \frac{s_{0} t}{s_{1}} = F_{1} Λ_{M}^{τ} (t) .

For simplicity of notation, we can write the above equations as

u (t) = K_{M} {E^{(1)}}_{M \times M} Λ_{M}^{τ} (t) + {F_{1}}_{M} Λ_{M}^{τ} (t) .

(14)

where

{E^{(1)}}_{M \times M} = H_{M \times M}^{τ, σ_{1}} (I_{M \times M} + R_{M \times M}^{\frac{1}{s_{1}}, 1, τ} - Q_{M \times M}^{\frac{m_{1}}{s_{1}}, 1, s (t)})

Repeating the same procedure for

v (t)

, we can get analogous representation as

v (t) = L_{M} {E^{(2)}}_{M \times M} Λ_{M}^{τ} (t) + {F_{2}}_{M} Λ_{M}^{τ} (t) .

(15)

where

E_{M \times M}^{(2)} = H_{M \times M}^{τ, σ_{2}} (I_{M \times M} + R_{M \times M}^{\frac{1}{r_{1}}, 1, τ} - Q_{M \times M}^{\frac{m_{2}}{r_{1}}, 1, r (t)})

Now, in view of (15), (14), (11), and Lemma 4, the equivalent matrix form for system (10) is given as

\begin{matrix} K_{M} Λ_{M}^{τ} (t) = & K_{M} E_{M \times M}^{(1)} (B_{M \times M}^{c_{0}, 0} + B_{M \times M}^{c_{1}, γ_{1}}) Λ_{M}^{τ} (t) + \\ {F_{1}}_{M} (B_{M \times M}^{c_{0}, 0} + B_{M \times M}^{c_{1}, γ_{1}}) Λ_{M}^{τ} (t) + \\ L_{M} E_{M \times M}^{(2)} (B_{M \times M}^{c_{2}, 0} + B_{M \times M}^{c_{3}, γ_{2}}) Λ_{M}^{τ} (t) + \\ {F_{2}}_{M} (B_{M \times M}^{c_{2}, 0} + B_{M \times M}^{c_{3}, γ_{2}}) Λ_{M}^{τ} (t) + {Z_{1}}_{M} Λ_{M}^{τ} (t), \\ L_{M} Λ_{M}^{τ} (t) = & K_{M} E_{M \times M}^{(1)} (B_{M \times M}^{d_{0}, 0} + B_{M \times M}^{d_{1}, γ_{1}}) Λ_{M}^{τ} (t) + \\ {F_{1}}_{M} (B_{M \times M}^{d_{0}, 0} + B_{M \times M}^{d_{1}, γ_{1}}) Λ_{M}^{τ} (t) + \\ L_{M} E_{M \times M}^{(2)} (B_{M \times M}^{d_{2}, 0} + B_{M \times M}^{d_{3}, γ_{2}}) Λ_{M}^{τ} (t) + \\ {F_{2}}_{M} (B_{M \times M}^{d_{2}, 0} + B_{M \times M}^{d_{3}, γ_{2}}) Λ_{M}^{τ} (t) + {Z_{2}}_{M} Λ_{M}^{τ} (t) \end{matrix}

(16)

Canceling out the common term and after simplification of notation, we can write

\begin{matrix} K_{M} = & K_{M} E_{M \times M}^{(1)} (B_{M \times M}^{c_{0}, 0} + B_{M \times M}^{c_{1}, γ_{1}}) + \\ L_{M} E_{M \times M}^{(2)} (B_{M \times M}^{c_{2}, 0} + B_{M \times M}^{c_{3}, γ_{2}}) + {Y_{1}}_{M} \\ L_{M} = & K_{M} E_{M \times M}^{(1)} (B_{M \times M}^{d_{0}, 0} + B_{M \times M}^{d_{1}, γ_{1}}) + \\ L_{M} E_{M \times M}^{(2)} (B_{M \times M}^{d_{2}, 0} + B_{M \times M}^{d_{3}, γ_{2}}) + {Y_{2}}_{M}, \end{matrix}

(17)

where

{Y_{1}}_{M} = {F_{2}}_{M} (B_{M \times M}^{c_{2}, 0} + B_{M \times M}^{c_{3}, γ_{2}}) +

{F_{1}}_{M} (B_{M \times M}^{c_{0}, 0} + B_{M \times M}^{c_{1}, γ_{1}}) + {Z_{1}}_{M},

and

{Y_{2}}_{M} = {F_{2}}_{M} (B_{M \times M}^{d_{2}, 0} + B_{M \times M}^{d_{3}, γ_{2}}) +

{F_{1}}_{M} (B_{M \times M}^{d_{0}, 0} + B_{M \times M}^{d_{1}, γ_{1}}) + {Z_{2}}_{M} .

The above equations can also be written as

\begin{matrix} [K_{M} L_{M}] = [K_{M} L_{M}] \\ [\begin{matrix} E^{(1)} (B^{c_{0}, 0} + B^{c_{1}, γ_{1}}) & E^{(1)} (B^{d_{0}, 0} + B^{d_{1}, γ_{1}}) \\ E^{(2)} (B^{c_{2}, 0} + B^{c_{3}, γ_{2}}) & E^{(2)} (B^{d_{2}, 0} + B^{d_{3}, γ_{2}}) \end{matrix}] \\ + [\begin{matrix} {Y_{1}}_{M} & {Y_{2}}_{M} \end{matrix}] \end{matrix}

(18)

We see that (18) is a linear system of matrix equations. By solving (17), we will get the required coefficients vector

K_{M}

and

L_{M}

, which can be used in (14) and (15) to get approximation to the solution of the main problem.

4.2. Nonlinear FDEs

Nonlinear FDEs cannot be directly solved using the OP method; however, combining it with the quasilinearization method makes it easy to recursively solve nonlinear FDEs. The procedure of this technique is as given as follows.

Approximate the initial solution, the solution of the linear part, by the method presented in previous section and name it $u_{0} (t)$ and $v_{0} (t) .$
Linearize the nonlinear part at $u_{0} (t)$ and $v_{0} (t)$ . This will convert the system of nonlinear FDEs into a system of linear FDEs that is easily solvable with the method devolved. Solve it and name the solution as $u_{1} (t)$ and $v_{1} (t)$ .
Repeat step 1.

Consider the following nonlinear FDEs.

\begin{matrix} ^{c} D^{σ_{1}} u (t) = f (u, v, u^{(γ_{1})}, v^{(γ_{2})}), \\ ^{c} D^{σ_{2}} v (t) = g (u, v, u^{(γ_{1})}, v^{(γ_{2})}), \\ u (0) = a_{0}, u (τ) = \int_{0}^{ω_{1}} s (t) u (t) d t, 0 < ω_{1} \leq τ, \\ v (0) = b_{0}, v (τ) = \int_{0}^{ω_{2}} r (t) v (t) d t, 0 < ω_{2} \leq τ . \end{matrix}

(19)

separating the linear and nonlinear parts, we get

\begin{matrix} ^{c} D^{σ_{1}} u (t) = L_{1} (u, v, u^{(γ_{1})}, v^{(γ_{2})}) + N_{1} (u, v, u^{(γ_{1})}, v^{(γ_{2})}), \\ ^{c} D^{σ_{2}} v (t) = L_{2} (u, v, u^{(γ_{1})}, v^{(γ_{2})}) + N_{2} (u, v, u^{(γ_{1})}, v^{(γ_{2})}), \end{matrix}

(20)

First solve the linear part:

\begin{matrix} ^{c} D^{σ_{1}} u (t) = L_{1} (u, v, u^{(γ_{1})}, v^{(γ_{2})}), \\ ^{c} D^{σ_{2}} v (t) = L_{2} (u, v, u^{(γ_{1})}, v^{(γ_{2})}), \end{matrix}

(21)

Its solution is named

u_{0} (t)

and

v_{0} (t)

. The next step is to linearize the nonlinear part.

\begin{matrix} ^{c} D^{σ_{1}} u_{1} (t) & = L_{1} (u_{1}, v_{1}, u_{1}^{(γ_{1})}, v_{1}^{(γ_{2})}) + N_{1} (u_{0}, v_{0}, u_{0}^{(γ_{1})}, v_{0}^{(γ_{2})}) \\ + (u_{1} - u_{0}) \frac{\partial N_{1}}{\partial u_{0}} + (v_{1} - v_{0}) \frac{\partial N_{1}}{\partial v_{0}} + \\ (u_{1}^{(γ_{1})} - u_{0}^{(γ_{1})}) \frac{\partial N_{1}}{\partial u_{0}^{(γ_{1})}} + (v_{1}^{(γ_{2})} - v_{0}^{(γ_{2})}) \frac{\partial N_{1}}{\partial v_{0}^{(γ_{2})}}, \\ ^{c} D^{σ_{2}} v_{1} (t) & = L_{2} (u_{1}, v_{1}, u_{1}^{(γ_{1})}, v_{1}^{(γ_{2})}) + N_{2} (u_{0}, v_{0}, u_{0}^{(γ_{1})}, v_{0}^{(γ_{2})}) \\ + (u_{1} - u_{0}) \frac{\partial N_{2}}{\partial u_{0}} + (v_{1} - v_{0}) \frac{\partial N_{2}}{\partial v_{0}} + \\ (u_{1}^{(γ_{1})} - u_{0}^{(γ_{1})}) \frac{\partial N_{2}}{\partial u_{0}^{(γ_{1})}} + (v_{1}^{(γ_{2})} - v_{0}^{(γ_{2})}) \frac{\partial N_{2}}{\partial v_{0}^{(γ_{2})}} . \end{matrix}

(22)

We get a system of FDEs with variable coefficients. The whole process can be expressed as

\begin{matrix} ^{c} D^{σ_{1}} u_{r + 1} (t) & = L_{1} (u_{r + 1}, v_{r + 1}, u_{r + 1}^{(γ_{1})}, v_{r + 1}^{(γ_{2})}) + \\ N_{1} (u_{r}, v_{r}, u_{r}^{(γ_{1})}, v_{r}^{(γ_{2})}) \\ + (u_{r + 1} - u_{r}) \frac{\partial N_{1}}{\partial u_{r}} + (v_{r + 1} - v_{r}) \frac{\partial N_{1}}{\partial v_{r}} + \\ (u_{r + 1}^{(γ_{1})} - u_{r}^{(γ_{1})}) \frac{\partial N_{1}}{\partial u_{r}^{(γ_{1})}} + (v_{r + 1}^{(γ_{2})} - v_{r}^{(γ_{2})}) \frac{\partial N_{1}}{\partial v_{r}^{(γ_{2})}} \\ ^{c} D^{σ_{2}} v_{r + 1} (t) & = L_{2} (u_{r + 1}, v_{r + 1}, u_{r + 1}^{(γ_{1})}, v_{r + 1}^{(γ_{2})}) + \\ N_{2} (u_{r}, v_{r}, u_{r}^{(γ_{1})}, v_{r}^{(γ_{2})}) + (u_{r + 1} - u_{r}) \frac{\partial N_{2}}{\partial u_{r}} \\ + (v_{r + 1} - v_{r}) \frac{\partial N_{2}}{\partial v_{r}} + (u_{r + 1}^{(γ_{1})} - u_{r}^{(γ_{1})}) \frac{\partial N_{2}}{\partial u_{r}^{(γ_{1})}} \\ + (v_{r + 1}^{(γ_{2})} - v_{r}^{(γ_{2})}) \frac{\partial N_{2}}{\partial v_{r}^{(γ_{2})}} \end{matrix}

(23)

The boundary conditions can be written as

\begin{matrix} u_{r + 1} (0) = a_{0}, u_{r + 1} (τ) = \int_{0}^{ω_{1}} s (t) u_{r + 1} (t) d t, 0 < ω_{1} \leq τ, \\ v_{r + 1} (0) = b_{0}, v_{r + 1} (τ) = \int_{0}^{ω_{2}} r (t) v_{r + 1} (t) d t, 0 < ω_{2} \leq τ . \end{matrix}

(24)

It can be easily noted that (23) is fractional differential equation with variable coefficients.

5. Error Bound of the Approximate Solution and Convergence

In this section, we calculate a upper bound for error of approximation of solution with the proposed method.

5.1. Error Bound for Single Differential Equation

Consider the following fractional differential equation.

^{c} D^{σ} u (t) = c_{0} (t) u (t) + c_{1} (t) u {(t)}^{(γ)} + h (t),

(25)

subject to the following initial and boundary conditions

u (0) = u_{0}, u (τ) = m_{1} \int_{0}^{τ} s (t) u (t) d t .

Our aim is to derive an upper bound for the proposed method. We have to calculate

R_{M}

defined as

R_{M} = |^{c} D^{σ} u (t) - K_{M} G_{M \times M}^{τ, σ} Λ_{M}^{τ} (t) | .

(26)

The solution of the above problem can be written in terms of shifted Legendre series such that

^{c} D^{σ} u (t) = \sum_{k = 0}^{\infty} u_{k} ρ_{k}^{τ} (t) = K_{M} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} u_{k} ρ_{k}^{τ} (t) .

(27)

Applying a fractional integral of order

σ

, using an operational matrix of integration and using Corollary (1), we can write

u (t) - c_{0} - c_{1} t = K_{M} H_{M \times M}^{τ, σ} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} c_{k} \sum_{i = 0}^{m} Θ_{i, k, τ} ρ_{i}^{τ} (t) .

(28)

Which can be simplified as

u (t) = K_{M} H_{M \times M}^{τ, σ} Λ_{M}^{τ} (t) + c_{0} + c_{1} t + \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} Θ_{i, k, τ} ρ_{i}^{τ} (t),

(29)

We know from (14) and the integral type boundary conditions that we can conclude,

u (t) = K_{M} E_{M \times M} Λ_{M}^{τ} (t) + {F_{1}}_{M} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} Θ_{i, k, τ} ρ_{i}^{τ} (t),

Assume that

\hat{X} = K_{M} E_{M \times M} + F_{1} .

Using Lemma (3) and Corollary (2), we can write

c_{l} (t) u {(t)}^{(γ)} = {\hat{X}}_{M}^{T} B_{c_{l}}^{γ} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(l)}} ρ_{i^{'}}^{τ} (t) .

(30)

Approximating

h (t) = F_{2} Λ_{M}^{τ} (t) + \sum_{k^{'} = m + 1}^{\infty} f_{k^{'}} ρ_{k^{'}}^{τ} (t)

and using (30) in (25) we get

\begin{matrix} K_{M} Λ_{M}^{τ} (t) - {\hat{X}}_{M}^{T} B_{c_{0}}^{0} Λ_{M}^{τ} (t) - {\hat{X}}_{M}^{T} B_{c_{1}}^{γ} Λ_{M}^{τ} (t) - F_{2} Λ_{M}^{τ} (t) = R_{M} (t) . \end{matrix}

where

R_{M} (t)

is defined by relation

\begin{matrix} R_{M} (t) = & \sum_{k = m + 1}^{\infty} u_{k} ρ_{k}^{τ} (t) + \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} ρ_{i^{'}}^{τ} (t) + \\ \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}} ρ_{i^{'}}^{τ} (t) + \sum_{k^{'} = m + 1}^{\infty} f_{k^{'}} ρ_{k^{'}}^{τ} (t), \end{matrix}

\begin{matrix} R_{M} (t) = & \sum_{k = m + 1}^{\infty} u_{k} [ρ_{k}^{τ} (t) + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} ρ_{i^{'}}^{τ} (t) + \\ \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}} ρ_{i^{'}}^{τ} (t)] + \sum_{k^{'} = m + 1}^{\infty} f_{k^{'}} ρ_{k^{'}}^{τ} (t) . \end{matrix}

Using the bounded property of Legendre polynomail, it follows that

\begin{matrix} | R_{M} (t) | \leq & \sum_{k = m + 1}^{\infty} | u_{k} | | [1 + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}}] | + \sum_{k^{'} = m + 1}^{\infty} | f_{k^{'}} | . \end{matrix}

(31)

In view of Theorem (1), it is evident that if the function

u (t)

and

f (t)

are sufficiently smooth functions, then the sequence that defines their coefficient is convergent to zero. Hence, we conclude that as

m \to \infty

the coefficients

u_{m} \to 0

and

f_{m} \to 0

. Hence it can be easily observed that the error

| R_{M} (t) | \to 0 .

Equation (31) also establishes an upper bound of the error between the exact and approximate solution.

5.2. Error Bound for Coupled System of Fractional Differential Equations

Consider the following system of FDEs.

\begin{matrix} ^{c} D^{σ} u (t) = c_{0} (t) u (t) + c_{1} (t) v (t) + c_{2} (t) u {(t)}^{(γ)} + c_{3} (t) v {(t)}^{(γ)} + h (t), \\ ^{c} D^{σ} v (t) = d_{0} (t) u (t) + d_{1} (t) v (t) + d_{2} (t) u {(t)}^{(γ)} + d_{3} (t) v {(t)}^{(γ)} + g (t), \end{matrix}

(32)

subject to the following initial and boundary conditions

\begin{matrix} u (0) = u_{0}, u (τ) = m_{1} \int_{0}^{τ} s (t) u (t) d t . \\ v (0) = v_{0}, v (τ) = m_{2} \int_{0}^{τ} r (t) v (t) d t . \end{matrix}

Our aim is to derive an upper bound for the proposed method. We have to calculate

R_{M}^{u}

and

R_{M}^{v}

, defined as

\begin{matrix} R_{M}^{u} = |^{c} D^{σ} u (t) - K_{M} G_{M \times M}^{τ, σ} Λ_{M}^{τ} (t) |, \\ R_{M}^{v} = |^{c} D^{σ} v (t) - L_{M} G_{M \times M}^{τ, σ} Λ_{M}^{τ} (t) | . \end{matrix}

(33)

We know from (14) that

\begin{matrix} u (t) = K_{M} E_{M \times M} Λ_{M}^{τ} (t) + {F_{1}}_{M} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} Θ_{i, k, τ} ρ_{i}^{τ} (t), \\ v (t) = L_{M} E_{M \times M} Λ_{M}^{τ} (t) + {F_{2}}_{M} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} v_{k} \sum_{i = 0}^{m} Θ_{i, k, τ} ρ_{i}^{τ} (t), \end{matrix}

Assume,

\hat{X} = K_{M} E_{M \times M} + F_{1}

and

\hat{Y} = L_{M} E_{M \times M} + F_{2} .

Using Lemma (3) and Corollary (2), we can write

\begin{matrix} c_{l} (t) u {(t)}^{(γ)} = {\hat{X}}_{M}^{T} B_{c_{l}}^{γ} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(l)}} ρ_{i^{'}}^{τ} (t), \\ d_{l} (t) v {(t)}^{(γ)} = {\hat{Y}}_{M}^{T} B_{d_{l}}^{γ} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} v_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(l)}} ρ_{i^{'}}^{τ} (t) . \end{matrix}

(34)

Approximating

h (t) = D_{1} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty} d_{k} ρ_{k}^{τ} (t)

and

g (t) = D_{2} Λ_{M}^{τ} (t) + \sum_{k = m + 1}^{\infty}

d_{k}^{'} ρ_{k}^{τ} (t)

and using (30) in (25) we get

\begin{matrix} K_{M} Λ_{M}^{τ} (t) - {\hat{X}}_{M}^{T} B_{c_{0}}^{0} Λ_{M}^{τ} (t) - {\hat{Y}}_{M}^{T} B_{c_{1}}^{0} Λ_{M}^{τ} (t) - {\hat{X}}_{M}^{T} B_{c_{2}}^{γ} Λ_{M}^{τ} (t) - {\hat{Y}}_{M}^{T} B_{c_{3}}^{γ} Λ_{M}^{τ} (t) - D_{1} Λ_{M}^{τ} (t) = R_{M}^{u} (t), \\ L_{M} Λ_{M}^{τ} (t) - {\hat{X}}_{M}^{T} B_{d_{0}}^{0} Λ_{M}^{τ} (t) - {\hat{Y}}_{M}^{T} B_{d_{1}}^{0} Λ_{M}^{τ} (t) - {\hat{X}}_{M}^{T} B_{d_{2}}^{γ} Λ_{M}^{τ} (t) - {\hat{Y}}_{M}^{T} B_{d_{3}}^{γ} Λ_{M}^{τ} (t) - D_{2} Λ_{M}^{τ} (t) = R_{M}^{v} (t) . \end{matrix}

where

R_{M}^{u} (t)

and

R_{M}^{v} (t)

is defined by the relation

\begin{matrix} R_{M}^{u} (t) = & \sum_{k = m + 1}^{\infty} u_{k} ρ_{k}^{τ} (t) + \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} ρ_{i^{'}}^{τ} (t) + \sum_{k = m + 1}^{\infty} v_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} ρ_{i^{'}}^{τ} (t) + \\ \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}} ρ_{i^{'}}^{τ} (t) + \sum_{k = m + 1}^{\infty} v_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}} ρ_{i^{'}}^{τ} (t) + \sum_{k = m + 1}^{\infty} d_{k} ρ_{k}^{τ} (t), \\ R_{M}^{v} (t) = & \sum_{k = m + 1}^{\infty} v_{k} ρ_{k}^{τ} (t) + \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} ρ_{i^{'}}^{τ} (t) + \sum_{k = m + 1}^{\infty} v_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} ρ_{i^{'}}^{τ} (t) + \\ \sum_{k = m + 1}^{\infty} u_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}} ρ_{i^{'}}^{τ} (t) + \sum_{k = m + 1}^{\infty} v_{k} \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}} ρ_{i^{'}}^{τ} (t) + \sum_{k = m + 1}^{\infty} d_{k}^{'} ρ_{k}^{τ} (t), \end{matrix}

Using the bounded property of the Legendre polynomial, it follows that

\begin{matrix} | R_{M}^{u} (t) | \leq & \sum_{k = m + 1}^{\infty} | u_{k} | | [1 + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}}] | + \\ \sum_{k = m + 1}^{\infty} | v_{k} | | [\sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}}] | + \sum_{k = m + 1}^{\infty} | d_{k} |, \\ | R_{M}^{v} (t) | \leq & \sum_{k = m + 1}^{\infty} | v_{k} | | [1 + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}}] | + \\ \sum_{k = m + 1}^{\infty} | u_{k} | | [\sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(0)}} + \sum_{i = 0}^{m} \sum_{i^{'} = 0}^{m} Θ_{i, k, τ} \overset{⏞}{Θ_{k, i^{'}, τ}^{(1)}}] | + \sum_{k = m + 1}^{\infty} | d_{k}^{'} | . \end{matrix}

(35)

The above equation establishes an error bound for the solution

u (t)

and

v (t) .

It also ensures the convergence of the proposed method for the solution of coupled system of FDEs.

6. Test Problems

We solve one single equation, three systems of linear FDEs with variable equations, and two systems of nonlinear problems, and analyze the convergence of the approximate solution by measuring the following error norms.

| | E_{u} {| |}_{2} = \frac{1}{τ} \int_{0}^{τ} U_{M} (t) - u (t) d t

and

| | E_{u} {| |}_{m a x} = m a x_{x \in [0, τ]} | U_{M} (t) - u (t) | .

We check the accuracy of the boundary condition by measuring the following error norms.

| | E_{u} {| |}_{b} = | U_{M} (τ) - m_{1} \int_{0}^{τ} U_{M} (t) | .

In the above bounds, the quantity

U_{M} (t)

represents the

m -

term approximation to the solution

u (t)

Test Problem 1.

\begin{matrix} ^{c} D^{1.2} u (t) = (t^{2} + s i n (t)) u (t) & + (2 t - t^{3}) u^{0.7} (t) + f (t), \\ u (0) = 4, & u (1) = 1.1216 \int_{0}^{1} c o s (t) u (t) d t, \end{matrix}

where the exact solution

u (t) = t^{3} + t^{2} + t + 4

, and the source term

f (t) = \frac{38683084397149375 t^{\frac{4}{5}} (50 t^{2} - 35 t + 21)}{378302368699121664} - (s i n (t) + t^{2}) (t^{4} - t^{3} + t^{2} + 4) - \frac{150543064388819875 t^{\frac{13}{10}} (2 t - t^{3}) (400 t^{2} - 330 t + 253)}{22218508761632342016} .

Test Problem 2.

\begin{matrix} ^{c} D^{σ} u (t) = (t + 1) u (t) & + (1 - t) u^{σ_{1}} (t) + (2 t) v (t) \\ + (t^{2}) v^{'} (t) + f (t), \\ ^{c} D^{σ} v (t) = (t^{2} + 1) u (t) & + (1 - t^{2}) u^{σ_{1}} (t) + (3 t) v (t) \\ + (t^{3}) v^{'} (t) + g (t), \\ u (0) = 1, & u (1) = 2.1270 \int_{0}^{1} s i n (t) u (t) d t, \\ v (0) = - 1, & v (1) = - 1.8925 \int_{0}^{1} c o s (t) v (t) d t . \end{matrix}

where the exact solution

u (t) = t^{2} + t^{3} + e^{t}

and

v (t) = t^{2} - t^{3} - e^{t}

, and the source term

f (t) = 6 t + e^{t} + 2 t (e^{t} - t^{2} + t^{3}) - (t + 1) (e^{t} + t^{2} + t^{3}) + (t - 1) (2 t + e^{t} + 3 t^{2}) + t^{2} (e^{t} - 2 t + 3 t^{2}) + 2

and

g (t) = 3 t (e^{t} - t^{2} + t^{3}) - e^{t} - 6 t + t^{3} (e^{t} - 2 t + 3 t^{2}) - (t^{2} + 1) (e^{t} + t^{2} + t^{3}) + (t^{2} - 1) (2 t + e^{t} + 3 t^{2}) + 2

Test Problem 3.

\begin{matrix} ^{c} D^{σ} u (t) = & s i n (t) u (t) + c o s (t) u^{σ_{1}} (t) + \\ (s i n (t) + c o s (t)) v (t) + s i n (2 t) v^{'} (t) + f (t), \\ ^{c} D^{σ} v (t) = & (c o s (2 t))) u (t) + (t s i n (t)) u^{σ_{1}} (t) \\ + (t c o s (t)) v (t) + (t^{2} s i n (t)) v^{'} (t) + g (t), \\ u (0) = 0, & u (1) = 0.8 \int_{0}^{1} e^{2 t} s i n (t) u (t) d t, \\ v (0) = 1, & v (1) = 0.71 \int_{0}^{1} e^{- 2 t} c o s (t) v (t) d t . \end{matrix}

where the exact solution

u (t) = e^{t} s i n (t)

and

v (t) = e^{- t} c o s (t)

, and the source term

f (t) = 2 e^{t} c o s (t) - e^{t} {s i n (t)}^{2} - c o s (t) (e^{t} c o s (t) + e^{t} s i n (t)) + s i n (2 t) (e^{- t} c o s (t) + e^{- t} s i n (t)) - e^{- t} c o s (t) (c o s (t) + s i n (t))

and

g (t) = 2 e^{- t} s i n (t) - t e^{- t} {c o s (t)}^{2} - t s i n (t) (e^{t} c o s (t) + e^{t} s i n (t)) - c o s (2 t) e^{t} s i n (t) + t^{2} s i n (t) (e^{- t} c o s (t) + e^{- t} s i n (t)) .

Test Problem 4.

\begin{matrix} ^{c} D^{σ} u (t) = e^{(t)} u (t) & + e^{(- t)} u^{σ_{1}} (t) + \\ (e^{(t)} + e^{(- t)}) v (t) + e^{(2 t)} v^{'} (t) + f (t), \\ ^{c} D^{σ} v (t) = (e^{(- 2 t)}) & u (t) + (t e^{(t)}) u^{σ_{1}} (t) \\ + (t e^{(t)}) v (t) + (t^{2} e^{(t)}) v^{'} (t) + g (t), \\ u (0) = 1, & u (1) = 0.5 \int_{0}^{1} (1 - t) e^{2 t} u (t) d t, \\ v (0) = 1, & v (1) = 10 \int_{0}^{1} (1 - t^{2}) e^{- 2 t} v (t) d t . \end{matrix}

where the exact solution

u (t) = t^{4} {(2 - t)}^{3}

and

v (t) = (t^{3}) (3 - t)

, and the source term

f (t) = e^{- t} (4 t^{3} {(t - 2)}^{3} + 3 t^{4} {(t - 2)}^{2}) - 12 t^{2} {(t - 2)}^{3} - 24 t^{3} {(t - 2)}^{2} - e^{2 t} (t^{3} (2 t - 6) + 3 t^{2} {(t - 3)}^{2}) - 3 t^{4} (2 t - 4) + t^{4} e^{t} {(t - 2)}^{3} - t^{3} (e^{- t} + e^{t}) {(t - 3)}^{2}

and

g (t) = 6 t {(t - 3)}^{2} + 6 t^{2} (2 t - 6) + 2 t^{3} - t^{4} e^{t} {(t - 3)}^{2} - t^{2} s i n (t) (t^{3} (2 t - 6) + 3 t^{2} {(t - 3)}^{2}) + t^{4} e^{- 2 t} {(t - 2)}^{3} + t e^{t} (4 t^{3} {(t - 2)}^{3} + 3 t^{4} {(t - 2)}^{2}) .

Test Problem 5.

\begin{matrix} ^{c} D^{σ} u (t) = u (t) + v (t) + u^{2} (t) + v^{'} (t) u^{'} (t) + f (t), \\ ^{c} D^{σ} v (t) = u^{'} (t) + v^{'} (t) + v {(t)}^{2} + u (t) u^{'} (t) + g (t), \\ u (0) = 0, u (1) = 5.4069 \int_{0}^{1} s i n (t) u (t) d t, \\ v (0) = 1, v (1) = 0 . \end{matrix}

where the exact solution

u (t) = (t^{2} + 1) t^{2}

and

v (t) = t^{3} (1 - t)

, and the source term

f (t) = t^{3} (t - 1) + (3 t^{2} (t - 1) + t^{3}) (2 t (t^{2} + 1) + 2 t^{3}) - t^{2} (t^{2} + 1) - t^{4} {(t^{2} + 1)}^{2} + 12 t^{2} + 2

and

g (t) = 3 t^{2} (t - 1) - 2 t (t^{2} + 1) - t^{6} {(t - 1)}^{2} - 6 t (t - 1) - 6 t^{2} - t^{3} - t^{2} (t^{2} + 1) (2 t (t^{2} + 1) + 2 t^{3}) .

Test Problem 6.

\begin{matrix} ^{c} D^{σ} u (t) = & (t + 1) u (t) + (1 - t) u^{'} (t) + 2 t v (t) \\ + t^{2} v^{'} (t) + u {(t)}^{2} + v^{'} (t) u^{'} (t) - v^{3} (t) + f (t), \\ ^{c} D^{σ} v (t) = & (t^{2} + 1) u (t) + (1 - t^{2}) u^{'} (t) + 3 t v (t) \\ + t^{3} v^{'} (t) + v {(t)}^{2} + u (t) u^{'} (t) - u^{4} (t) + g (t), \\ u (0) = 0, & u (1) = 0.6173 \int_{0}^{1} (2 t + 1) u (t) d t, \\ v (0) = 1, & v (1) = 0.7311 \int_{0}^{1} (3 t + 1) v (t) d t . \end{matrix}

where the exact solution

u (t) = s i n (t)

and

v (t) = e^{t}

, and the source term

f (t) = e^{3 t} - s i n (t) + c o s (t) (t - 1) - t^{2} e^{t} - e^{t} c o s (t) - s i n (t) (t + 1) - {s i n (t)}^{2} - 2 t e^{t}

and

g (t) = e^{t} - e^{2 t} - t^{3} e^{t} - c o s (t) s i n (t) + {sin (t)}^{4} + c o s (t) (t^{2} - 1) - s i n (t) (t^{2} + 1) - 3 t e^{t} .

7. Results and Discussion

The first test problem is solved using the proposed method. The problem is linear and is relatively easy to solve. We compare the approximate solution with the exact solution of the problem. We observe that by increasing the scale level of approximation, the approximate solution draws closer to the exact solution as expected; see for example Figure 1a. At

M = 8

the approximate solution (black line) coincides with the exact solution (red dots). In the second part of the same figure, we plot the absolute difference of the exact and approximate solution considering different scales. It is observed that at

M = 11

, the absolute error is less than

10^{- 9}

. This means accuracy up to the ninth decimal place is achieved. We also calculate all the three error norms at different scales. The results are displayed in Table 1. One can easily note that at

M = 5

, the value of is

| | E_{u} {| |}_{2} = 0.1314 \times 10^{- 1}

| | E_{u} {| |}_{m a x} = 0.6431 \times 10^{- 1}

and

| | E_{u} {| |}_{b} = 8.0032 \times 10^{- 17}

. While increasing the scale levels, these values start to decrease with great speed. At

M = 13

, these values become

4.5147 \times 10^{- 10},

2.4318 \times 10^{- 9}

and

2.4362 \times 10^{- 17}

, respectively.

The analysis of the second problem is given in Table 2. We fix the orders of equation to

σ = 1.8

and

γ = 0.8

and solve the problem using scale level ranges from

M = 5

M = 13 .

We observe that the error norm decreases rapidly with the increase in scale level. For example, the value

| | E_{u} {| |}_{2}

M = 5

1.1044 \times 10^{- 1}

and at

M = 13

, this value drops to

4.5177 \times 10^{- 10},

which is very high accuracy. At the same scale, the value of

| | E_{v} {| |}_{2}

1.3237 \times 10^{- 10}

. The value of norms

| | E_{u} {| |}_{m a x}

and

| | E_{v} {| |}_{m a x}

are

1.4501 \times 10^{- 9}

and

2.2449 \times 10^{- 10}

, respectively. Similarly the approximate solution satisfy the boundary condition with high accuracy. The errors in the boundary condition are

9.7194 \times 10^{- 17}

and

6.2399 \times 10^{- 16}

. The error in the boundary condition is observed to be constant, that is, we have the same accuracy at all scale level. The accuracy of the proposed method for all possible values of

σ_{1}

and

σ_{2}

are analyzed by calculating the error norm

| | E_{v} {| |}_{2}

. The error norm for the solution

u (t)

and

v (t)

are displayed in Figure 2. We observe that the method produces excellent approximation to the solution at almost all values of parameters.

The same analysis is performed for Test Problem 3 and Test Problem 4. The results are shown in Table 3 and Table 4. The same conclusion is reached for these two examples. The error bounds are shown in Figure 3 and Figure 4. We observe that the method yields an almost accurate solution for all values of these parameters. In Figure 2, we can easily see that the error norm is less than

10^{- 3}

. Note that for these problems we have set

M = 5

. It is always possible to get a more accurate solution by selecting the highest choices of scale level.

The nonlinear Test Problem 5 is solved with the proposed iterative scheme. We use three different choices of the parameters

σ

and

γ

and calculate the the error norms

| | E_{u} {| |}_{2}

at different iterations. We use seven iterations for the approximation of solution. The results are displayed in Figure 5 and Figure 6. We observe that the error decreases with respect to the number of iterations and is highly convergent. Note that for this example, we fix the scale level

M = 10

. It is clear form the figure that the proposed method is highly efficient, especially in the solution of nonlinear equations. The same phenomenon is observed for Test Problem 6. The results are displayed in Figure 7.

8. Conclusions and Future Work

This article presents a new algorithm for the solution of fractional order differential equations. The Newton Raphson method is combined with the operational matrix method for the solution of these problems. The convergence of the proposed method is checked analytically and is confirmed by solving several test problems. It is found that the approximate solution is highly accurate, and one can get high accuracy by using high scale levels. The mathematical proof of convergence and error analysis is our future plan of research.

Author Contributions

Conceptualization, H.K., M.K. and P.A.; investigation, H.K., M.K. and P.A.; writing—original draft, H.K., M.K., I.H. and P.A.; writing—review and editing, H.K., M.K., I.H. and P.A.; funding acquisition, I.H. All authors have read and agreed to the published version of the manuscript.

Funding

This study is supported by Universiti Kebangsaan Malaysia (Grant No: DIP-2021-018).

Data Availability Statement

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

Aytac, A.; Ozkol, I. Solution of fractional differential equations by using differential transform method. Chaos Solitons Fractals 2007, 34, 1473–1481. [Google Scholar]
Muhammet, K.; Bayram, M. Approximate analytical solution for the fractional modified KdV by differential transform method. Commun. Nonlinear Sci. Numer. Simul. 2010, 15, 1777–1782. [Google Scholar]
Che, H.C.H.; Kilicman, A. On the solution of fractional order nonlinear boundary value problems by using differential transformation method. Eur. J. Pure Appl. Math. 2011, 4, 174–185. [Google Scholar]
Shaher, M.; Odibat, Z. Homotopy perturbation method for nonlinear partial differential equations of fractional order. Phys. Lett. A 2007, 365, 345–350. [Google Scholar]
Zaid, O.; Momani, S. Modified homotopy perturbation method: Application to quadratic Riccati differential equation of fractional order. Chaos Solitons Fractals 2008, 36, 167–174. [Google Scholar]
Mehdi, D.; Manafian, J.; Saadatmandi, A. Solving nonlinear fractional partial differential equations using the homotopy analysis method. Numer. Methods Partial Differ. Equ. Int. J. 2010, 26, 448–479. [Google Scholar]
Saha, R.S.; Bera, R.K. An approximate solution of a nonlinear fractional differential equation by Adomian decomposition method. Appl. Math. Comput. 2005, 167, 561–571. [Google Scholar] [CrossRef]
Varsha, D.-G.; Jafari, H. Solving a multi-order fractional differential equation using Adomian decomposition. Appl. Math. Comput. 2007, 189, 541–548. [Google Scholar]
Shaher, M.; Odibat, Z. Numerical approach to differential equations of fractional order. J. Comput. Appl. Math. 2007, 207, 96–110. [Google Scholar]
Bhrawy, A.H.; Alofi, A.S.; Ezz-Eldien, S.Z. A quadrature tau method for variable coefficients fractional differential equations. Appl. Math. Lett. 2011, 24, 2146–2152. [Google Scholar] [CrossRef] [Green Version]
Bhrawy, A.H.; Alofi, A.S. A Jacobi–Gauss collocation method for solving nonlinear Lane-Emden type equations. Commun. Non. Sci. Numer. Simul. 2012, 17, 62–70. [Google Scholar] [CrossRef]
Bhrawy, A.H. A shifted Legendre spectral method for fractional-order multi-point boundary value problems. Adv. Differ. Equ. 2012, 8, 2012. [Google Scholar] [CrossRef] [Green Version]
Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. Spectral Methods in Fluid Dynamics; Springer: New York, NY, USA, 1988. [Google Scholar]
Canuto, C.; Hussaini, M.Y.; Quarteroni, A.; Zang, T.A. Spectral Methods: Fundamentals in Single Domains; Springer: New York, NY, USA, 2006. [Google Scholar]
Doha, E.H.; Bhrawy, A.H.; Hafez, R.M. A Jacobi-Jacobi dual-Petrov-Galerkin method for third- and fifth-order differential equations. Math. Comput. Model. 2011, 53, 1820–1832. [Google Scholar] [CrossRef]
Doha, E.H.; Bhrawy, A.H.; Saker, M.A. Integrals of Bernstein polynomials: An application for the solution of high even-order differential equations. Appl. Math. Lett. 2011, 24, 559–565. [Google Scholar] [CrossRef] [Green Version]
Doha, E.H.; Bhrawy, A.H.; Saker, M.A. On the Derivatives of Bernstein Polynomials: An application for the solution of higher even-order differential equations. Bound. Value Probl. 2011, 16, 2011. [Google Scholar] [CrossRef] [Green Version]
Doha, E.H.; Bhrawy, A.H.; Ezz-Eldien, S.S. Efficient Chebyshev spectral methods for solving multi-term fractional orders differential equations. Appl. Math. Model. 2011, 35, 5662–5672. [Google Scholar] [CrossRef]
Doha, E.H.; Bhrawy, A.H.; Ezz-Eldien, S.S. A Chebyshev spectral method based on operational matrix for initial and boundary value problems of fractional order. Comput. Math. Appl. 2011, 62, 2364–2373. [Google Scholar] [CrossRef] [Green Version]
Doha, E.H.; Bhrawy, A.H.; Hafez, R.M. A Jacobi dual-Petrov-Galerkin method for solving some odd-order ordinary differential equations. Abstr. Appl. Anal. 2011, 2011, 947230. [Google Scholar] [CrossRef] [Green Version]
Khalil, H.; Khan, R.A. A new method based on Legendre polynomials for solutions of the fractional two-dimensional heat conduction equation. Comput. Math. Appl. 2014, 67, 1938–1953. [Google Scholar] [CrossRef]
Khalil, H.; Khan, R.A. New operational matrix of integration and coupled system of fredholm integral equations. Chin. J. Math. 2014, 2014, 146013. [Google Scholar] [CrossRef]
Khalil, H.; Khan, R.A. A new method based on legender polynomials for solution of system of fractional order partial differential equation. Int. J. Comput. Math. 2014, 91, 2554–2567. [Google Scholar] [CrossRef]
Khalil, H.; Khan, R.A. The use of Jacobi polynomials in the numerical solution of coupled system of fractional differential equations. Int. J. Comput. Math. 2014, 92, 1452–1472. [Google Scholar] [CrossRef]
Khalil, H.; Khan, R.A. New Operational Matrix For Shifted Legendre Polynomials and Fractional Differential Equations with Variable Coefficients. J. Math. 2020, 47, 81–103. [Google Scholar]
Saadatmandi, A.; Deghan, M. A new operational matrix for solving fractional-order differential equation. Comput. Math. Appl. 2010, 59, 1326–1336. [Google Scholar] [CrossRef] [Green Version]
Almomani, R.; Almefleh, H. On Heat Conduction Problem with Integral Boundary Condition. J. Em. T. Eng. Sci. 2012, 3, 977–979. [Google Scholar]
Ma, R. A survey on nonlocal boundary value problems. Appl. Math. E-Notes 2007, 7, 257–279. [Google Scholar]
Eziani, D.G.; Shvili, G.A. Investigation of the nonlocal initial boundary value problems for some hyperbolic equations. Hiroshima Math. J. 2001, 31, 345–366. [Google Scholar]
Berikelashvili, G.; Khomeriki, N. On a numerical solution of one nonlocal boundary-value problem with mixed Dirichlet–Neumann conditions. Lith. Math. J. 2013, 53, 367–380. [Google Scholar] [CrossRef]
Sajavicius, S. Radial basis function method for a multidimensional linear elliptic equation with nonlocal boundary conditions. Comput. Math. Appl. 2014, 67, 1407–1420. [Google Scholar] [CrossRef]
Vaqueroa, J.M.; Aguiar, J. On the numerical solution of the heat conduction equations subject to nonlocal conditions. Appl. Numer. Math. 2009, 59, 2507–2514. [Google Scholar] [CrossRef]
Siddique, M. Smoothing of cranknicolson scheme for the two-dimensional diffusion with an integral condition. Appl. Math. Comput. 2009, 214, 512–522. [Google Scholar]
Yang, A.M.; Zhang, Y.Z.; Cattani, C.; Xie, G.N.; Rashidi, M.M.; Zhou, Y.J.; Yang, X.J. Application of Local Fractional Series Expansion Method to Solve Klein-Gordon Equations on Cantor Sets. Abstr. Appl. Anal. 2014, 2014, 372741. [Google Scholar] [CrossRef]
Kumar, S.; Kumar, D.; Abbasbandy, S.; Rashidi, M.M. Analytical solution of fractional Navier–Stokes equation by using modified Laplace decomposition method. Ain Shams Eng. J. 2014, 5, 569–574. [Google Scholar] [CrossRef] [Green Version]
Liu, Y. Numerical solution of the heat equation with nonlocal boundary condition. J. Comput. Appl. Math. 1999, 110, 115–127. [Google Scholar] [CrossRef] [Green Version]
Ang, W. A method of solution for the one-dimensional heat equation subject to nonlocal condition. Southeast Asian Bull. Math. 2002, 26, 185–191. [Google Scholar] [CrossRef]
Dehghan, M. The one-dimensional heat equation subject to a boundary integral specification. Chaos Solitons Fractals 2007, 32, 661–675. [Google Scholar] [CrossRef]
Noye, K.H.B.J. Explicit two-level finite difference methods for the two-dimensional diffusion equation. Int. J. Comput. Math. 1992, 42, 223–236. [Google Scholar] [CrossRef]
Avalishvili, G.; Avalishvili, M.; Gordeziani, D. On integral nonlocal boundary value problems for some partial differential equations. Bull. Georgian Natl. Acad. Sci. 2011, 5, 31–37. [Google Scholar]
Petras, I. A note on the fractional-order Chuas system. Chaos Solitons Fractals 2008, 38, 140–147. [Google Scholar] [CrossRef]
Podlubny, I. Fractional Differential Equations. In Mathematics in Science and Engineering; Academic Press Inc.: San Diego, CA, USA, 1999; Volume 198. [Google Scholar]
Samko, S.G.; Kilbas, A.A.; Marichev, O.I. Fractional Integrals and Derivatives: Theory and Applications; Gordon and Breach: Yverdon, Switzerland, 1993. [Google Scholar]
Ma, R. Multiple positive solutions for nonlinear m-point boundary value problems. Appl. Math. Comput. 2004, 148, 249–262. [Google Scholar] [CrossRef]
Ahmad, B. Approximation of solutions of the forced Duffing equation with m-point boundary conditions. Commun. Appl. Anal. 2009, 13, 11–20. [Google Scholar]
Agarwal, R.P.; Chow, Y.M. Iterative methods for a fourth order boundary value problem. J. Comput. Appl. Math. 1984, 10, 203–217. [Google Scholar] [CrossRef] [Green Version]
AkyuzDascioglu, A.; Isler, N. Bernstein collocation method for solving nonlinear differential equations. Math. Comput. Appl. 2013, 18, 293–300. [Google Scholar]
Bellman, R.E.; Kalaba, R.E. Quasilinearization and Non-Linear Boundary Value Problems; Elsevier: New York, NY, USA, 1965. [Google Scholar]
Charles, A.; Baird, J. Modified quasilinearization technique for the solution of boundary-value problems for ordinary differential equations. J. Optim. Theory Appl. 1969, 3, 227–242. [Google Scholar]
Mandelzweig, V.B.; Tabakin, F. Quasilinearization approach to nonlinear problems in physics with application to nonlinear ODEs. Comput. Phys. Commun. 2001, 141, 268–281. [Google Scholar] [CrossRef] [Green Version]
Gupta, C.P. Solvability of a three-point nonlinear boundary value problem for a second order ordinary differential equation. J. Math. Anal. Appl. 1992, 168, 540–551. [Google Scholar] [CrossRef] [Green Version]
Saeed, U.; ur Rehman, M. Wavelet-Galerkin Quasilinearization Method for Nonlinear Boundary Value Problems. Abstr. Appl. Ana. 2014, 2014, 868934. [Google Scholar] [CrossRef] [Green Version]
Stanley, E.L. Quasilinearization and Invariant Imbedding; Academic Press: New York, NY, USA, 1968. [Google Scholar]

Figure 1. (a) Comparison of exact and approximate solution at different scale levels of Test Problem 1. (b) Absolute difference in the exact and approximate solutions at different scale levels of Test Problem 1.

Figure 2. Errornorm for different values of

σ^{1}

and

σ^{2}

in Test Problem 2.

Figure 2. Errornorm for different values of

σ^{1}

and

σ^{2}

in Test Problem 2.

Figure 3. Errornorm for different values of

σ^{1}

and

σ^{2}

in Test Problem 3.

Figure 3. Errornorm for different values of

σ^{1}

and

σ^{2}

in Test Problem 3.

Figure 4. Error norm on boundary for different values of

σ^{1}

and

σ^{2}

in Test Problem 4.

Figure 4. Error norm on boundary for different values of

σ^{1}

and

σ^{2}

in Test Problem 4.

Figure 5. Error norm of Test Problem 4.

Figure 6. Error norm of Test Problem 5.

Figure 7. Error norm of Test Problem 6.

Table 1. Error norms of Test Problem 1 at scale level

M = 5 : 13

Table 1. Error norms of Test Problem 1 at scale level

M = 5 : 13

M	$\| \| E_{u} {\| \|}_{2}$	$\| \| E_{u} {\| \|}_{\max}$	$\| \| E_{u} {\| \|}_{b}$
5	$0.1314 \times 10^{- 1}$	$0.6431 \times 10^{- 1}$	$8.0032 \times 10^{- 17}$
6	$1.1355 \times 10^{- 2}$	$9.3415 \times 10^{- 1}$	$0.4367 \times 10^{- 16}$
7	$2.3631 \times 10^{- 3}$	$8.3640 \times 10^{- 3}$	$3.0669 \times 10^{- 16}$
8	$3.9474 \times 10^{- 5}$	$3.6641 \times 10^{- 4}$	$9.2291 \times 10^{- 17}$
9	$1.5963 \times 10^{- 5}$	$0.5933 \times 10^{- 5}$	$8.3432 \times 10^{- 17}$
10	$1.6057 \times 10^{- 6}$	$0.5364 \times 10^{- 6}$	$7.7554 \times 10^{- 17}$
11	$1.5823 \times 10^{- 7}$	$1.0592 \times 10^{- 7}$	$6.9857 \times 10^{- 17}$
12	$7.4855 \times 10^{- 9}$	$9.0092 \times 10^{- 8}$	$5.4743 \times 10^{- 17}$
13	$4.5147 \times 10^{- 10}$	$2.4318 \times 10^{- 9}$	$2.4362 \times 10^{- 17}$

Table 2. Error norms of Test Problem 2 at scale level

M = 5 : 13

Table 2. Error norms of Test Problem 2 at scale level

M = 5 : 13

M	$\| \| E_{u} {\| \|}_{2}$	$\| \| E_{v} {\| \|}_{2}$	$\| \| E_{u} {\| \|}_{\max}$	$\| \| E_{v} {\| \|}_{\max}$	$\| \| E_{u} {\| \|}_{b}$	$\| \| E_{v} {\| \|}_{b}$
5	$1.1044 \times 10^{- 1}$	$4.3598 \times 10^{- 2}$	$3.6499 \times 10^{- 1}$	$6.5235 \times 10^{- 2}$	$7.9467 \times 10^{- 17}$	$6.2402 \times 10^{- 16}$
6	$1.9005 \times 10^{- 2}$	$9.3636 \times 10^{- 3}$	$1.6085 \times 10^{- 1}$	$1.6796 \times 10^{- 2}$	$9.7023 \times 10^{- 16}$	$6.2370 \times 10^{- 16}$
7	$2.3621 \times 10^{- 3}$	$7.6605 \times 10^{- 4}$	$7.5980 \times 10^{- 3}$	$1.6380 \times 10^{- 3}$	$9.7237 \times 10^{- 17}$	$6.2396 \times 10^{- 16}$
8	$3.9834 \times 10^{- 5}$	$4.9422 \times 10^{- 5}$	$1.0331 \times 10^{- 4}$	$1.0113 \times 10^{- 4}$	$9.7191 \times 10^{- 17}$	$6.2400 \times 10^{- 16}$
9	$1.4293 \times 10^{- 5}$	$5.3436 \times 10^{- 6}$	$4.5571 \times 10^{- 5}$	$8.6999 \times 10^{- 6}$	$9.7195 \times 10^{- 17}$	$6.2400 \times 10^{- 16}$
10	$1.9227 \times 10^{- 6}$	$5.0884 \times 10^{- 7}$	$6.1295 \times 10^{- 6}$	$7.3992 \times 10^{- 7}$	$9.7194 \times 10^{- 17}$	$6.2399 \times 10^{- 16}$
11	$1.2515 \times 10^{- 7}$	$3.1891 \times 10^{- 8}$	$4.0112 \times 10^{- 7}$	$6.2640 \times 10^{- 8}$	$9.7194 \times 10^{- 17}$	$6.2399 \times 10^{- 16}$
12	$7.1330 \times 10^{- 9}$	$1.9583 \times 10^{- 9}$	$2.2596 \times 10^{- 8}$	$3.6441 \times 10^{- 9}$	$9.7194 \times 10^{- 17}$	$6.2399 \times 10^{- 16}$
13	$4.5177 \times 10^{- 10}$	$1.3237 \times 10^{- 10}$	$1.4501 \times 10^{- 9}$	$2.2449 \times 10^{- 10}$	$9.7194 \times 10^{- 17}$	$6.2399 \times 10^{- 16}$

Table 3. Error norms of Test Problem 3 at scale level

M = 5 : 13

Table 3. Error norms of Test Problem 3 at scale level

M = 5 : 13

N	$\| \| E_{u} {\| \|}_{2}$	$\| \| E_{v} {\| \|}_{2}$	$\| \| E_{u} {\| \|}_{\max}$	$\| \| E_{v} {\| \|}_{\max}$	$\| \| E_{u} {\| \|}_{b}$	$\| \| E_{v} {\| \|}_{b}$
5	$6.5505 \times 10^{- 3}$	$5.1204 \times 10^{- 4}$	$1.5519 \times 10^{- 2}$	$8.6002 \times 10^{- 4}$	$7.0670 \times 10^{- 17}$	$1.8084 \times 10^{- 18}$
6	$3.2583 \times 10^{- 4}$	$2.8746 \times 10^{- 5}$	$7.7139 \times 10^{- 4}$	$4.7281 \times 10^{- 5}$	$7.0336 \times 10^{- 17}$	$1.6297 \times 10^{- 18}$
7	$4.6407 \times 10^{- 6}$	$4.3808 \times 10^{- 7}$	$1.0949 \times 10^{- 5}$	$7.1483 \times 10^{- 7}$	$7.0320 \times 10^{- 17}$	$1.6203 \times 10^{- 18}$
8	$4.6870 \times 10^{- 7}$	$4.6599 \times 10^{- 8}$	$1.1065 \times 10^{- 6}$	$7.8519 \times 10^{- 8}$	$7.0320 \times 10^{- 17}$	$1.6201 \times 10^{- 18}$
9	$4.3761 \times 10^{- 8}$	$3.4302 \times 10^{- 9}$	$1.0365 \times 10^{- 7}$	$5.7780 \times 10^{- 9}$	$7.0320 \times 10^{- 17}$	$1.6201 \times 10^{- 18}$
10	$1.3030 \times 10^{- 9}$	$1.1448 \times 10^{- 10}$	$3.0748 \times 10^{- 9}$	$1.8688 \times 10^{- 10}$	$7.0320 \times 10^{- 17}$	$1.6201 \times 10^{- 18}$
11	$1.2252 \times 10^{- 11}$	$1.1526 \times 10^{- 12}$	$2.8927 \times 10^{- 11}$	$1.8785 \times 10^{- 12}$	$7.0320 \times 10^{- 17}$	$1.6201 \times 10^{- 18}$
12	$7.2561 \times 10^{- 13}$	$7.2734 \times 10^{- 14}$	$1.7118 \times 10^{- 12}$	$1.1828 \times 10^{- 13}$	$7.0320 \times 10^{- 17}$	$1.6201 \times 10^{- 18}$
13	$4.6262 \times 10^{- 14}$	$3.6748 \times 10^{- 15}$	$1.0959 \times 10^{- 13}$	$6.1876 \times 10^{- 15}$	$7.0320 \times 10^{- 17}$	$1.6201 \times 10^{- 18}$

Table 4. Error norms of Test Problem 4 at scale level

M = 5 : 13

Table 4. Error norms of Test Problem 4 at scale level

M = 5 : 13

N	$\| \| E_{u} {\| \|}_{2}$	$\| \| E_{v} {\| \|}_{2}$	$\| \| E_{u} {\| \|}_{\max}$	$\| \| E_{v} {\| \|}_{\max}$	$\| \| E_{u} {\| \|}_{b}$	$\| \| E_{v} {\| \|}_{b}$
5	$0.0108$	$0.0036249$	$0.030961$	$0.0095518$	$- 3.5187 \times 10^{- 17}$	$3.7849 \times 10^{- 15}$
6	$0.0048582$	$0.0020882$	$0.01299$	$0.0075678$	$3.4932 \times 10^{- 17}$	$3.78 \times 10^{- 15}$
7	$0.00055459$	$0.00024661$	$0.001525$	$0.00089155$	$3.4845 \times 10^{- 17}$	$3.7773 \times 10^{- 15}$
8	9.6413 × $10^{- 16}$	$7.5388 \times 10^{- 16}$	$1.8345 \times 10^{- 15}$	$1.2652 \times 10^{- 15}$	$3.486 \times 10^{- 17}$	$3.7777 \times 10^{- 15}$
9	$6.8784 \times 10^{- 16}$	$1.4702 \times 10^{- 15}$	$1.3601 \times 10^{- 15}$	$2.8364 \times 10^{- 15}$	$3.4861 \times 10^{- 17}$	$3.7777 \times 10^{- 15}$
10	$4.5115 \times 10^{- 15}$	$5.6863 \times 10^{- 16}$	$1.302 \times 10^{- 14}$	$8.6167 \times 10^{- 16}$	$3.4861 \times 10^{- 17}$	$3.7777 \times 10^{- 15}$
11	$8.4757 \times 10^{- 17}$	$2.7414 \times 10^{- 16}$	$1.6027 \times 10^{- 16}$	$4.7286 \times 10^{- 16}$	$3.4861 \times 10^{- 17}$	$3.7777 \times 10^{- 15}$
12	$1.5906 \times 10^{- 14}$	$3.1092 \times 10^{- 15}$	$4.6821 \times 10^{- 14}$	$6.3783 \times 10^{- 15}$	$3.4861 \times 10^{- 17}$	$3.7777 \times 10^{- 15}$
13	$4.926 \times 10^{- 15}$	$1.0866 \times 10 - 15$	$1.4378 \times 10^{- 14}$	$2.0953 \times 10^{- 15}$	$3.4861 \times 10^{- 17}$	$3.7777 \times 10^{- 15}$

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Khalil, H.; Khalil, M.; Hashim, I.; Agarwal, P. Extension of Operational Matrix Technique for the Solution of Nonlinear System of Caputo Fractional Differential Equations Subjected to Integral Type Boundary Constrains. Entropy 2021, 23, 1154. https://doi.org/10.3390/e23091154

AMA Style

Khalil H, Khalil M, Hashim I, Agarwal P. Extension of Operational Matrix Technique for the Solution of Nonlinear System of Caputo Fractional Differential Equations Subjected to Integral Type Boundary Constrains. Entropy. 2021; 23(9):1154. https://doi.org/10.3390/e23091154

Chicago/Turabian Style

Khalil, Hammad, Murad Khalil, Ishak Hashim, and Praveen Agarwal. 2021. "Extension of Operational Matrix Technique for the Solution of Nonlinear System of Caputo Fractional Differential Equations Subjected to Integral Type Boundary Constrains" Entropy 23, no. 9: 1154. https://doi.org/10.3390/e23091154

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Extension of Operational Matrix Technique for the Solution of Nonlinear System of Caputo Fractional Differential Equations Subjected to Integral Type Boundary Constrains

Abstract

1. Introduction

2. Preliminaries

The Shifted Legendre Polynomials (LP)

3. Operational Matrices (OP)

4. Application of Operational Matrices

4.1. Linear FDEs with Variable Coefficients

4.2. Nonlinear FDEs

5. Error Bound of the Approximate Solution and Convergence

5.1. Error Bound for Single Differential Equation

5.2. Error Bound for Coupled System of Fractional Differential Equations

6. Test Problems

7. Results and Discussion

8. Conclusions and Future Work

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI